Metodos numericos aplicados a engenharia de potencia

1 UNIVERSIDADE FEDERAL DO MARANHÃO PÓS-GRADUAÇÃO EM ENGENHARIA DE ELETRICIDADE MÉTODOS NUMÉRICOS APLICADOS À ENGENHARIA DE SISTEMAS DE ENERGIA ELÉTRICA PROF. JOSÉ EDUARDO O. PESSANHA 2014 2 CONTEÚDO PARTE 1 Introdução Solução de Equações Lineares Associadas ao Problema de Fluxo de Carga    O Subproblema Linear de Fluxo de Carga Características Estruturais da Matriz Jacobiana Propriedades Espectrais da Matriz Jacobiana de Fluxo de Carga Métodos Diretos     O método de Eliminação de Gauss-Jordan Eliminação Gaussiana com Substituição Regressiva Estratégias de Pivotamento Fatoração Triangular LU Métodos Iterativos do Subespaço de Krylov     A Escolha do Método Iterativo Características dos Principais Métodos Iterativos Vantagens e Desvantagens dos Principais Métodos Iterativos Avaliação do Desempenho dos Principais Métodos Iterativos Sem Estratégias O Método do Resíduo Mínimo Generalizado – GMRES    Exemplo Ilustrativo do GMRES (Sem Estratégias) Estratégia de Reinicialização Para O Método GMRES Exemplo Ilustrativo – GMRES(m) Reinicializado Estratégias de Reordenamento da Matriz Jacobiana   Características Principais das Técnicas de Reordenamento MD e RCM Avaliação do Desempenho do Método GMRES Puro Usando Reordenamento da Matriz de Coeficientes 3 PARTE II Estratégias de Pré-condicionamento   O Pré-condicionador Formas de Pré-condicionamento Efeitos do Pré-condicionamento no Subespaço de Krylov do Método GMRES  Efeitos da Estratégia de Pré-condicionamento nas Propriedades Espectrais da Matriz Jacobiana Pré-condicionadores de Fatoração Incompleta     Regras de Preenchimento em Pré-condicionadores de Fatoração Incompleta Pré-condicionador ILU(0) Pré-condicionador ILU(k) Pré-condicionador ILUT(,ρ) Pré-condicionadores Aplicados a Problemas de Sistemas Elétricos de Potência  Características dos Principais Pré-condicionadores Avaliação do Desempenho do Método GMRES Reordenamento da Matriz Jacobiana Evolução do Desempenho do GMRES 4 1 Introdução Muitos componentes dos sistemas de energia elétrica, dependendo do fenômeno investigado, são matematicamente modelados por sistemas de equações algébricas e diferenciais (dinâmicas rápidas e lentas associadas à estabilidade de tensão, angular e de frequência), equações diferenciais ordinárias (transitórios eletromagnéticos), ou sistemas puramente algébricos lineares (regime-permanente). Essencialmente, quando se investiga um sistema de energia elétrica real de grandes dimensões, tem-se associado à primeira e a terceira categoria, sistemas algébricos lineares também de grande porte, e em geral, a matriz de coeficientes resultante é esparsa com uma estrutura apropriada de elementos não nulos (zeros), e a solução pode ser obtida por um método direto. Entretanto, métodos diretos podem ser bem dispendiosos em termos de espaço de memória e tempo de CPU, visto que durante o processo de solução podem aparecer elementos não nulos, em posições anteriormente ocupadas por elementos nulos, mesmo com estratégias de reordenamento para reduzir estas entradas. Nestes casos, métodos iterativos são mais atrativos em comparação aos métodos diretos, mesmo se a matriz for densa; seus requisitos de espaço de memória e tempo de CPU são geralmente menores da ordem de magnitudes. Ao invés de solucionar o sistema original do tipo (1.1)Erro! Fonte de referência não encontrada., métodos iterativos solucionam um sistema equivalente menos complexo, ou seja, o sistema original é substituído por outro cuja solução é mais fácil de obter. Além disso, métodos iterativos não requerem o armazenamento de matrizes, apenas a habilidade de realizar multiplicações matriz-vetor. Ax =b (1.1) 5 2 Problema de Fluxo de Carga 2.1 Introdução A solução do problema de fluxo de carga em sistemas elétricos de potência, assim como em diversas áreas da engenharia e da ciência, requer a solução de sistemas de equações lineares esparsos (1.1) e, em muitos casos, mal condicionados. Nesta relação, A é uma matriz quadrada, conhecida como matriz de coeficientes; b é um vetor coluna de valores conhecidos e, x é um vetor de incógnitas. Já em (2.1), x* é a solução calculada com uma determinada precisão numérica. Para a solução de (1.1), é possível classificar os métodos em duas categorias: diretos e iterativos. Um método direto, por exemplo, decompõe a matriz de coeficientes em um produto de dois fatores triangulares, L e U, sendo posteriormente resolvidos mediante um processo conhecido como substituição direta-inversa, resultando na solução do sistema original. Este é conhecido como decomposição LU. x* = A -1  b (2.1) No caso dos métodos iterativos, a ideia central é que, a partir de uma solução inicial x 0, seja gerada uma sequência de soluções {xk}, isto é, uma solução aproximada x k em cada iteração k, tal que depois de um determinado número finito de iterações resulte na solução x*. Inicialmente, a aplicação dos primeiros métodos iterativos (estacionários) em estudos envolvendo sistemas de energia elétrica, tais como Jacobi, Gauss-Seidel e sobrerelaxação sucessiva (SOR), foi tratada com certa cautela pelos pesquisadores sob o argumento de baixa confiabilidade e lentidão. Atualmente, tanta cautela não se justifica, devido ás características dos métodos iterativos mais modernos, baseados na teoria de Subespaços de Krylov, como por exemplo, o Resíduo Mínimo Generalizado (GMRES). Portanto, métodos numéricos deste tipo podem ser vistos como eficientes alternativas de solução, já que exploram muito bem a estrutura esparsa da matriz de coeficientes, visando sempre economizar memória e operações de pontos flutuantes, o que se refletirá no tempo de execução. Ainda, uma vez que a robustez do método está associada ao grau de influência no resultado dos erros acumulados nas diversas operações efetuadas ao longo do processo de solução, os métodos iterativos neste aspecto apresentam uma vantagem sobre arredondamento. os métodos diretos, visto que acumulam menos erros de 6 A escolha de um método iterativo apropriado para a solução do subproblema linear de fluxo de carga é o objetivo principal aqui, sendo uma das tarefas mais importantes no desenvolvimento de um solucionador iterativo, dado que esta definição afeta decisões futuras associadas com a adequação de estratégias numéricas no solucionador, tais como, pré-condicionamento e reordenamento. Não menos importante é determinar o procedimento adequado para a seleção do método iterativo. Considerando-se que o desempenho dos métodos varia em função das características do problema, estudam-se aqui os métodos cujas propriedades sejam mais vantajosas para solucionar o subproblema linear de fluxo de carga e que explore de forma eficiente as características da matriz de coeficientes. Portanto, são estudadas as principais características dos métodos iterativos, destacandose suas vantagens e desvantagens. Adicionalmente, são estabelecidas típicas matrizes de coeficientes associadas a diferentes sistemas elétricos e níveis de carga, para avaliação e comparação destes métodos. A análise dos resultados, em termos de robustez e eficiência computacional, é apresentada aqui, visando validar a escolha do método. Entre os métodos considerados, é dada maior ênfase ao Resíduo Mínimo Generalizado (GMRES) por apresentar características mais convenientes. Faz-se uma abordagem teórica dos conceitos associados a este método, explicando-se de forma ilustrativa os processos envolvidos no algoritmo do GMRES. Adicionalmente, são apresentadas e avaliadas as estratégias de reordenamento e reinicialização, muito usadas para melhorar o desempenho dos solucionadores iterativos. Nos próximos itens são apresentados conceitos associados ao problema de interesse, as características numéricas, estruturais e as propriedades espectrais de suas matrizes de coeficientes. 2.2 O Subproblema Linear de Fluxo de Carga O cálculo de fluxo de carga em uma rede de energia elétrica consiste essencialmente na determinação do estado da rede, da distribuição dos fluxos e de algumas outras grandezas de interesse. Neste tipo de problema, a modelagem do sistema é estática, significando que as variações com o tempo são suficientemente lentas para não incluir os efeitos transitórios. A solução do problema de fluxo de carga requer a solução de um sistema de equações algébricas não-lineares, deduzidas pela aplicação das leis de Kirchhoff, resultando em duas equações por cada barra k da rede elétrica, dadas por: 7 Pk  Pkesp  Vk  Qk  Qkesp  Vk   Vm Gkm  cos(k  m )  Bkm  sen(k  m ) mK  mK Vm  Gkm  sen(k  m )  Bkm  cos(k  m )  (2.2) Para k=1,...,n; sendo n o número de barras totais da rede elétrica; k e m são os subíndices indicativos de que a variável está associada às barra k e m, respectivamente. Por exemplo, ΔPk é o erro de potência ativa na barra k (barra do tipo PQ ou PV), ΔQ k é o erro de potência reativa na barra k (barra do tipo PQ); P kesp e Qkesp são as potências ativa e reativa especificadas na barra k, [V 1, θ1, ..., Vk, θk, ..., Vm, θm , ..., Vn, θn] são as incógnitas do sistema de equações não-lineares associadas com os módulos e ângulos das tensões nas barras da rede, finalmente, G km e Bkm são a condutância e susceptância, respectivamente, associadas com as posições k e m da matriz de admitância nodal. A variável K representa o conjunto de barras conectadas fisicamente à barra k, incluindo-se ela mesma. Este sistema de equações algébricas não-lineares é solucionado por meio de linearizações sucessivas (através da Série de Taylor), usando-se o método NewtonRaphson. O algoritmo do método Newton-Raphson pode ser expresso na iteração i, como a solução de um subproblema linear (sistema linearizado) apresentado em (2.3), no qual a matriz de coeficientes (ou Jacobiana J) é recalculada a cada iteração, seguido da atualização dos valores dos módulos e ângulos das tensões das barras, como em (2.4). A matriz Jacobiana é definida como a matriz de derivadas parciais dos erros de potência ativa e reativa, calculada da forma apresentada em (2.5), onde H, N, M e L são as submatrizes de derivadas parciais.  Δθi   ΔPi  J  =   ΔVi   ΔQi  (2.3)  θi+1   θi   Δθi       Vi1   Vi   ΔVi  (2.4)  ΔP  θ J   ΔQ  θ ΔP   H N V      ΔQ  M L  V  (2.5) 8 O sistema de equações lineares apresentado em (2.3) representa o problema principal que será resolvido em cada iteração de fluxo de carga pelo solucionador proposto neste trabalho. Quando usados métodos iterativos para solução deste tipo de equações, é conveniente caracterizar a matriz de coeficientes do problema que está sendo solucionado a fim de selecionar apropriadamente o método e suas estratégias numéricas, principalmente a de pré-condicionamento. A seguir são apresentadas as principais características estruturais, numéricas e espectrais da matriz Jacobiana de fluxo de carga. 2.2.1 Características Estruturais da Matriz Jacobiana Conforme apresentado em (2.5), a matriz Jacobiana é formada por quatro submatrizes H, N, M e L estruturalmente idênticas e simétricas; consequentemente, a matriz Jacobiana também é estruturalmente simétrica, embora seja numericamente assimétrica. Outra característica relevante da matriz Jacobiana é sua esparsidade. Esta característica é inerente às redes elétricas dado que cada barra tem apenas alguns poucos ramos conectados. Visto que o número médio de ramos conectados às barras é o mesmo, independentemente do tamanho do sistema, então, quanto maior for o sistema, maior será o número de elementos nulos na matriz Jacobiana, tornando-a cada vez mais esparsa. Nas Figuras 2.1 e 2.2 são apresentadas estruturas de típicas matrizes Jacobianas para diferentes sistemas elétricos, observando-se as características estruturais mencionadas acima. 9 a) Sistema IEEE30-barras b) Sistema IEEE118-barras c) Sistema IEEE145-barras d) Sistema IEEE162-barras Figura 2.1 – Estruturas típicas de matrizes Jacobianas de sistemas hipotéticos. 10 a) Sistema IEEE300-barras b) Sistema Norte-Nordeste de 274 barras c) Sistema brasileiro de 2.256 barras d) Sistema brasileiro de 3.513 barras Figura 2.2 – Estruturas típicas de matrizes Jacobianas de sistemas hipotético e reais. Em todos os gráficos das figuras apresentadas acima, pode-se notar as estruturas similares, esparsas e simétricas das submatrizes H, N, M e L. Dados numéricos destas matrizes, como por exemplo, dimensão, número de elementos não-nulos, número médio e máximo de elementos por linha da matriz Jacobiana e grau de esparsidade são apresentados na Tabela 2.1. Observa-se que, conforme aumenta a dimensão do sistema elétrico, aumentam também as dimensões das matrizes Jacobianas, o número de elementos não-nulos e o grau de esparsidade. A esparsidade aumenta porque o número de elementos por linha da matriz Jacobiana (médio e máximo) não aumenta significativamente. O número médio de elementos por linha é calculado como o número total de elementos dividido pelo número de linhas da matriz Jacobiana. Já o número máximo de elementos por linha corresponde à linha com maior número de elementos, estando associada com a barra do sistema elétrico com maior número de ramos conectados. Por exemplo, no caso da configuração do sistema 11 interligado nacional de 3.513 barras, a barra com maior número de conexões é a barra 180, identificada como JACAREP—138. Esta barra está conectada a 18 barras do sistema. Tabela 2.1 – Dados numéricos de matrizes Jacobianas de sistemas elétricos de potência. Sistema Dimensão da Matriz Jacobiana Número de Elementos Não-nulos Número de Elementos por Linha Médio Máximo Grau de Esparsidade IEEE30-barras 60 447 7 16 87,58% IEEE118-barras 236 1.900 8 20 96,59% IEEE145-barras 290 3.947 14 42 95,31% IEEE162-barras 324 2.884 9 20 97,25% IEEE300-barras 600 4.472 7 24 98,76% 548 3.316 6 27 98,90% 4.512 30.918 7 32 99,85% 7.026 44.300 6 38 99,91% Norte-Nordeste de 274 barras Brasil de 2.256 barras Brasil de 3.513 barras 2.2.2 Propriedades Espectrais da Matriz Jacobiana de Fluxo de Carga As propriedades espectrais são usadas para identificar maus condicionamentos e singularidades em sistemas de equações lineares. Estas propriedades são inerentes às matrizes de coeficientes e variam em função do problema que está sendo solucionado. Neste trabalho, estas propriedades são usadas para caracterizar típicas matrizes Jacobianas associadas ao problema de fluxo de carga para sistemas hipotéticos do IEEE e configurações dos sistemas norte-nordeste e brasileiro. As propriedades espectrais observadas são os autovalores e o número de condicionamento da matriz de coeficientes. O conjunto de todos os autovalores da matriz quadrada A é chamado de espectro de A. Esta propriedade pode ser um indicador de singularidade da matriz de coeficientes, de cenários de difícil convergência para métodos iterativos e de difícil solução para métodos diretos. Por exemplo, quando usado o método Newton-Raphson para solucionar sistemas não-lineares, normalmente uma matriz Jacobiana que apresente autovalores tanto com parte real negativa como também positiva (matriz indefinida), representa um dos piores cenários para um método iterativo, isto é, nestes casos, o solucionador requer o uso de pré-condicionadores para que o método iterativo consiga convergir. 12 Uma matriz é singular se admite um autovalor igual à zero, portanto, autovalores muito próximos de zero podem ser aparentemente indicadores de singularidade. No entanto, os autovalores, algumas vezes, podem ocultar singularidades em matrizes escalonadas, já que seus valores dependem muito do escalonamento da matriz, portanto, nestes casos, recomenda-se também o uso do número de condicionamento, pois esta propriedade independe do escalonamento da matriz. O número de condicionamento, definido em (2.6), segundo a teoria da perturbação para sistemas lineares relaciona as variações que acontecem na solução x* com o tamanho das possíveis perturbações (por exemplo, erros de arredondamento) que podem ocorrer durante o processo de solução, como visto em (2.7). Portanto, se o número de condicionamento é grande, então o problema A∙x = b é dito mal condicionado, isto é, pequenos erros em A ou b causam grandes erros em x. O número de condicionamento de norma-2, usado aqui, indica se a matriz está próxima da singularidade. k(A)  A 2  A x x  k(A) 1 2 b b (2.6) (2.7) Ambas as propriedades são muito usadas em estudos de sistemas elétricos mal condicionados, principalmente para identificar singularidades e também para avaliar as melhorias obtidas quando utilizadas estratégias numéricas como, por exemplo, o précondicionamento e o reordenamento. Nas Figuras 2.3 e 2.4 apresentam-se os gráficos dos autovalores das matrizes Jacobianas associadas aos sistemas IEEE30-barras, IEEE118-barras, IEEE162-barras, IEEE145-barras, IEEE300-barras, Norte-Nordeste de 274 barras e duas configurações do sistema brasileiro de 2.256 e 3.513 barras, respectivamente. Cada sistema possui dois cenários, chamados Caso Base e ponto mais Próximo do Máximo Carregamento (PMC), representados pelos autovalores azuis e vermelhos, respectivamente. Os casos base dos sistemas IEEE30-barras, IEEE118-barras, IEEE162-barras, IEEE145-barras e IEEE300barras correspondem aos pontos padrões de operação aceitos em relatórios e forçastarefas realizados pela IEEE. O ponto PMC foi obtido aumentando-se gradativamente o consumo de potência ativa e reativa em todas as barras de carga, mantendo-se o fator de potência constante utilizando-se o programa de Fluxo de Potência Continuado do ANAREDE. Portanto, o 13 último ponto de operação convergido de fluxo de carga Continuado foi considerado como o ponto mais próximo do máximo carregamento (PMC). Nas Figuras 2.3 e 2.4 observa-se que, à medida que o tamanho dos sistemas elétricos aumenta, os autovalores ocupam uma região maior do plano complexo, aumentando também seus módulos e a magnitude de suas partes reais. Observa-se também que, nos sistemas IEEE300-barras, Norte-Nordeste e brasileiros aparecem autovalores com parte real negativa, isto é, para estes sistemas, as matrizes Jacobianas são indefinidas. Nestes casos, os métodos iterativos podem apresentar convergência deficiente ou até falhar, sendo, portanto, necessário o uso de pré-condicionadores para que o solucionador iterativo consiga convergir. Tabela 2.2, onde são apresentados os autovalores com menor e maior parte real, tanto para o caso base como também para o ponto PMC. a) Sistema IEEE30-barras b) Sistema IEEE118-barras 14 c) Sistema IEEE145-barras d) Sistema IEEE162-barras Figura 2.3 – Autovalores da matriz Jacobiana de sistemas de pequeno porte, no caso base e no ponto mais próximo do máximo carregamento (PMC). Nos gráficos espectrais pode-se observar que, tanto no caso base como no PMC, a maior concentração de autovalores encontra-se próxima ao zero do plano complexo. A aproximação do zero é maior ainda nos sistemas de maiores dimensões indicando que todos os sistemas testados apresentam matrizes Jacobianas próximas da singularidade. A proximidade do zero é maior para as matrizes Jacobianas obtidas no PMC (autovalores em vermelho). No sistema IEEE145-barras, a diferença espectral entre o caso base e o PMC não é nítida, dado que ambos os pontos de operação encontram-se próximos, ou seja, o caso base deste sistema apresenta um nível de carregamento similar ao ponto PMC. a) Sistema IEEE300-barras b) Sistema Norte-Nordeste de 274 barras 15 c) Sistema Brasil de 2.256 barras d) Sistema Brasil de 3.513 barras Figura 2.4 – Autovalores da matriz Jacobiana de sistemas de maiores dimensões no caso base e no ponto mais próximo do máximo carregamento (PMC). Tabela 2.2 – Propriedades espectrais de típicas matrizes Jacobianas de sistemas elétricos de potência. Autovalores no Caso Base Sistema Com a Menor Parte Real Com a Maior Parte Real Autovalores no Ponto Mais Próximo do Máximo Carregamento Com a Menor Parte Real Com a Maior Parte Real IEEE30-barras 0,23309 104,32±28,87i 0,00798 95,70±26,42i IEEE118-barras 0,18924 373,78±169,73i 0,00983 360,60±164,28i IEEE145-barras 0,20514 5.683,93±210,43i 0,00519 5.621,68±202,72i IEEE162-barras 0,11977 2.187,74±216,50i 0,01871 2.122,98±209,50i IEEE300-barras -1,45621 4.456,39±603,11i -1,45015 4.164,81±673,07i -1.873,07495 20.662,45 -1.873,08000 20.703,42 -6.761,09±14,30i 43.191,30±51,67i -6.036,07±12,42i 42.284,53±11,23i Norte-Nordeste de 274 barras Brasil de 2.256 barras Brasil de 3.515 barras -6.317,46±32,81i 307.962,64±30,09i -6.286,24±33,34i 304.360,52±30,86i Com a finalidade de confrontar os resultados anteriores, na Figura 2.5 apresentam-se em escala logarítmica os valores do número de condicionamento das matrizes Jacobiana de todos os sistemas considerados aqui. Observe-se que o número de condicionamento 16 aumenta diretamente com as dimensões dos sistemas, sendo o sistema brasileiro de 3.513 barras o pior em termos de condicionamento, seguido do sistema brasileiro de 2.256 barras. Portanto, estes dois sistemas estão mais susceptíveis à problemas de convergência quando forem utilizados métodos iterativos. Nota-se também que o número de condicionamento é maior no ponto mais próximo do máximo carregamento, já que este ponto de operação está mais próximo da singularidade. Ambas as propriedades comprovam que, o aumento na dimensão dos sistemas elétricos bem como nos seus níveis de carga dificultam, de fato, a solução do subproblema linear de fluxo de carga e consequentemente do problema global. Baseado no exposto, o desenvolvimento de solucionadores e a adaptação de métodos iterativos robustos e eficientes para superar dificuldades associadas à solução do subproblema linear é um desafio interessante e, por esta razão, se tornou o tema principal de investigação neste trabalho. Antes de apresentar a caracterização dos principais métodos iterativos, suas vantagens e desvantagens, é interessante entender o significado do Subespaço de Krylov, como é gerado e como os métodos baseados nesta teoria matemática procuram soluções aproximadas dentro do subespaço. A seguir apresentam-se todos estes conceitos de forma ilustrativa e didática. Figura 2.5 – Número de condicionamento para o caso base e para o ponto próximo do máximo carregamento (PMC). 2.3 Métodos Diretos O método de Eliminação de Gauss-Jordan Pode ser aplicado para resolver sistemas de equações lineares da forma 2.1, para um ou mais vetores a esquerda [b], como também determinar a inversa da matriz A (matriz cheia) se necessário. Para atividades especificas de inversão da matriz de coeficientes, a 17 técnica Gauss-Jordan apresenta uma aceitável eficiência quando comparada à outras técnicas. Entretanto, existem desvantagens associadas a esse método, sendo estas:   Armazenamento e manipulação simultânea de todos os vetores a esquerda, Caso o cálculo da inversa da matriz A não seja necessário, o tempo computacional do método é reduzido em pelo menos três vezes se comparado aos similares quando aplicados para o mesmo fim, e  Quando a matriz inversa é determinada explicitamente, aplica-se um novo vetor [b], resultando em uma solução diferente da encontrada no processo inicial. A razão está no erro de arredondamento da aritmética computacional associada. Uma vantagem interessante e importante é a confiabilidade do método baseada na estabilidade apresentada quando comparado à outros métodos diretos, aumentando sua robustez se uma estratégia de pivotamento completo for aplicada. A sequência de operações executadas na eliminação Gauss-Jordan constitui-se como uma base de novas alternativas, como é o caso dos métodos de eliminação Gaussiana e decomposição triangular. Eliminação Gaussiana com Substituição Regressiva Para ilustrar os mecanismos algébricos deste método, considere o simples sistema algébrico representado na Equação (2.8a e 2.8b). E1 : x1  x 2  3x 4  4 E2 : 2x1  x 2 - x3  x 4  1 E3 : 3x1 - x 2 - x 3  2x 4  -3 (2.8a) E4 : -x1  2x 2  3x3 - x 4  4 E2 - 2E1   E2  E3 - 3E1   E3  E4  E1   E4  E1 : x1  x 2  E2 : E3 : E4 :  3x 4  4 - x 2 - x3 - 5x 4  -7 3x3  13x 4  13 (2.8b) - 13x 4  -13 O sistema 2.8b está na forma triangular e pode ser facilmente resolvido pelo processo de substituição regressiva. A Equação 2.9 mostra a forma geral para solução através do método de Eliminação Gaussiana com Substituição Regressiva. 18 k     E - aik E   E i  i ak  k  kk   (2.9) Se akk for zero, não significa necessariamente que o sistema não tenha solução, mas que o método deva ser alterado. Para ilustrar a forma de resolver este problema, considere o exemplo a seguir. Ex. Considere o sistema linear E1 : 1x1  1x 2  2x3  x 4  8 E2 : 2x1  2x 2  3x3  3x 4  20 E3 : 1x1  1x 2  1x3  0x 4  2 E4 : 1x1  1x 2  4x3  3x 4  4 Formando a matriz aumentada e operando (E 2-2E1)(E2); (E3-E1)(E3) e (E4-E1)(E4), tem-se: A (2)  1 1 2 1 8  0 0 1 1 4    0 2 1 1 6    0 0 2 4 12 Uma vez que a22(2)=0, chamado de pivô, o procedimento não pode continuar na presente forma, mas a operação (Ei)(Ej) é permitida, então uma procura é feita dos elementos a32(2) e a42(2) pelo primeiro elemento não nulo. Uma vez que a 32(2) 0 a operação (E2)  (E3) é realizada, resultando em: A (2)'  1 1 2 1 8  0 2 1 1 6    0 0 1 1 4    0 0 2 4 12 Uma vez que x2 é eliminado de (E3) e (E4), A(3) será A(2) ' e o cálculo poderá continuar com a operação (E4+2E3)(E4), resultando em: 19 A (4)  1 1 2 1 8  0 2 1 1 6    0 0 1 1 4    0 0 0 2 4  Enfim, a substituição regressiva é aplicada, encontrando a solução x  [2,2,3, 7]t . Estratégias de Pivotamento A troca de linhas é necessária quando um elemento pivô a kk(k) é zero. Esta troca de linhas tem a forma (Ei)(Ep) onde p é o menor inteiro maior que k com a pk(k)=0. De fato, é geralmente necessário realizar trocas de linhas mesmo quando os elementos pivô não forem zero, a fim de se reduzir os erros de arredondamento associados a aritmética de dígitos finitos. Para se evitar este problema, o processo de pivotamento é realizado selecionando-se um elemento grande apq(k) para pivô e trocar as k-ésimas e p-ésimas linhas, seguido das trocas das k-ésimas e p-ésimas colunas, se necessário. A estratégia mais simples é selecionar o elemento na mesma coluna que esteja abaixo da diagonal e que possua o maior valor absoluto. Um exemplo ilustrativo que segue pode ser útil para se entender este procedimento. Ex. Considere o sistema linear: E1:0,003000x1+59,140x 2 =59,17 E2 : 5,291x1- 6,130x2 =46,78 Cuja solução exata é x1 = 10,00 e x2 = 1,000. Para ilustrar as dificuldades do erro de arredondamento, usa-se eliminação Gaussiana para solucionar o sistema considerandose aritmética de quatro dígitos com arredondamento. O primeiro elemento pivô é um número pequeno, a11(1)=0,003000 , e o seu multiplicador associado é: m21= 5,291 =1763,66 0,003000 O resultado considerando-se aritmética de quatro dígitos com arredondamento é um número grande, ou seja, 1764. Fazendo (E2-m21E1)(E2) e o arredondamento apropriado tem-se: 0,003000x1  59,14x 2  59,17 - 104300x 2  -104400 Ao invés dos valores precisos: 20 0,003000 x1+ 59,14 x 2 = 59,17 - 104309,376 x 2 = -104309,376 A disparidade nas amplitudes de m21a13 e a23 introduziu erros de arredondamento, mas ainda não foram propagados. Substituição regressiva leva a: x2  1,001 Que é uma boa aproximação para o valor atual, x2 = 1,000. Entretanto, devido ao pequeno pivô a11, x1 » 59,17 -  59,14 1,001 0,003000 = -10,00 Arruinando desta forma a aproximação ao valor exato, que é x 1 = 10,00. O procedimento   de pivotamento é iniciado encontrando-se: 1 1 1 max a11 , a21 = max  0,003000 , 5,291 = 5,291 = a21 A operação (E2)(E1) é realizada resultando no seguinte sistema: E1 : 0,003000 x1+ 59,14 x 2 = 59,17 E2 : 5,291x1- 6,130 x 2 = 46,78 O multiplicador para este sistema é: 1 a21 m21 = = 0,0005670 1 a11 e a operação (E2-m21E1)(E2) reduz o sistema a: 5,291x1 -6,130x2  46,78 59,14x2  59,14 As soluções de quatro dígitos resultantes da substituição regressiva são os valores corretos x1=10,00 e x2=1,000. Esta técnica é conhecida como pivotamento parcial, existindo outras, mas o objetivo aqui não é discuti-las, mas sim ilustrar o porquê, quando, onde aplicá-las e o resultado final. Tanto as técnicas diretas quando as iterativas fazem uso desse artifício para evitar problemas de precisão. Fatoração Triangular LU Considere a hipótese de se representar a matriz A como um produto de duas matrizes, 21 LU=A (2.10) Onde L é uma matriz triangular inferior, e U uma matriz triangular superior. O sistema linear de interesse pode ser representado na seguinte forma matricial:  11 0  21 22  31 32   41  42 0 0 33  43 0  11 12 13  0   0 22 23 0  0 0 33   44   0 0 0 14   a11   24  a21  34  a31   44  a41 a12 a22 a13 a23 a32 a33 a42 a43 a14   a24  a34   a44  (2.11) Considerando-se a decomposição (Equação 2.10) para resolver o conjunto linear: Ax=(LU)x=L(Ux)=b (2.12) Pode-se observar por inspeção a possibilidade de dividir o sistema linear em dois sistemas triangulares lineares, resultando na seguinte simplificação: Ly=b (2.13) Ux=y (2.14) E depois resolver: A vantagem está na facilidade em se resolver um sistema de equações triangulares, simplificando o processo de solução. Assim, a Equação 2.13 pode ser resolvida através do processo substituição progressiva, representada pelas Equações 2.15 - 2.16. Já o sistema representado pela Equação 2.14 pode ser resolvido pelo processo de substituição regressiva, representada pelas Equações 2.17 - 2.18. É interessante notar que, uma vez que determinada a decomposição LU de A é possível resolver mais rapidamente o sistema algébrico ao variar o vetor a esquerda, sendo esta uma vantagem relevante sobre outros métodos. yi  b1 / 11 (2.15) 22 i1 yi  1/ i,i bi   j1i,j yi    i=2,3,..n xn  yn / n,n (2.16) (2.17) n xi  1/ i,i  yi   ji1i,j y j    i = n-1, n-2,..,1 (2.18) Para matrizes altamente esparsas e de grande porte, o processo de solução não é, a principio, diferente do método geral da fatoração triangular LU; o problema está no preenchimento de elementos não nulos e grandes (preenchimento) durante o processo de solução. Estes esquemas de decomposição (ou de eliminação) devem minimizar o número de elementos indesejados não zero. Estes novos elementos, inicialmente zeros, tornam-se valores representativos durante o processo de solução, e, como consequência, o espaço de armazenagem requerido é bem maior. As n entradas conhecidas da matriz A (n  n) podem ser usadas para determinar parcialmente as m incógnitas de L e o mesmo número em U. Se uma única solução é desejada, entretanto, certo número de condições adicionais nas entradas de L e U são necessárias. Vários métodos podem ser usados para estabelecer essas condições, como o Doolittle’s. 2.4 Métodos Iterativos do Subespaço de Krylov Um subespaço Krylov κk é definido por uma sequência de vetores base, como apresentado em (2.19), no qual A é a matriz de coeficientes e r0 (2.20) o resíduo inicial. κ k (A,r0 ) = subespaço(r0 ,A  r0 ,A 2  r0 ,...,Ak-1  r0 ) (2.19) r0 = b - A  x0 (2.20) Um método do subespaço Krylov para resolver sistemas lineares na forma (2.21) procura uma solução aproximada xk dentro de um subespaço x0+κk, impondo alguma condição que facilite a busca, como por exemplo, a condição da norma mínima do erro, apresentada em (2.21), onde x* é a solução exata do sistema linear. x * -xk 2  = mín x * -x 2 : x  x0 + κk  (2.21) 23 Para ilustrar como um método iterativo em uma iteração k pode usar a condição da norma mínima do erro na busca de uma aproximação xk dentro de um subespaço x0+κk, considera-se um sistema linear de três incógnitas. Neste caso, para fins de simplificação, é assumido k igual a 2, ou seja, o método iterativo está na segunda iteração e a primeira aproximação é zero (x0 = 0). Então, o subespaço onde é procurada a aproximação é dado por: x0 + κ 2 = x0 + subespaço(r0 ,A  r0 ) = subespaço(r0 ,A  r0 ) (2.22) Convenientemente, o subespaço está sendo definido unicamente por dois vetores (r 0 e A·r0) e a sua representação geométrica corresponde ao plano da Figura 2.6. Qualquer vetor x pertencente ao plano pode ser definido como uma combinação linear dos vetores base r0 e A·r0 , da forma: x = c0∙r0 + c1∙A·r0. Figura 2.6 – Representação geométrica do subespaço de Krylov pelos seus vetores base Na Figura 2.7, observa-se que a aproximação xk pode ser qualquer vetor x (em cinza) pertencente ao plano, produzindo um vetor de erro x*-x (linha pontilhada) cujo módulo pode aumentar ou diminuir em função da escolha de x. Entretanto, a condição da norma mínima do erro exige minimizar x*-x e por uma simples inspeção na figura pode-se notar que o erro mínimo x*-x corresponde ao vetor representado pelas linhas tracejadas, sendo perpendicular (ortogonal) ao plano (subespaço) que une o plano com a solução exata x*. Finalmente, a interseção deste vetor com o plano determina a aproximação xk (vetor cinza em linha tracejada) que satisfaz a condição da norma mínima do erro e cujo módulo é igual à projeção ortogonal da solução exata x* sobre o plano. As aproximações obtidas a partir dos métodos do subespaço de Krylov são da forma apresentada em (2.23), onde qk−1 é um polinômio de grau k−1, apresentado em (2.24). 24 Figura 2.7 – Método do subespaço de Krylov usando a condição da norma mínima do erro. A -1  b  xk = x0 + qk-1(A)  r0 qk-1(A) = c0 + c1  A + c 2  A 2 + ... + ck-1  Ak-1 (2.23) (2.24) É prática comum definir a primeira aproximação como um vetor de zeros (x 0 = 0), consequentemente, o resíduo associado (r0) é o vetor b. Neste caso o subespaço de Krylov é da forma (2.25) e a inversa da matriz A pode ser representada apenas pelo polinômio qk−1, como mostrado em (2.26). κ k (A,b) = subespaço(b,A  b,A2  b,...,,Ak-1  b) A -1  b  xk = qk-1(A)  b A -1  qk-1(A) (2.25) (2.26) Os vetores base b, A·b, A2·b,..., Ak-1·b, normalmente são quase linearmente dependentes, o que dificulta a busca da aproximação. Portanto, não é recomendável realizar a busca da aproximação no subespaço sem antes ortonormalizar suas bases, para torná-las mais linearmente independentes, facilitando o processo de busca como realizado, por exemplo, no algoritmo do método iterativo GMRES. Métodos que fornecem melhores aproximações no subespaço Krylov são conhecidos como métodos de projeção de Krylov devido às classes de projeções (ortogonais ou bi-ortogonais) usadas normalmente para reduzir o resíduo. Esses métodos podem ser classificados de acordo com quatro diferentes condições de busca da solução, sendo estas: i. Condição de Ritz-Galerkin. Constrói xk e o resíduo é ortogonal com relação ao mais recente subespaço construído: (b - A∙xk) ┴ КK(A, r0). 25 ii. Condição da norma mínima residual. Identifica x k cuja norma Euclidiana ║b A∙xk║2 é mínima sobre КK(A, r0). iii. Condição de Petrov-Galerkin. Encontra um xk de forma que o resíduo b - A∙xk seja ortogonal a algum subespaço k-dimensional aceitável. iv. Condição da norma mínima do erro. Determina x k em AT∙КK(AT, r0) pelo qual a norma Euclidiana ║x* - xk║2 é mínima (ilustrada na Figura 2.7). Como mostrado na Tabela 2.3, a primeira abordagem, Ritz-Galerkin, conduz ao método Gradiente Conjugado (CG) e o Método de Ortogonalização Completa (FOM). A abordagem da norma mínima residual produz métodos como o Mínimo Resíduo (MINRES) e o Resíduo Mínimo Generalizado (GMRES). A terceira abordagem, PetrovGalerkin, leva a métodos como o Gradiente Bi-conjugado (BICG) e o Resíduo Quase Mínimo (QMR). Na quarta abordagem leva a métodos como o LQ Simétrico (SYMMLQ). Além dessas quatro abordagens, existem métodos híbridos, resultado da mistura das abordagens anteriores, tais como o Gradiente Conjugado Quadrado (CGS), o Gradiente Bi-conjugado Estabilizado (BICGSTAB) e o Resíduo Quase Mínimo Livre de Transpostas (TFQMR). Com o intuito de selecionar um método iterativo para ser parte do solucionador proposto, são apresentadas a seguir as principais características dos métodos iterativos mencionados acima, destacando-se suas vantagens e desvantagens. Contudo, enfatizam-se aquelas características mais convenientes para solução dos sistemas lineares cujas matrizes Jacobianas apresentam estruturas e propriedades espectrais típicas de fluxo de carga. Adicionalmente, são avaliados os principais métodos iterativos em termos de robustez e eficiência computacional. 2.5 A Escolha do Método Iterativo Conforme exposto anteriormente, a escolha do método iterativo é fundamental para o desenvolvimento de um solucionador adequado para o problema de interesse. Para tal, foram investigados muitos dos avanços ocorridos no âmbito dos métodos de Subespaço de Krylov durante os últimos anos. A informação apresentada a seguir, geralmente não disponível nos livros, provém de publicações realizadas pelos próprios autores dos métodos iterativos e de estudos sobre recentes desenvolvimentos computacionais realizados em métodos de subespaços de Krylov. 26 2.5.1 Características dos Principais Métodos Iterativos Na Tabela 2.3 apresentam-se os principais métodos iterativos e suas características mais relevantes. Na segunda coluna é indicado o tipo de matriz de coeficientes que cada método é capaz de solucionar. Esta primeira característica descarta os métodos CG, SYMMLQ e MINRES, dado que a Jacobiana do problema de fluxo de carga não é numericamente simétrica. Embora os métodos restantes sejam capazes de solucionar sistemas não-simétricos, apenas o método GMRES não apresenta maiores dificuldades quando as matrizes de coeficientes são indefinidas. Deve-se lembrar que, nos sistemas elétricos de maior-porte (IEEE300-barras, Norte-Nordeste e brasileiros) aparecem autovalores com parte real negativa e também outros com parte real positiva e para estes sistemas, as matrizes Jacobianas são indefinidas. Portanto, teoricamente o método GMRES deve apresentar melhor comportamento durante a solução dos sistemas lineares de fluxo de carga reais de médio e maior porte. Tabela 2.3 – Principais características dos métodos iterativos. Método Iterativo CG Tipo de Matriz de Algoritmo de Construção Condições de Busca da Coeficientes do Subespaço de Krylov Solução Lanczos Ritz-Galerkin Simétrica definida positiva BICG Não-simétrica Bi-Lanczos Petrov-Galerkin SYMMLQ Simétrica indefinida Lanczos Norma mínima do erro MINRES Simétrica indefinida Lanczos Norma mínima residual 27 FOM GMRES Não-simétrica Não-simétrica indefinida Arnoldi Ritz-Galerkin Arnoldi Norma mínima residual CGS Não-simétrica Lanczos Não-simétrico Híbrido QMR Não-simétrica Lanczos antecipado Petrov-Galerkin BICGSTAB Não-simétrica Bi-Lanczos Modificado Híbrido TFQMR Não-simétrica Lanczos Não-simétrico Híbrido Outra característica importante, apresentada na terceira coluna da Tabela 2.3, está associada com o algoritmo adotado pelo método iterativo para que, a partir da matriz de coeficientes A e do resíduo r0, sejam construídas as bases ortogonais do Subespaço de Krylov κk, dado por (2.19). Basicamente, o Subespaço de Krylov pode ser construído usando-se o algoritmo de Arnoldi para matrizes não-simétricas, o algoritmo de Lanczos para matrizes simétricas ou usando-se alguma das várias versões do algoritmo de Lanczos não-simétrico, por exemplo, biortogonalização de Lanczos (Bi-Lanczos ou também chamado Two-side Lanczos) e Lanczos antecipado (look-ahead Lanczos). Em cada iteração destes algoritmos é calculada uma base ortogonal do subespaço de Krylov. O tipo de algoritmo adotado pelos métodos iterativos influencia bastante na sua robustez, podendo algumas vezes interromper o processo iterativo devido a divisões por zero. Tanto o algoritmo de Arnoldi (usado pelos métodos FOM e GMRES) quanto o algoritmo de Lanczos (usado pelos métodos CG, SYMMLQ e MINRES) não apresentam interrupções (breakdowns) nem problemas de estabilidade numérica durante sua execução, consequentemente, os métodos FOM, GMRES, CG, SYMMLQ e MINRES podem ser considerados os mais robustos. No entanto, no caso dos métodos que usam variações ou adaptações do algoritmo de Lanczos não-simétrico (BICG, CGS, QMR, BICGSTAB e TFQMR), podem-se apresentar dois tipos de interrupções conhecidas como: interrupções de Lanczos (Lanczos-breakdowns ou também chamadas ghostbreakdown) e interrupções pivô (pivot-breakdown ou também chamadas truebreakdowns). Em ambos os casos o algoritmo usado para calcular as bases ortogonais do subespaço de Krylov é interrompido devido a uma divisão por zero. Quando acontece uma interrupção do tipo Lanczos, a divisão por zero se produz nas relações de recorrência usadas no algoritmo, porém, a interrupção do tipo pivô ocorre quando a base do Subespaço de Krylov, que está sendo calculada, não existe. Quando a divisão não é por zero, mas por um número muito pequeno (próximo de zero), ocorrem as chamadas quase-interrupções (near-breakdown), as quais introduzem erros 28 de arredondamento no algoritmo de Lanczos e, consequentemente, problemas de estabilidade numérica no processo de solução do método iterativo. De acordo com o tipo de interrupção à qual estão associadas, as quase-interrupções são conhecidas na literatura como quase-interrupção de Lanczos e quase-interrupção pivô. Quando ocorre uma quase-interrupção o algoritmo não é interrompido, mas a cada iteração são introduzidos erros que, quando acumulados produzem instabilidade numérica. O método iterativo normalmente para quando atinge o máximo número de iterações sem ter encontrado a solução do sistema de equações lineares. Uma vez apresentadas as principais características dos métodos iterativos, é necessário determinar se estas favorecem ou não os respectivos desempenhos quando usados para solucionar o subproblema linear de fluxo de carga, principalmente em termos de robustez. A seguir, a partir das características expostas na Tabela 2.3, são apresentadas as principais vantagens e desvantagens dos métodos iterativos mencionados anteriormente, na Tabela 2.4. 2.5.2 Vantagens e Desvantagens dos Principais Métodos Iterativos Conforme mencionado anteriormente, os métodos CG, SYMMLQ e MINRES não são adequados para solucionar os sistemas lineares de fluxo de carga, já que são aptos apenas para solucionar sistemas simétricos. A primeira tentativa para solucionar este impasse é apresentada em, com o método conhecido como Gradiente Bi-conjugado (BICG) sendo posteriormente melhorado por. Este método, na sua versão final, consegue solucionar sistemas não-simétricos, porém apresenta muitas desvantagens associadas com seu péssimo desempenho, principalmente em termos de robustez. O BICG utiliza uma das formas do algoritmo Lanczos não-simétrico, conhecida como biortogonalização de Lanczos ou simplesmente Bi-Lanczos, para construir as bases ortogonais do Subespaço de Krylov. No entanto, este algoritmo pode apresentar todos os tipos de interrupções e quase-interrupções descritas anteriormente, e adicionalmente pode também apresentar outro tipo de interrupção (divisão por um pivô zero) associado com a fatoração LDU requerida internamente pelo algoritmo Bi-Lanczos. Todos estes problemas podem ocasionar a falha do BICG ou provocar instabilidades numéricas que impedem sua convergência. Além disso, o BICG precisa realizar, em cada iteração, uma multiplicação matriz-vetor com cada matriz A e AT, aumentando o custo computacional. Contudo, durante o processo iterativo a taxa de convergência varia drasticamente 29 (aumenta e diminui) de uma iteração para outra. Portanto, o BICG apresenta uma convergência irregular e oscilatória. Pode-se considerar o BICG como a primeira tentativa de se solucionar problemas nãosimétricos. Porém, devido à sua baixa robustez e ineficiência computacional, este método está descartado para o solucionador proposto. Os problemas que o BICG produz foram posteriormente motivo de intenso estudo pelos pesquisadores, cujo resultado propõe duas alternativas, sendo estas: o método Gradiente Conjugado Quadrado (CGS) e o método Resíduo Quase Mínimo (QMR). No caso do CGS, o método requer apenas uma operação matriz-vetor por iteração. Esta mudança torna o método CGS computacionalmente mais eficiente que o BICG. Porém, as interrupções e quase-interrupções ainda podem ocorrer de forma similar ao método BICG. Além disso, as mudanças realizadas no CGS amplificam mais ainda as irregularidades e oscilações da sua convergência. O CGS apenas melhorou a eficiência computacional do BICG, portanto, os problemas de robustez não solucionados são determinantes para que o método CGS também não seja adotado como o método do solucionador proposto. No caso do QMR, a atenção dos pesquisadores esteve direcionada a solucionar os problemas de robustez, sendo proposto um novo algoritmo de Lanczos não-simétrico conhecido como Lanczos antecipado (look-ahead Lanczos). Este algoritmo reduz a ocorrência de interrupções e problemas de estabilidade numérica, mas a robustez do método vai depender das características e propriedades espectrais dos sistemas lineares do problema que está sendo solucionado. Embora as mudanças realizadas no algoritmo de Lanczos tornem este método mais robusto, a ocorrência de interrupções é ainda possível. Estas mudanças eliminam qualquer tipo de irregularidade ou oscilações da sua convergência (convergência suave). Contudo, o QMR apresenta uma taxa de convergência similar ao BICG e é pouco eficiente, já que precisa realizar duas operações de multiplicação matriz-vetor em cada iteração. 30 Tabela 2.4 – Principais vantagens e desvantagens dos métodos iterativos, quando usados para solucionar o subproblema linear de fluxo de carga. Método Iterativo Vantagens CG Robusto em problemas simétricos. Falha em sistemas simétricos indefinidos ou não-simétricos. BICG Soluciona sistemas não-simétricos. Realiza multiplicações matriz-vetor com ambas as matrizes A e A . Pode parar quando apresentar interrupções do tipo Lanczos e pivô. O algoritmo Bi-Lanczos também pode parar se apresentar divisões por zero, dado que implicitamente realiza uma fatoração triangular do tipo LDU. Apresenta instabilidade numérica devido às possíveis quase-interrupções. Convergência oscilatória e irregular, tornando-o não confiável. SYMMLQ Robusto em problemas simétricos. Falha quando soluciona sistemas não-simétricos. MINRES Robusto em problemas simétricos. Falha quando soluciona sistemas não-simétricos. FOM Robusto inclusive em sistemas não-simétricos. Computacionalmente menos eficiente que o GMRES. GMRES Robusto, garantindo a convergência ainda quando soluciona sistemas mal-condicionados. Melhor controle do erro que o método FOM. Muito mais eficiente computacionalmente que os métodos FOM, BICG, CGS, QMR, TFQMR e BICSTAB. CGS Computacionalmente mais eficiente que o BICG. QMR A convergência não é oscilatória nem irregular. O algoritmo Lanczos antecipado reduz a ocorrência de interrupções. BICGSTAB Computacionalmente mais eficiente que os métodos BICG, QMR, CGS e TFQMR. Reduz as possibilidades de ocorrência de interrupções e quase-interrupções apresentadas no BICG e CGS. Não possui convergência oscilatória nem irregular. TFQMR Apenas realiza multiplicações matriz-vetor com a matriz A. Desvantagens T Menos eficiente computacionalmente que os métodos diretos, quando usado sem nenhuma estratégia numérica adicional. A convergência irregular e oscilatória do método BICG é amplificada. Apresenta taxa de convergência similar ao BICG. Computacionalmente pouco eficiente. Ainda pode apresentar interrupções no algoritmo Lanczos antecipado. T Realiza multiplicações matriz-vetor com ambas as matrizes A e A . Ainda pode apresentar interrupções e quase-interrupções que podem parar o processo iterativo ou impedir a convergência (instabilidade numérica). Perde eficiência quando o espectro da matriz de coeficientes apresenta autovalores com grandes quantidades imaginárias, como as encontradas nos espectros dos sistemas elétricos de médio e maior-porte. Susceptível a apresentar interrupções no algoritmo Lanczos. Computacionalmente pouco eficiente e pouco robusto. 48 30 31 Este método também não é recomendável, dado que não garante robustez. Ambos os métodos CGS e QMR não conseguem satisfazer completamente os requerimentos de eficiência e robustez. Portanto, são propostos dois novos métodos para amenizar as deficiências dos métodos CGS e QMR. O primeiro método, proposto em é conhecido como o Gradiente Bi-conjugado Estabilizado (BICGSTAB) e o segundo método, proposto em, como o Resíduo Quase Mínimo Livre de Transpostas (TFQMR). O BICGSTAB surge como uma variante do CGS, modificando o algoritmo Bi-Lanczos com a intenção de controlar e minimizar melhor o resíduo. Com isso, consegue-se uma convergência mais suave e rápida que a obtida pelo CGS, evitando-se irregularidades e oscilações. As mudanças realizadas conseguem também reduzir as interrupções e instabilidades numéricas, embora ainda sejam relatados casos onde o processo de convergência deste método é interrompido. Neste método, como acontece no CGS, não é necessário calcular explicitamente a transposta da matriz de coeficientes, o que permite obter um custo computacional aceitável. No entanto, o BICGSTAB ainda requer realizar dois produtos matriz-vetor em cada iteração. O método BICGSTAB pode ser considerado como a evolução mais atual e bem sucedida dos métodos baseados nos algoritmos de Lanczos. Em relação aos métodos BICG, CGS e QMR, o método BICGSTAB é superior, já que melhora o desempenho em termos de robustez e eficiência computacional. Porém, este método não é confiável devido às possíveis interrupções que não garantem robustez. Este método ainda apresenta outra séria desvantagem, perde eficiência quando o espectro da matriz de coeficientes apresenta autovalores com grandes quantidades imaginárias, como as encontradas nos espectros dos sistemas elétricos de maior porte. O método TFQMR, embora tenha sido melhorado reduzindo-se a ocorrência de interrupções do algoritmo Lanczos e evitando-se a multiplicação matriz-vetor com a transposta da matriz de coeficientes, ainda apresenta interrupções e instabilidades numéricas, portanto, também não será considerado como o método iterativo do solucionador proposto. Conforme mencionado anteriormente, os métodos baseados no algoritmo de Arnoldi (FOM e GMRES) não apresentam problemas de interrupções nem quase-interrupções, como os apresentados nos métodos baseados nos algoritmos de Lanczos. Embora o método FOM seja robusto, não supera o desempenho do método GMRES. O método GMRES, além de garantir a convergência, ainda quando os sistemas lineares são malcondicionados, é computacionalmente muito mais eficiente que todos os métodos 31 32 anteriores. Portanto, é o método mais indicado para ser usado no solucionador iterativo proposto. 2.5.3 Avaliação do Desempenho dos Principais Métodos Iterativos Sem Estratégias Nos parágrafos anteriores foram apresentadas as principais características numéricas dos métodos iterativos. As informações foram consideradas para determinar os efeitos que estas podem produzir na robustez e eficiência computacional dos métodos iterativos, quando usados para solução de sistemas de equações lineares de fluxo de carga. Depois de realizada esta análise, foi possível apresentar as vantagens e desvantagens de cada método iterativo. A seguir são apresentadas simulações computacionais para avaliar e comparar o desempenho destes métodos iterativos. Os métodos iterativos considerados nestas simulações estão disponíveis no programa computacional SPARSKIT, sob a forma de solucionador de domínio público, sendo estes o BICG, BICGSTAB, TFQMR, FOM e GMRES. Os testes computacionais considerando estes métodos envolvem simulações para resolver os sistemas lineares de fluxo de carga associados com as matrizes Jacobianas da primeira iteração Newton-Raphson, apresentadas na Tabela 2.2, para dois cenários de carga: caso base e ponto próximo do máximo carregamento (PMC), ou seja, um sistema linear para cada cenário, resultando em dois sistemas lineares por cada sistema elétrico. Para todos os métodos iterativos foi usado o mesmo tipo de teste de convergência e as mesmas tolerâncias (absoluta e relativa), ambas iguais a 10 -5. As matrizes Jacobianas são armazenadas no formato de armazenamento esparso conhecido como compressão esparsa por linha (CSR) para serem compatíveis com o solucionador SPARSKIT. Para estas simulações foram considerados alguns dos sistemas de energia apresentados e os métodos CGS e QMR não foram usados por serem versões anteriores aos métodos BICGSTAB e TFQMR, respectivamente. O objetivo principal é verificar a eficiência dessas técnicas numéricas sob dois aspectos: eficiência computacional e robustez. A análise dos resultados permite confirmar, na prática, que o método GMRES, o qual apresenta as melhores vantagens, é realmente o mais apropriado para solucionar o subproblema linear de fluxo de carga. Na Tabela 2.5 apresenta-se o número de operações de ponto flutuante em MFLOPS (10 6 FLOPS) resultantes da simulação computacional de ambos os casos (base e PMC) referentes a cada método. 32 33 Tabela 2.5 – Número de operações de ponto flutuante para os métodos BICG, FOM, GMRES, BICGSTAB e TFQMR (MFLOPS). Sistema Cenário BICG FOM GMRES BICGSTAB TFQMR C. Base 1,9206 0,8421 0,5846 4,0508 4,6651 PMC 2,8020 1,2726 1,0801 Interrupção 2,9828 C. Base 253,6703 233,7633 134,4423 523,4859 549,7364 PMC 366,8282 260,5156 193,1021 548,7201 557,2733 56,2609 51,1462 IEEE30-barras IEEE118-barras C. Base Interrupção IEEE145-barras PMC Interrupção 58,8182 55,6415 Quase- Quase- interrupção Quase- interrupção Quase- interrupção interrupção C. Base 39,3701 33,8631 24,3895 40,4352 45,4439 PMC Quase- 48,0227 41,0450 54,4859 69,0345 705,9867 624,7670 Quase- Quase- interrupção Quase- interrupção Quase- interrupção interrupção 457,1362 Interrupção Interrupção 472,7717 Interrupção Interrupção IEEE162-barras C. Base interrupção Quase- interrupção QuasePMC 749,3342 interrupção QuaseC. Base 493,7071 Norte-Nordeste interrupção de 274 barras QuasePMC 524,7766 interrupção IEEE300-barras 635,0290 A primeira (e a mais importante) observação está associada à robustez dos métodos FOM e GMRES. Ambos os métodos, baseados no algoritmo de Arnoldi, solucionaram todos os sistemas lineares sem ter nenhum tipo de interrupção ou problema de estabilidade numérica. No entanto, o método GMRES apresenta melhor eficiência computacional (menores tempos de processamento). Todos os outros métodos baseados nos algoritmos de Lanczos não-simétrico falharam, apresentando as já comentadas interrupções e quase-interrupções responsáveis pela parada imediata do processo de convergência do método iterativo, enquanto que as quase-interrupções realizam iterações indefinidamente. Neste caso, o processo iterativo só para quando é atingido o número máximo de iterações pré-estabelecido, sem ter encontrado a solução. Observa-se que estes problemas são mais frequentes no ponto PMC, quando as matrizes estão próximas da singularidade; também acontecem quando solucionados sistemas lineares associados com redes elétricas de maior porte. Estes resultados também mostram que os sistemas lineares do problema de fluxo de carga são de difícil solução (mal condicionados), inclusive para o método BICGSTAB, 33 34 que provou eficiência e robustez em outras aplicações da engenharia e que também é adotado por alguns pesquisadores da área de sistemas elétricos. Na Figura 2.8 apresentam-se os valores de taxa de convergência dos cinco métodos iterativos avaliados, agrupados em barras de cores diferentes, para cada sistema elétrico no caso base. Observa-se que, os métodos GMRES (barra verde) e o FOM (barra azul) apresentam maiores taxas de convergência que os outros métodos, por isso a convergência é mais rápida. Portanto, são computacionalmente mais eficientes. Neste trabalho, a correlação entre o aumento da taxa de convergência e da eficiência computacional será muito utilizada para análise dos resultados. Observa-se também que as taxas de convergência negativas correspondem aos casos não convergentes, quando os métodos experimentaram interrupções (I) ou quase-interrupções (QI). Nota-se que estes problemas são mais frequentes a medida que a dimensão do sistema aumenta. Estes resultados mostram que o GMRES além de ser eficiente, também se mostrou mais robusto. Figura 2.8 – Taxa de convergência dos métodos BICG, FOM, GMRES, BICGSTAB e TFQMR no CASO BASE. 34 35 Na Figura 2.9 apresentam-se os valores de taxa de convergência dos cinco métodos iterativos avaliados, agrupados em barras de cores diferentes para cada sistema elétrico no ponto PMC. De forma similar ao caso base, observa-se que os métodos GMRES e FOM apresentam maiores taxas de convergência que os outros métodos. Observa-se também que, no PMC, os métodos baseados no algoritmo de Lanczos nãosimétrico experimentaram maior número de interrupções e quase-interrupções que no caso base. Nota-se que tanto no caso base como no ponto PMC as taxas de convergência diminuem quando a dimensão do sistema aumenta. Contudo, o método GMRES não mostrou nenhum tipo de dificuldade para solucionar os sistemas próximos ao ponto de máximo carregamento, próximos da singularidade, mal condicionados e de tamanhos diversos. Finalmente, investigadas as principais características dos métodos iterativos, destacando suas vantagens e desvantagens, realizadas avaliações individuais de desempenho para solução de típicos sistemas lineares de fluxo de carga, conclui-se que o GMRES deve integrar o solucionador proposto. A opção por este método baseou-se nos seguintes fatos: Figura 2.9 – Taxa de convergência dos métodos BICG, FOM, GMRES, BICGSTAB e TFQMR no PMC. 35 36 • Outros métodos (CG, MINRES e o SYMMLQ) tão robustos quanto o GMRES são descartados devido à não-simetria numérica da matriz Jacobiana de fluxo de carga. Outro fator determinante é que alguns dos sistemas elétricos considerados apresentaram matrizes Jacobianas indefinidas. Uma vez que o método GMRES é o único que garante a convergência quando se soluciona sistemas lineares com matrizes indefinidas, então, os métodos restantes foram descartados devido à baixa confiabilidade. • O algoritmo de Arnoldi usado pelo GMRES para gerar as bases do subespaço de Krylov é robusto e não apresenta nenhum tipo de interrupção, nem quase-interrupções. Esta característica torna o método GMRES muito robusto quando soluciona sistemas mal condicionados. As interrupções e quase-interrupções que aconteceram nos métodos baseados nos algoritmos de Lanczos não-simétrico descartam os métodos BICG, BICGSTAB e TFQMR. • A análise dos resultados, em termos de robustez e eficiência computacional, demonstrou que o método GMRES apresenta características mais adequadas para a solução do subproblema linear de fluxo de carga. O método GMRES apresentou a melhor (maior) taxa de convergência e o menor número de operações de ponto flutuante quando comparado com os outros métodos, sendo, portanto, o método iterativo mais eficiente. Estudos preliminares encontrados na literatura especializada comparam o GMRES com outros métodos e concluem que o GMRES é muito robusto e de melhor desempenho quando aplicado ao problema de fluxo de carga. Contudo, os métodos iterativos sem nenhum tipo de estratégia numérica ainda são considerados computacionalmente pouco eficientes quando comparados com os métodos diretos. Um método iterativo começa a ser tão eficiente quanto um método direto quando uma estratégia eficiente de pré-condicionamento for usada. Vale ressaltar que não é objetivo aqui propor um novo método iterativo (sem estratégias), já que provavelmente este seria tão ineficiente quanto os métodos iterativos existentes. Também não é objetivo modificar o algoritmo de um método iterativo existente, pois as melhorias na sua eficiência computacional seriam mínimas comparadas com aquelas que podem ser conseguidas usando-se estratégias numéricas. Como já mencionado, o objetivo é propor um solucionador iterativo (método iterativo e estratégias numéricas) para solução do subproblema linear de fluxo de carga, sendo que o foco principal do trabalho é reunir e propor adequadamente estratégias numéricas, visando aproveitar melhor os recursos e características do método iterativo a fim de melhorar significativamente sua eficiência computacional e ainda mais sua robustez. A 36 37 seguir, faz-se uma abordagem teórica dos conceitos associados ao método GMRES, explicando-se de forma ilustrativa os processos envolvidos no seu algoritmo. 2.6 O Método do Resíduo Mínimo Generalizado – GMRES Em é introduzido o método iterativo do Resíduo Mínimo Generalizado (GMRES) para resolver sistemas de equações lineares não simétricos. O método GMRES baseia a busca da aproximação na condição da norma residual mínima, identificando x k de tal forma que sua correspondente norma Euclidiana do resíduo seja mínima sobre o subespaço de Krylov, cujas bases são ortonormalizadas. Neste método estima-se uma solução aproximada xk utilizando-se a expressão (1.17), onde V é uma matriz cujas colunas são os vetores (1.18), base do subespaço de Krylov após sua ortonormalização. O GMRES, em uma primeira etapa, usa o processo de Arnoldi para ortonormalizar os vetores base do subespaço de Krylov КK(A,r0). xk = x0 +V  yk  V  nxm = v1f (2.27) v 2f ... vkf   (2.28) A Equação (2.29) é usada pelo algoritmo de Arnoldi para ortonormalizar o vetor v ik+1 em relação aos outros k vetores (vf1, vf2, ..., vfk), previamente ortonormalizados. f vk+1 = k T   i   v gf  v ik+1-     v gf   v k+1      g=1   T   i   v gf  v ik+1-     v gf   v k+1      g=1   (2.29) k 2 Na Figura ilustra-se o processo de ortonormalização de Arnoldi (também conhecido como método Gram-Schmidt modificado), onde os vetores (vi1, vi2, ..., vik) pertencem ao subespaço original de Krylov (r0, A·r0, A2·r0, ..., Ak·ro). Uma vez estimados V e seus vetores base, o GMRES em uma segunda etapa determina o valor de yk baseando-se na condição da norma mínima residual, apresentada na seção 2.3. Deve-se minimizar a norma euclidiana do resíduo para conseguir uma boa aproximação da solução x* (ver equação 2.30). Utilizando-se (2.27) e (2.29) pode-se demonstrar que a norma euclidiana do resíduo pode também ser representada por (2.31). 37 38 Figura 2.20 – Ortonormalização do Subespaço de Krylov. Em cada iteração k do GMRES, o valor mínimo esperado para a norma do resíduo é zero ou muito próximo de zero, sendo suficiente para satisfazer o teste de convergência. Portanto, os valores dos elementos do vetor yk devem ser aqueles que tornam a expressão (2.31) igual a zero (ou muito próximo de zero); para conseguir isso, a expressão ||r0||2·ê1-Hk·yk deve ser igual a zero. O sistema linear resultante a ser solucionado é apresentado em (2.32) e representa um típico problema de mínimos quadrados de pequeno porte, solucionado usando-se um processo de fatorações ortogonais baseadas em rotações no plano ou método de Rotação de Givens. rk 2 = b-A  xk rk 2 = Vk+1   ro 2 (2.30) 2  eˆ 1-Hk  yk  Hk  yk = ro 2  eˆ 1 2 (2.31) (2.32) Onde ê1 é o vetor canônico de dimensão k; Hk é a matriz de Hessenberg superior (k+1)×k apresentada em (2.33), cujos elementos são calculados durante a ortonormalização pelo processo de Arnoldi, usando (2.34); e ||r0||2∙ê1 é conhecido na literatura como o vetor do lado direito. 38 39    h(1,1)   h(1,2)    h    h(2,2)  (2,1)    h    (3,2)  Hk                   0   0   h h2 1    h(1,k-1)   h(1,k)                  h    (k,k-1)   h(k,k)       0   h(k+1,k)    hk hk-1  h(i,k) =  v gf    A  vkf      (2.33) T k T   h(k+1,k) = A  vkf -    v gf    A  vkf    v gf      i=1  (2.34) Uma vez determinado o subespaço V e o vetor y k, em uma terceira e última etapa são calculados a aproximação xk e o resíduo rk, usando-se as expressões (2.27) e (2.30), respectivamente. As três etapas são repetidas em cada iteração k dentro de um laço que termina quando o teste de convergência for satisfeito. Na Figura 2.10, apresentam-se em forma de fluxograma e de forma sucinta, todas as etapas do processo iterativo do método GMRES. Figura 2.10 – Fluxograma do método iterativo GMRES. Duas grandezas podem ser usadas no teste de convergência dado por (2.35), sendo a norma-2 do resíduo ||r||2 e a norma-2 da variação da solução aproximada ||xk-xk-1||2. Aqui, o teste de convergência é baseado no cálculo da norma-2 do resíduo, conforme definido em (2.30) e usado no passo 25 do algoritmo. Os parâmetros rtol e atol são as tolerâncias relativa e absoluta do método iterativo, respectivamente. r 2  rtol  b  atol 2 (2.35) Para ilustrar melhor o significado de V e y k, pode-se expressar (2.27) de outra forma, como apresentado em (2.36). Como se pode observar, o método iterativo procura a aproximação xk, partindo da aproximação inicial x0 (normalmente x0=0), mediante a soma 39 40 de sucessivos vetores (ε 1·vf1, ε2·vf2, ..., εk·vfk) em cada iteração k, onde (ε 1, ε2, ..., εk) são os elementos do vetor coluna yk apresentado em (2.37). Portanto, pode-se afirmar que a solução xk é uma combinação linear dos vetores da base ortonormalizada V, ou seja, a solução xk é a resultante da soma desses vetores, como ilustrado na Figura 2.11. Neste caso, para um sistema linear de duas incógnitas (neste gráfico x 0≠ 0) foram necessárias apenas duas iterações para encontrar a solução exata x*. xk = x0 +  v1f v 2f   ε1    f   ε2  ... v k  = ε   v f  + ε   v f  + ... + εk   v kf     1  1 2  2      εm  yk = ε1 ε 2 ... εk  T (2.36) (2.37) Figura 2.11 – Processo de aproximação à solução x*. A computação de yk requer a solução de um problema de mínimos quadrados de dimensão (k+1)  k, onde k é o número da iteração GMRES. Portanto, estimar y k pode se tornar uma tarefa difícil ao solucionar sistemas lineares de grande-porte quando o número de iterações aumentar. Dentre algumas das características apresentadas pelo GMRES, destaca-se a convergência em, no máximo, n iterações (onde n é o número de incógnitas e a dimensão da matriz de coeficientes). 40 41 Obviamente, se são consideradas matrizes em que n é grande, no pior caso (convergência na n-ésima iteração), essa garantia de convergência não será vantagem, pois o processo de Arnoldi irá requerer o armazenamento de n vetores de dimensão-n, além de uma matriz de dimensão (n+1)  n, o que exigirá n produtos matriz-vetor, com uma complexidade espacial (espaço de memória) O(n 2) e temporal (tempo de execução) O(n3). Isso poderá fazer com que não seja possível obter a solução do sistema num determinado computador, pois a memória disponível não seria suficiente; ou a resposta não seria obtida em tempo hábil, causando dificuldades no desempenho do GMRES. As dificuldades associadas ao baixo desempenho do GMRES (sem estratégias) quando comparado com os métodos diretos são exemplificadas a seguir. 2.6.1 Exemplo Ilustrativo do GMRES (Sem Estratégias) Para ter uma ideia real das dificuldades associadas ao método GMRES quando usado sem nenhuma estratégia numérica, considere o sistema-teste IEEE118-barras. Deve-se solucionar um sistema linear com uma matriz de coeficientes (Jacobiana) de dimensão igual a (236  236) em cada iteração Newton-Raphson do fluxo de carga. Neste caso, para cada sistema linear é necessário realizar no máximo 236 iterações do GMRES, o que significaria que na última iteração seriam armazenados 236 vetores de dimensão 236, além de uma matriz de Hessenberg (2.33) de dimensão 237  236. Na Figura 2.12 são apresentados dois gráficos de barras. As primeiras 4 barras (lado esquerdo) correspondem ao número de operações realizadas ao longo da solução do problema (em milhões de operações de ponto flutuante - MFLOPS). As segundas barras (lado direito) correspondem ao número de iterações lineares para cada iteração do método Newton-Raphson (1, 2 e 3). As barras pretas e a última linha da tabela apresentam o número total de MFlops e o número total de iterações da solução do problema. O tempo computacional (CPU) requerido para realizar 534,1593 milhões de operações é aproximadamente 6,02 segundos. 41 42 Figura 2.12 – Desempenho computacional do método GMRES. Na Figura 2.13 apresenta-se em gráfico de barras o número de operações realizadas pelo método direto (MA28), ao longo da solução do problema (em milhões de operações de ponto flutuante - MFLOPS), para cada iteração do método Newton-Raphson (1, 2 e 3). A barra preta e a última linha da tabela apresentam o número total de MFLOPS realizados. O tempo computacional requerido para realizar 89.432 operações é insignificante, aproximadamente 0,00001 segundos. Figura 2.13 – Desempenho computacional do método direto. Comparando-se os gráficos da Figura 2.12 e Figura 2.13, observa-se que, embora este seja um sistema considerado de pequeno porte, o método GMRES quando usado na sua forma convencional apresenta um desempenho abaixo do esperado, se comparado a um método direto. Sob estas condições, o método GMRES pode se tornar ainda mais 42 43 ineficiente se comparado com os métodos diretos, principalmente se usado para sistemas de maior porte. São apresentadas a seguir estratégias para aumentar a eficiência do GMRES, começando pela de reinicialização que, quando usada, reduz a quantidade e a dimensão das bases do subespaço, o número de operações de ponto flutuante nos produtos matriz-vetor, a dimensão do problema de mínimos-quadrados e o número de iterações necessárias para alcançar a convergência. 2.7 Estratégia de Reinicialização para o método GMRES Uma desvantagem do GMRES é que o esforço computacional e os requisitos de memória por iteração crescem linearmente com processo iterativo e o custo associado se torna rapidamente excessivo. A forma usual de superar esta condição indesejada é usar uma versão do GMRES com estratégia de reinicialização, identificada como GMRES(m). Neste caso, o valor de m limitará a dimensão da base do subespaço de Krylov, como apresentado em (2.38). O parâmetro m se torna relevante quando for muito menor que o valor k de cada iteração linear, pois limita o tamanho e quantidade dos vetores utilizados no algoritmo de Arnoldi e durante a solução do problema de mínimos quadrados. κ m (A,r0 ) = subespaço(r0 , A  r0 , A 2  r0 ,..., Am-1  r0 ) (2.38) A cada iteração k, calcula-se o resíduo e os m vetores ortonormais através do processo de Arnoldi e, a partir deles, resolve-se o problema de mínimos-quadrados obtendo-se uma aproximação xk+1 para a solução do sistema linear. Caso o resíduo r k+1 não seja suficientemente pequeno, faz-se uma nova iteração. Portanto, quanto menor for o valor de m, menor será a quantidade e a dimensão dos vetores bases do subespaço de Krylov a serem ortonormalizadas. Isto resulta em um menor o número de operações de ponto flutuante realizado nos produtos matriz-vetor e durante a solução do problema de mínimos-quadrados. A escolha de um valor apropriado para m não é uma tarefa fácil, já que se m for muito pequeno, o processo associado ao GMRES pode ser lento para convergir, ou falhar completamente e não convergir. Por outro lado se o valor de m for grande, acima do necessário, um considerável espaço de memória será usado. Não existem regras específicas para a escolha desse parâmetro, portanto a escolha apropriada do instante de reinicialização é uma questão de prática. Se estratégias de reinicialização não forem consideradas, o GMRES irá convergir em não mais do que n (dimensão da matriz) 43 44 passos, e isto não é prático para grandes valores de n. Além do mais os requisitos de memória e tempo de CPU na ausência de estratégias de reinicialização são excessivos. 2.7.1 Exemplo Ilustrativo – GMRES(m) Reinicializado A fim de ilustrar as vantagens da estratégia de reinicialização no método GMRES, soluciona-se aqui o mesmo problema do exemplo anterior. O procedimento adotado para estimar o valor ótimo de m a fim de tornar o GMRES mais eficiente envolveu diversas simulações computacionais (Figura 2.15). Em problemas de maior-porte, este procedimento pode não ser eficiente, mas aqui como o sistema é de pequeno porte, esta procura é viável. Foram realizadas várias simulações para calcular o número de MFLOPS associado à solução do sistema linear para valores de m compreendidos entre 2 e 236, encontrando-se o menor número de MFlops (20,83 milhões de operações) para m igual a 53. Figura 2.14 – Simulações para busca do tamanho ótimo do subespaço de Krylov. Para este valor de m é necessário realizar apenas 16 iterações do GMRES, o que significa que na última iteração seria necessário armazenar 16 vetores de dimensão 236, além de uma matriz de Hessenberg de dimensão 17  16. Isto resulta numa solução do problema de mínimos quadrados menos complexa se comparada ao exemplo anterior. Deve-se mencionar que o valor de m é estimado apenas para o sistema linear da primeira iteração Newton-Raphson, sendo este considerado como o valor ótimo para o problema. Na Figura 2.15, ilustra-se o desempenho do GMRES(m), na qual é possível observar uma redução significativa no número de operações e no número de iterações lineares (comparar com valores da Figura 2.15). 44 45 Figura 2.15 – Desempenho computacional do método GMRES(m). Embora esta estratégia tenha melhorado o desempenho computacional do método GMRES, os resultados são inferiores aos conseguidos com o método direto. Portanto, o desempenho do solucionador (GMRES e estratégia de reinicialização) continua abaixo do esperado. Como já mencionado, uma desvantagem da estratégia de reinicialização está associada com a dificuldade da busca de um valor apropriado do parâmetro m, para cada sistema linear em cada iteração Newton-Raphson. A escolha do valor ótimo de m é um dos pontos mais críticos para a eficiência do processo iterativo. Outra forma de mitigar o problema seria reduzir diretamente o número de iterações k necessárias para atingir a convergência do método GMRES, uma vez que o tamanho do subespaço de Krylov e o aumento das operações estão associados a este parâmetro. Na primeira iteração do primeiro exemplo, se fosse possível reduzir o número de iterações lineares totais, por exemplo, de 107 para menos de 5 iterações, então seria necessário armazenar apenas 5 vetores de dimensão 236, além de uma matriz de Hessenberg de dimensão 6 x 5. A dificuldade para resolver o problema de mínimos quadrados seria bastante amenizada, melhorando a eficiência computacional do método GMRES. Embora pareça difícil conseguir tal melhoria, existe uma estratégia para diminuir a complexidade do problema e reduzir significativamente o número de iterações do GMRES, conhecida como pré-condicionamento do método iterativo. A seguir, apresenta-se outra estratégia normalmente usada junto às técnicas de précondicionamento, conhecida como reordenamento da matriz de coeficientes (Jacobiana). Esta estratégia costuma melhorar significativamente a convergência dos métodos iterativos reduzindo, consequentemente, o número de iterações e o tempo gasto de processamento. 45 46 2.8 Estratégias de Reordenamento da Matriz Jacobiana Estratégias de reordenamento têm sido muito usadas pelos métodos diretos para reduzir o número de elementos não-nulos e para melhorar a estabilidade numérica durante a fatoração triangular LU. O reordenamento de uma matriz esparsa é realizado mediante permutações de linhas, colunas ou ambas. Entende-se como permutação o produto da matriz de coeficientes com outras matrizes P e Q, onde estas matrizes são permutações da matriz identidade. Se a matriz de permutação P é multiplicada pela matriz de coeficientes, as linhas são trocadas, e caso a multiplicação envolva a matriz Q, as colunas são trocadas. No caso de ambas as multiplicações, tanto as linhas como as colunas são trocadas. Uma vez determinadas as matrizes de permutação P e Q, o sistema em (1.1) pode ser substituído pelo sistema equivalente (reordenado) de (2.39). Se A é estruturalmente simétrica, (isto é, aij  0 se e somente se aji  0 para todos i,j  [1,n]), então é desejável que as permutações sejam simétricas, ou seja, Q = PT para preservar as propriedades da matriz original. Neste caso, o sistema em (1.1) pode ser substituído pelo sistema equivalente de (2.30) (reordenado simetricamente). Este tipo de reordenamento é adotado neste curso para aproveitar a simetria estrutural da matriz Jacobiana de fluxo de carga. P  A  Q  y  P  b, x  Q y (2.39) P  A  PT  y  P  b, x  PT  y (2.40) Quando realizado o reordenamento simétrico, o conjunto de elementos da diagonal da matriz é preservado e os mesmos são deslocados somente dentro da diagonal. O reordenamento simétrico não afeta os autovalores de A e não modifica o número de condicionamento. Se a matriz A não for simétrica, as permutações devem ser realizadas com base na estrutura da matriz simétrica A+AT. No caso dos métodos iterativos, o reordenamento simétrico é aplicado antes da execução do método, alterando estruturalmente a matriz de coeficientes sem modificar os valores de seus elementos e mantendo sempre o sistema resultante equivalente ao sistema linear original. Os benefícios ou vantagens das estratégias de reordenamento sobre a convergência dos métodos de subespaço Krylov só podem ser observados quando usadas às técnicas de pré-condicionamento. Quando não usadas técnicas de précondicionamento, o reordenamento não produz nenhum tipo beneficio para o método iterativo, já que não altera as propriedades espectrais da matriz de coeficientes, ou seja, 46 47 o número de condicionamento e os autovalores da matriz não são modificados pelo reordenamento. A influência das estratégias de reordenamento no desempenho dos métodos iterativos pré-condicionados tem sido estudada por uma série de autores, principalmente de forma experimental, ainda objeto de discussão. O reordenamento melhora a qualidade do pré-condicionador e reduz o tempo necessário para sua construção, gerando subespaços de Krylov mais próximos da solução exata que quando não usado, chegando à convergência em poucas iterações. Portanto, a influência do reordenamento na convergência do método iterativo se dá de forma indireta, através do pré-condicionador. Estratégias clássicas de reordenamento para redução de largura de banda, tais como Reverse Cuthill-McKee (RCM), ordenamento de Sloan, o ordenamento de Gibbs-PooleStockmeyer, além de outras estratégias usadas para redução do número de elementos não-nulos, tais como as variantes do ordenamento de Mínimo Grau (Minimum Degree MD), e o ordenamento Nested Dissection, têm sido muito estudadas, pesquisadas e avaliadas em diferentes áreas da engenharia e da ciência. Porém, tradicionalmente a estratégia de ordenamento do tipo Mínimo Grau, conhecida como Tinney 2, onde o propósito principal é reduzir o número de elementos não-nulos, tem sido preferida pelos métodos diretos para aplicações em sistemas de energia elétrica, apresentando excelentes resultados em termos de eficiência computacional. Com relação ao método RCM, onde além da redução do número de elementos não-nulos também se reduz a largura de banda, baseado na literatura consultada tem sido a técnica mais popular entre os métodos iterativos e aplicado com sucesso em diferentes áreas da engenharia e aplicações matemáticas. Embora o reordenamento RCM apresente resultados interessantes, sua influência na convergência dos métodos iterativos na área de sistemas de energia elétrica tem sido pouco explorada e estudada. Considerando-se que a influência das técnicas de reordenamento no processo iterativo através do pré-condicionador pode ser favorável em termos de eficiência computacional, então, ambos os métodos (RCM e MD) são considerados a fim de determinar qual o mais apropriado para o problema de fluxo de carga. A seguir, uma análise comparativa entre os reordenamentos RCM e MD é realizada para determinar as características, vantagens e desvantagens de cada um. 47 48 2.8.1 Características Principais das Técnicas de Reordenamento MD e RCM São comparadas as principais características dos reordenamentos RCM e MD quando usados para reordenar típicas matrizes Jacobianas de fluxo de carga. As principais características consideradas neste trabalho são: o número de novos elementos não-nulos e a largura de banda. Na Figura 2.16 compara-se, para ambos os tipos de reordenamento RCM (barra cinza) e MD (barra preta), o número de novos elementos não-nulos gerados depois de realizada a fatoração completa das matrizes Jacobianas de fluxo de carga dos sistemas: IEEE118-barras, Norte-Nordeste de 274 barras, sistemas brasileiros de 2256 e 3515 barras. Figura 2.16 – Número e percentagem de novos elementos não-nulos - RCM e MD. Para entender os gráficos de barras da Figura 2.16 assume-se que os novos elementos não-nulos encontrados, quando usado o reordenamento RCM representam sempre 100%. Os valores mostrados nas barras pretas representam a percentagem de novos elementos não-nulos encontrados quando usado o reordenamento MD, em relação aos encontrados quando usado RCM. Na Figura 2.16 apresenta-se também o número de novos elementos não-nulos gerados para cada matriz Jacobiana dos sistemas-teste. Observa-se que cada vez que aumenta a dimensão dos sistemas elétricos, a eficiência do método MD aumenta, reduzindo significativamente o número de novos elementos nãonulos. Para o sistema brasileiro de 3.513 barras, enquanto que o método RCM gerou 87.648 novos elementos não-nulos, o MD gerou apenas 15.307, representando 17% do valor encontrado pelo RCM, 83% menos elementos que o RCM. 48 49 O número de elementos não-nulos pode não ser uma característica importante para os métodos iterativos não pré-condicionados, já que estes métodos não realizam fatorações LU. Porém, é relevante quando consideradas as estratégias de pré-condicionamento baseadas em fatoração triangular incompleta (ILU). Na Figura 2.17, assume-se que as larguras de banda máximas encontradas, quando usado o reordenamento MD representam sempre 100%, portanto, os valores mostrados nas barras cinza representam a percentagem de largura de banda máxima encontrada quando usado o reordenamento RCM, em relação aos encontrados quando usado MD. Na Figura 2.17 é informada a largura de banda máxima para cada matriz Jacobiana dos sistemas-teste. Observa-se que o método RCM reduz significativamente a largura de banda, por exemplo, 89% menos largura de banda que quando ordenada pelo MD, para o caso do sistema brasileiro de 3.513 barras. Figura 2.17 – Largura de banda - RCM e MD. Segundo a literatura, tanto a redução da largura de banda como do número de novos elementos não-nulos podem ser vantajosos ou não, quando usadas em métodos iterativos pré-condicionados. Uma pode ser mais relevante que a outra dependendo das características da matriz de coeficientes associada ao problema que está sendo solucionado. Nas Figuras 2.20 e 2.21 apresentam-se as estruturas de típicas matrizes Jacobiana reordenadas por MD e RCM, tanto para sistemas elétricos de pequeno porte como também para sistemas de maior porte. 49 50 a) Sistema IEEE118-barras (RCM) b) Sistema IEEE118-barras (MD) c) Sistema Norte-Nordeste de 274 barras d) Sistema Norte-Nordeste de 274 barras (RCM) (MD) Figura 2.18 – Estruturas típicas de matrizes Jacobianas reordenadas de sistemas elétricos de pequeno porte. Com o objetivo de mostrar que o reordenamento não produz nenhum beneficio no desempenho dos métodos iterativos quando não são usadas as estratégias de précondicionamento (GMRES puro), são apresentadas a seguir simulações considerando-se os reordenamentos MD e RCM. 50 51 a) Sistema brasileiro de 2.256 barras (RCM) b) Sistema brasileiro de 2.256 barras (MD) c) Sistema brasileiro de 3.513 barras (RCM) d) Sistema brasileiro de 3.513 barras (MD) Figura 2.19 – Estruturas típicas de matrizes Jacobianas reordenadas de sistemas elétricos de maior porte. 2.8.2 Avaliação do Desempenho do Método GMRES Puro Usando Reordenamento da Matriz de Coeficientes São apresentados experimentos numéricos para avaliar o desempenho do método GMRES quando usadas estratégias de reordenamento sem pré-condicionamento. Os testes envolvem simulações para resolver os sistemas lineares de fluxo de carga associados com as matrizes Jacobianas da primeira iteração Newton-Raphson, apresentadas na seção 2.2, dos sistemas hipotéticos do IEEE e o sistema NorteNordeste Brasileiro. As simulações foram realizadas para dois cenários de carga: caso base e ponto próximo do máximo carregamento (PMC), correspondendo a um sistema linear para cada cenário, resultando em dois sistemas lineares para cada sistema elétrico. 51 52 O objetivo principal é verificar qual a influência das estratégias de reordenamento MD e RCM sobre o desempenho do método GMRES puro (sem pré-condicionamento). Para tal, verifica-se a eficiência computacional e a robustez sem reordenamento (S.R.), com o MD e o RCM. Na Tabela 2.6 apresentam-se o número de operações de ponto flutuante em MFLOPS resultantes da simulação computacional de ambos os casos (base e PMC) referentes a cada técnica de reordenamento. No teste de convergência do método GMRES considera-se que tanto a tolerância absoluta, como a relativa, são iguais a 10-5. Observa-se que o reordenamento pode aumentar como também reduzir o número de operações de ponto flutuante realizadas pelo método GMRES. No entanto, em todos os casos a eficiência computacional é similar quer utilize-se reordenamento ou não. Portanto, pode-se concluir que a estratégia de reordenamento não resultou em nenhum tipo de beneficio para o processo iterativo em termos de eficiência computacional com o GMRES não pré-condicionado. Na Figura 2.20 apresentam-se os valores de taxa de convergência do método GMRES na primeira iteração Newton-Raphson, quando não utilizado reordenamento (barra branca), quando usado reordenamento RCM (barra cinza) e quando usado reordenamento MD (barra preta), para cada sistema elétrico no caso base. Observa-se que para todos os sistemas a partir do IEEE-118 barras, os valores das taxas de convergência, com ou sem reordenamento, são similares. Portanto, o reordenamento não produziu nenhum tipo de beneficio na convergência do GMRES puro. 52 53 Tabela 2.6 – Número de operações de ponto flutuante para o método GMRES (MFLOPS). Sistema GMRES Sem GMRES GMRES Reordenamento Com MD Com RCM C. Base 0,5846 0,3990 0,4872 PMC 1,0801 1,3105 1,2811 C. Base 134,4423 125,234 108,984 PMC 193,1021 152,553 184,324 C. Base 51,1462 59,0819 50,0736 PMC 55,6415 61,7359 53,8167 C. Base 24,3895 23,8688 24,0204 PMC 41,0450 31,3902 31,6158 C. Base 624,7670 602,8944 687,7626 PMC 635,0290 781,5147 700,4683 457,1362 451,7398 468,1108 472,7717 620,6830 598,8441 Cenário IEEE30-barras IEEE118-barras IEEE145-barras IEEE162-barras IEEE300-barras Norte-Nordeste de C. Base 274 barras PMC Figura 2.20 – Taxa de convergência do método GMRES no CASO BASE. Na Figura 2.22 apresentam-se os valores de taxa de convergência do método GMRES, quando não utilizado reordenamento (barra branca), quando usado reordenamento RCM (barra cinza) e quando usado reordenamento MD (barra preta), para cada sistema elétrico no ponto PMC. De forma similar ao caso base, observa-se que para todos os sistemas os valores das taxas de convergência, com ou sem reordenamento, são 53 54 similares. Portanto, o reordenamento não produziu nenhum tipo de beneficio na convergência do GMRES puro. Ambas as estratégias estudadas e avaliadas aqui se mostraram ineficazes quando usadas para melhorar o desempenho do método GMRES. No caso da estratégia de reinicialização, o desempenho do método GMRES aumenta, mas os resultados continuam sendo muito inferiores aos conseguidos pelo método direto. A estratégia de reordenamento, quando não usada junto a alguma estratégia de précondicionamento, não resulta em nenhum tipo de beneficio ao solucionador em termos de eficiência computacional. Figura 2.21 – Taxa de convergência do método GMRES no ponto PMC. As dificuldades associadas com o baixo desempenho do método GMRES, quando comparado com os métodos diretos, ainda estão presentes. No entanto, deve-se fazer uso eficiente da estratégia de pré-condicionamento, várias vezes mencionada em alguns pontos importantes aqui. O pré-condicionamento visa melhorar significativamente o desempenho do método GMRES, talvez superando os resultados conseguidos pelos métodos diretos em termos de eficiência computacional e principalmente de robustez, solucionando sistemas lineares associados ao subproblema linear de fluxo de carga, principalmente quando a matriz for mal-condicionada, o que normalmente está associado a um cenário de condição adversa de operação. Portanto, considera-se que a estratégia de pré-condicionamento é importante para processo iterativo, devido a: • A estratégia de pré-condicionamento explora tanto características estruturais como também as propriedades espectrais da matriz Jacobiana de fluxo de carga, melhorando o desempenho do processo iterativo. O beneficio em termos de eficiência 54 55 computacional é devido ao fato do pré-condicionador aproveitar a esparsidade e simetria estrutural das matrizes do problema. O beneficio em termos de robustez ocorre porque o pré-condicionador permite solucionar sistemas lineares com matrizes indefinidas (que apresentam autovalores tanto com parte real negativa como também positiva), típicas nos sistemas de maior porte, como o Sistema Interligado Nacional (SIN). Este tipo de sistema pode representar cenários de difícil convergência para um método iterativo. • A função principal da estratégia de pré-condicionamento é substituir o sistema linear a ser solucionado por outro equivalente, cuja matriz de coeficientes apresente menor número de condicionamento e autovalores mais afastados da origem do plano complexo. Dado que a matriz do problema de interesse, dependendo da sua dimensão e das condições de operação pode estar mal-condicionada e próxima da singularidade, como visto em exemplos anteriores, o pré-condicionamento é uma estratégia para ser considerada no processo iterativo. • O sistema equivalente conseguido pelo pré-condicionador permite construir subespaços de Krylov cuja aproximação da solução está muito próxima da solução exata. Consequentemente, o método GMRES precisa realizar poucas iterações (muitas vezes apenas uma) até satisfazer o teste de convergência. Uma convergência rápida (em poucas iterações) ameniza a dificuldade associada com a solução do problema de mínimos quadrados do algoritmo do GMRES, melhorando a sua eficiência computacional e tornando desnecessária a estratégia de reinicialização. • O pré-condicionamento permite aproveitar a estratégia de reordenamento da matriz Jacobiana. Os benefícios ou vantagens das estratégias de reordenamento sobre a convergência dos métodos de subespaço Krylov só podem ser observados quando usadas técnicas de pré-condicionamento. O reordenamento permite melhorar a qualidade do pré-condicionador e reduz o tempo gasto necessário para sua construção. O précondicionador construído após o reordenamento da matriz de coeficientes gera subespaços de Krylov mais próximos da solução exata que quando não usado o reordenamento, conseguindo atingir a convergência em um menor número de iterações. Portanto, a influência do reordenamento na convergência do método iterativo se dá de forma indireta, através do pré-condicionador. 55 56 3 3.1 Estratégias de Pré-condicionamento Introdução A solução de sistemas esparsos do subproblema linear do fluxo de carga foi enfocada no tópico anterior. O estudo das características e propriedades espectrais de matrizes Jacobinas típicas de fluxo de carga mostrou que, as matrizes de coeficientes destes sistemas lineares são normalmente indefinidas, mal condicionadas e muitas vezes próximas da singularidade, especialmente quando os sistemas elétricos são de maiorporte e operam próximos ao máximo carregamento. Estas características dificultam o processo iterativo de solução levando a uma convergência lenta, ou até mesmo à falha do processo, justificando a necessidade de se considerar estratégias que atenuem ou eliminem estas dificuldades. Neste tópico é estudada, avaliada e justificada a estratégia de pré-condicionamento para as dificuldades abordadas no tópico anterior, em especial para o subproblema o fluxo de carga. Esta estratégia explora tanto as características estruturais como também as propriedades espectrais da matriz Jacobiana, substituindo o sistema linear original por outro equivalente, com a matriz de coeficientes apresentando um menor número de condicionamento e autovalores mais afastados da origem do plano complexo, agrupados próximos à unidade. 3.2 O Pré-condicionador O termo ―pré-condicionamento‖ foi usado pela primeira vez para se referir às transformações realizadas em um sistema linear para reduzir o número de condicionamento da matriz de coeficientes. Posteriormente, foi visto que para melhorar a convergência de um método iterativo é necessário que a matriz de coeficientes apresente o menor número de condicionamento possível. Caso contrário, é aconselhável que o sistema seja preparado ou ―pré-condicionado‖ antes de ser solucionado. Portanto, inicialmente a estratégia de pré-condicionamento foi indicada para reduzir o número de condicionamento da matriz de coeficientes e consequentemente aumentando a taxa de convergência do método iterativo melhorando sua eficiência computacional. 56 57 A significativa melhoria do desempenho dos métodos pré-condicionados motivou o desenvolvimento de diversos algoritmos, estes agrupados em três categorias: précondicionadores de matriz descomposta, pré-condicionadores de aproximação esparsa da matriz inversa e, pré-condicionadores multinível (ou multigrid). As duas primeiras categorias já foram aplicadas em problemas de engenharia elétrica, especificamente para análises de fluxo de carga e no domínio do tempo. Já os pré-condicionadores multinível são apropriados para solucionar sistemas de equações lineares associados aos métodos de discretização do domínio dos Elementos Finitos e das Diferenças Finitas. Devido à extrema dificuldade da solução destes problemas, os algoritmos dos pré-condicionadores multinível, normalmente são desenvolvidos para implementação em ambiente computacional paralelo ou em arquiteturas computacionais com memória distribuída. Pré-condicionador de fatoração incompleta tem sido a técnica mais estudada e explorada entre os pré-condicionadores do tipo matriz descomposta, apresentando em diferentes aplicações os melhores desempenhos. Entre as estratégias de pré-condicionamento é dada maior ênfase aos pré-condicionadores de fatoração incompleta por apresentarem características que se adequam bem à solução do subproblema linear de fluxo de carga. 3.3 Formas de Pré-condicionamento O pré-condicionamento de um sistema linear é realizado usando-se um operador M conhecido como matriz de pré-condicionamento, de transformação, ou simplesmente précondicionador. No subproblema do fluxo de carga, para se construir um précondicionador deve-se decidir entre encontrar M como uma aproximação de A (précondicionadores de matriz descomposta) ou encontrar M como uma aproximação de A -1 (pré-condicionadores de aproximação esparsa da matriz inversa). Em ambos os casos, o sistema pré-condicionado é mais fácil de resolver do que o sistema original. Existem três formas de utilizar um pré-condicionador M, sendo estas: pré- condicionamento pelo lado esquerdo, pré-condicionamento pelo lado direito e précondicionamento por ambos os lados. As três formas produzem o mesmo sistema equivalente, com as mesmas propriedades espectrais e normalmente desempenhos similares do método iterativo. Porém, o pré-condicionamento pelo lado esquerdo é a forma mais simples e de fácil implementação computacional. Abaixo são apresentadas as três formas de se realizar o pré-condicionamento, para ambos os tipos: aproximação esparsa da matriz inversa e matriz descomposta. 57 58 Pré-condicionadores de aproximação esparsa da matriz inversa: M ≈ A-1 I. (Lado esquerdo) M·A·x = M·b II. (Lado direito) A·M·y = b, x = M·y III. (Ambos os lados) M2·A·M1·y = M2·b, x = M1·y Pré-condicionadores de matriz descomposta: M ≈ A I. (Lado esquerdo) M-1·A·x = M-1·b II. (Lado direito) A·M-1·y = b, x = M-1·y III. (Ambos os lados) M2-1·A·M1-1·y = M2-1·b, x = M1-1·y O pré-condicionador M possui duas partes fundamentais; algoritmo de pré- condicionamento e a regra de preenchimento. O algoritmo de pré-condicionamento é usado para construir M. Por exemplo, no caso dos pré-condicionadores baseados em fatoração incompleta, normalmente utiliza-se algum algoritmo de fatoração triangular (do tipo IKJ ou KIJ) para construir implicitamente M na forma de duas matrizes triangulares L’ e U’. Já a regra de preenchimento é o critério usado para decidir que elementos devem ser preenchidos ou não no pré-condicionador. As regras de preenchimento normalmente usam parâmetros escalares cujos valores são escolhidos antes de se iniciar o processo de solução. O significado de cada parâmetro depende do critério associado a cada regra de preenchimento. A ideia fundamental é escolher os valores dos parâmetros que permitam preencher o menor número de elementos sem prejudicar a qualidade do précondicionador, isto é, sem torná-lo muito diferente de A ou A-1. A grande quantidade de parâmetros associados aos métodos iterativos com pré-condicionadores é visto como uma desvantagem se comparados aos métodos diretos, já que estes não precisam de parâmetro nenhum. 3.4 Efeitos do Pré-condicionamento no Subespaço de Krylov do Método GMRES O método GMRES usa o algoritmo de Arnoldi para calcular os vetores base V que definem o subespaço de Krylov (vetores coluna vfk de V) e a técnica de rotações de Givens para calcular o vetor coluna y k. Tanto os vetores base como o vetor y k são usados para estimar em cada iteração k uma solução aproximada xk, como mostrado em (3.1). As informações numéricas para este procedimento são adquiridas da matriz de coeficientes A e do vetor do lado direito b. No entanto, quando o sistema é précondicionado, as informações são adquiridas da matriz de coeficientes pré-condicionada 58 59 M-1·A e do vetor do lado direito pré-condicionado M-1·b. Consequentemente, são gerados vetores base pré-condicionados VM e um vetor yMk pré-condicionado, cuja solução aproximada xMk, calculada de acordo com (3.2) está muito mais próxima da solução exata que a conseguida pelo sistema não pré-condicionado. xk = x0 + V  yk =  v1f   ε1    f f   ε2  v 2 ... v k  = ε1   v1f  + ε 2   v 2f  + ... + εk   v kf  (3.1)            V εm  yk M M M  Mf  M  Mf  M  Mf  xM k = x0 + V  yk = ε1   v1  + ε 2   v 2  + ... + εk   vk  (3.2) Ilustra-se através de um exemplo, como o pré-condicionamento influencia positivamente na criação de um subespaço de Krylov permitindo encontrar uma melhor aproximação da solução, consequentemente melhorar a convergência do método GMRES. O sistema linear (3.3) é solucionado de duas formas diferentes; a primeira usando-se o GMRES sem pré-condicionador e a segunda o GMRES com um pré-condicionador de fatoração incompleta apresentado em (3.4). Em (3.5) apresenta-se o subespaço de Krylov para o caso sem pré-condicionador (V e yk) e em (3.6) para o caso com pré-condicionador (VM e yMk), ambos calculados pelo algoritmo de Arnoldi e a técnica de Rotações de Givens do GMRES. 10 30 40  500   20 10 20   x  290         40 30 50  680  A (3.3) b 1  10 30 40   L '   1  , U'   50 60    1 2 (3.4) 0,56 -0,83 0,04  V  0,32 0,17 -0,93 0,76 0,53 0,37  (3.5) 11,9 yk  10,2 17,9 59 60  0,36 -0,93 0,01  V   0,71 0,28 0,64  -0,60 -0,23 0,77  M yM k -20,8    -6,9   9,0  (3.6) Na Figura 3.1(a) apresenta-se, para a primeira iteração do método GMRES, os subespaços de Krylov (retas), para os casos: sem pré-condicionamento (reta cinza) e com pré-condicionamento (reta preta). Nesta iteração, o erro para o caso sem précondicionamento é ||X*-X1|| = 20,5 e para o caso com pré-condicionamento é ||X*-XM1|| = 11,3. Na Figura 3.1(b) apresenta-se, para a segunda iteração do método GMRES, os subespaços de Krylov (planos). Nesta iteração, o erro para o caso sem précondicionamento é ||X*-X2|| = 17,8 e para o caso com pré-condicionamento é ||X*-XM2|| = 8,9. Observa-se que, em ambos os gráficos, o GMRES quando pré-condicionado aproximou-se mais da solução exata X*=[-1, -11, 21]T, apresentando os menores erros. a) Primeira iteração do GMRES b) Segunda iteração do GMRES Figura 3.1 – Subespaços de Krylov, suas soluções aproximadas e erros para cada iteração do GMRES, com e sem estratégia de pré-condicionamento. 3.5 Efeitos da Estratégia de Pré-condicionamento nas Propriedades Espectrais da Matriz Jacobiana As propriedades espectrais são úteis para avaliar os efeitos do pré-condicionamento no condicionamento e singularidade dos sistemas de equações lineares típicos de fluxo de carga. 60 61 São comparados os autovalores e número de condicionamento das matrizes Jacobianas de sistemas lineares pré-condicionados e não pré-condicionados. As matrizes Jacobianas estão associadas aos sistemas lineares da primeira iteração Newton-Raphson do problema de fluxo de carga no ponto mais Próximo do Máximo Carregamento (PMC), apresentados no tópico anterior para os sistemas hipotéticos do IEEE e para as configurações dos sistemas norte-nordeste e brasileiro. Na Figura 3.2 apresentam-se os gráficos dos autovalores das matrizes Jacobianas não pré-condicionadas (em vermelho) e pré-condicionadas (em azul). Observa-se que, em todos os gráficos os autovalores das matrizes Jacobianas pré-condicionadas estão agrupados próximos d a unidade e grande parte afastada da origem do plano complexo. Esta característica indica que os sistemas lineares pré-condicionados estão mais afastados da singularidade. Observa-se também que em todos os casos os autovalores possuem parte real sempre positiva, ou seja, os sistemas lineares que antes eram indefinidos (IEEE300-barras e Norte-Nordeste de 274 barras), agora são definidos. Nestes casos, os métodos iterativos que antes podiam apresentar um processo lento de convergência ou até mesmo falhar, agora devem apresentar uma rápida convergência livre de falhas. Todas estas informações podem também ser conferidas na Tabela 3.1, onde são apresentados os autovalores com menor e maior parte real, tanto para o caso não pré-condicionado como também para o caso pré-condicionado. Com a finalidade de comparar os resultados anteriores, na Figura 3.3 apresentam-se em escala logarítmica os valores do número de condicionamento das matrizes Jacobiana de todos os sistemas considerados aqui, para os casos com e sem pré-condicionamento. Observe-se que sempre o número de condicionamento é menor nos sistemas lineares com pré-condicionamento, comprovando que a estratégia de pré-condicionamento melhora significativamente as propriedades espectrais da matriz Jacobiana, substituindo o sistema linear original por outro sistema linear menos singular e não indefinido. Portanto, espera-se que o desempenho do método GMRES melhore significativamente durante a solução destes sistemas lineares pré-condicionados. 61 62 a) Sistema IEEE30-barras b) Sistema IEEE118-barras c) Sistema IEEE145-barras d) Sistema IEEE162-barras e) Sistema IEEE300-barras f) Sistema Norte-Nordeste de 274 barras Figura 3.2 – Autovalores da matriz Jacobiana pré-condicionada (azul) e não pré-condicionada (vermelho) de sistemas lineares de pequeno porte, no caso base e no ponto mais próximo do máximo carregamento (PMC). 62 63 Tabela 3.1 – Propriedades espectrais das matrizes Jacobianas Pré-condicionadas e Não Précondicionadas. Sistema Autovalores do Sistema Précondicionado Autovalores do Sistema Não Précondicionado Com a Menor Parte Com a Maior Parte Com a Menor Parte Com a Maior Parte Real Real Real Real IEEE30 0,19±0,13i 1,12±0,16i 0,00798 95,70±26,42i IEEE118 0,00924 1,14 0,00983 360,60±164,28i IEEE145 0,03460 1,24±0,01i 0,00519 5.621,68±202,72i IEEE162 0,00741 1,21 0,01871 2.122,98±209,50i IEEE300 0,00014 1,37 -1,45015 4.164,81±673,07i Norte-Nordeste 274 0,00187 1,20 -1.873,08000 20.703,42 Figura 3.3 – Número de condicionamento para o caso próximo do máximo carregamento (PMC) 3.6 Pré-condicionadores de Fatoração Incompleta Os pré-condicionadores de fatoração incompleta são aproximações da matriz de coeficientes, construídos na forma de fatores triangulares L’ e U’, sendo que, o preenchimento ou não de cada elemento dos fatores triangulares deve estar de acordo com uma ―regra de preenchimento‖. Consequentemente, os fatores triangulares são calculados de forma incompleta e esparsa, sendo que seu produto é igual à matriz de pré-condicionamento M, como apresentado em (3.7). Normalmente, o processo de 63 64 construção deste tipo de pré-condicionador é realizado utilizando-se algum dos diferentes tipos de algoritmos de eliminação de Gauss comumente usados nos métodos diretos, alguns identificados segundo a ordem dos três laços associados com o controle das variáveis i, j e k, usada nos algoritmos. M  L' U' (3.7) 3.7 Regras de Preenchimento em Pré-condicionadores de Fatoração Incompleta Quando usado um método direto, durante a fatoração LU calculam-se todos os elementos para encontrar a solução exata do sistema linear após substituição direta/inversa. Já um pré-condicionador de fatoração incompleta não precisa calcular todos os elementos dos fatores triangulares para que o método iterativo encontre a solução. Para reduzir o esforço computacional realizado durante a construção do précondicionador, são considerados apenas alguns elementos nos fatores triangulares. Este procedimento é conhecido como preenchimento de elementos não-nulos e produz um pré-condicionador M de fatores triangulares incompletos (L’ e U’). O procedimento é realizado utilizando-se algumas heurísticas conhecidas como regras de preenchimento para decidir quais elementos dos fatores triangulares devem ser preenchidos e quais devem ser descartados (substituídos por zero). A regra pode ser aplicada para preencher elementos em L’ e U’, tanto em posições ocupadas como também em posições antes não ocupadas da matriz de coeficientes, neste último caso, está-se referindo ao preenchimento ou não de novos elementos não-nulos. Observa-se que, por cada elemento descartado, evita-se realizar todas as operações matemáticas futuras onde este elemento poderia ser requerido, tanto durante a construção do pré-condicionador como também cada vez que é requerido no processo iterativo do GMRES, melhorando a eficiência computacional da sua construção. Simultaneamente, por cada elemento descartado se introduz um erro que pode chegar a ser considerável ou não, dependendo principalmente da quantidade de vezes que o elemento descartado seria usado nas próximas operações, e também dos valores numéricos de todos os outros elementos envolvidos em tais operações. Depois de terminada a fatoração incompleta, uma forma simples de estimar o erro total originado pelos elementos descartados é através de (3.8). O produto dos fatores triangulares incompletos é igual ao pré-condicionador M, e o produto dos fatores triangulares completos é igual à matriz de coeficientes A (3.9). Quando não usada regra de preenchimento, os fatores triangulares são calculados completamente e o pré64 65 condicionador é igual à matriz de coeficientes, e o GMRES converge em uma única iteração. R  A M (3.8) R  L  U  L ' U' (3.9) Portanto, uma regra de preenchimento eficiente deveria descartar o maior número possível de elementos, sem incrementar o erro até valores que afetem negativamente ou impeçam a convergência do método GMRES. A seguir ilustra-se como uma regra de preenchimento quando muito simplificada introduz muitos erros e pode ser inapropriada para problemas que requerem regras mais sofisticadas. 3.8 Pré-condicionador ILU(0) Uma das regras de preenchimento mais simples foi adotada pelo pré-condicionador ILU(0), permitindo elementos não-nulos apenas em posições antes ocupadas por elementos não-nulos da matriz de coeficientes. Este tipo de pré-condicionador pode ser útil para algumas matrizes do tipo banda de pequeno porte, associadas à discretização de equações diferencias parciais por diferenças finitas. Na Figura 3.4 apresentam-se os fatores triangulares incompletos L’ e U’ em (a) e (b) respectivamente, a matriz de coeficientes A em (c) e finalmente o produto L’·U’ em (d). Observa-se que o pré-condicionador ILU(0) criou fatores triangulares com a mesma estrutura que a matriz de coeficiente, no entanto, o produto L’ U’ possui uma estrutura diferente à estrutura da matriz de coeficientes e 51 elementos de valores diferentes (de cor cinza), sendo que estes elementos são os responsáveis pelo aumento do erro. Neste exemplo estas diferenças não chegam a prejudicar a convergência do GMRES e este tipo de pré-condicionador pode ser muito apropriado para este tipo de sistema linear. Na Figura 3.5 apresenta-se outro conjunto de matrizes (Jacobiana do problema de fluxo de carga - configuração do sistema elétrico brasileiro de 3.513 barras) com os fatores triangulares incompletos L’ e U’ apresentados em (a) e (b) respectivamente, a matriz de coeficientes A em (c) e o produto L’·U’ em (d). O pré-condicionador ILU(0) criou fatores triangulares com a mesma estrutura que a matriz de coeficientes, mas o produto L’ U’ possui muitos elementos de valores diferentes aos encontrados na matriz de coeficientes (elementos de cor cinza), sendo os responsáveis pelo aumento do erro. Neste caso, 65 66 estas diferenças podem prejudicar bastante a convergência do GMRES, sendo, portanto, inapropriado para este tipo de problema de dimensões e dados reais. Nos gráficos nz é o número total de elementos não-nulos. a) Matriz triangular inferior L’ c) Matriz de coeficientes A b) Matriz triangular superior U’ d) Produto L’·U’ (Pré-condicionador M) Figura 3.4 – Fatoração incompleta ILU(0) para uma matriz banda de pequeno porte. Normalmente, em problemas reais envolvendo matrizes de certas dimensões, é necessário preencher mais posições nos fatores triangulares do que os permitidos no ILU(0), pois este pode ficar muito distante dos fatores triangulares completos. Portanto, para estes casos são necessárias regras de preenchimento mais complexas e sofisticadas que permitam adicionar novos elementos não-nulos aos fatores triangulares. No entanto, a escolha de quais novos elementos não-nulos devem ser adicionados pela regra de preenchimento deve ser realizada de forma cuidadosa, a fim de preencher apenas aqueles elementos que vão diminuir significativamente o erro. Uma eventual desconsideração destes elementos poderia manter os erros introduzidos pela aproximação, e poderiam continuar prejudicando ou impedindo a convergência do método iterativo. Este problema está associado às deficiências das heurísticas das regras de preenchimento dos pré-condicionadores convencionais. 66 67 a) Matriz triangular inferior L’ b) Matriz triangular superior U’ c) Matriz de coeficientes A d) Produto L’·U’ (Pré-condicionador M) Figura 3.5 – Fatoração incompleta ILU(0) para a matriz Jacobiana do sistema brasileiro de 3.513 barras. 3.9 Pré-condicionador ILU(k) Fatorações triangulares incompletas mais precisas são geralmente mais eficientes e confiáveis por se aproximarem mais da matriz de coeficientes, e por produzirem uma taxa de convergência mais adequada no GMRES. Este tipo de pré-condicionador se diferencia da classe ILU(0) por permitir alguns preenchimentos de novos elementos não-nulos em posições antes nulas. Este é o caso da classe de pré-condicionadores ILU(k), onde k especifica qual o nível de preenchimento usado pelo pré-condicionador. 67 68 A cada elemento que é processado através da eliminação de Gauss é associado um nível de preenchimento, sendo que, dependendo do valor desse nível, o elemento pode ser preenchido ou descartado. O nível de preenchimento indica de forma aproximada o tamanho dos elementos; quanto maior o nível, menor a magnitude dos elementos. Os aspectos negativos associados ao pré-condicionador ILU(k) são os seguintes: 1. Devido a que o número de novos elementos não-nulos preenchidos depende do valor de k, então, para k > 0, não é possível estimar a quantidade de elementos não-nulos e o esforço computacional necessário para se obter a fatoração ILU(k). 2. O custo associado à atualização dos níveis de preenchimento de elementos nãonulos pode ser elevado. 3. O nível de preenchimento de elementos não-nulos para matrizes indefinidas pode não ser um bom indicador da magnitude dos elementos que estão sendo descartados e o algoritmo pode descartar elementos grandes resultando numa imprecisão da fatoração incompleta, no sentido que ||R|| = ||L·U-A|| pode não ser pequeno. A prática tem mostrado que geralmente isto leva a um grande número de iterações para se obter a convergência. 4. Não existe um método que permita escolher apropriadamente o valor de k, a escolha de k depende do problema que está sendo solucionado, normalmente k é escolhido de forma empírica para apenas um sistema linear específico. Em problemas que precisam solucionar vários sistemas lineares sequencialmente, como é o caso do problema de fluxo de carga é necessário encontrar um valor ou faixa de valores adequados para k tal que garantem a solução de todos os sistemas lineares no menor tempo de execução possível. Em análises mais complexas como, por exemplo, aquelas que precisam resolver vários fluxos de carga de diferentes pontos de operação, a escolha de um valor apropriado para k se torna ainda mais relevante. Estes aspectos negativos sugerem que o pré-condicionador ILU(k) com uma regra baseada nos níveis de preenchimento, pode não ser o mais indicado para análises de sistemas de potência como, por exemplo, o fluxo de carga. A seguir, apresenta-se outro pré-condicionador de fatoração incompleta com uma regra de preenchimento diferente, que considera os valores numéricos de cada elemento da matriz de coeficientes. 68 69 3.10 Pré-condicionador ILUT(,ρ) Pré-condicionadores cujas regras são baseadas apenas no nível de preenchimento consideram apenas a posição dos elementos (estrutura da matriz de coeficientes) sem levar em consideração o valor numérico dos elementos dos fatores triangulares. Portanto, este tipo de pré-condicionador pode provocar a perda da robustez e reduzir a velocidade de convergência do GMRES, causando problemas em aplicações práticas e reais. Outra estratégia de eliminação de elementos não-nulos considera a magnitude dos elementos ao invés de apenas a sua posição. Neste caso, a regra de preenchimento usa dois parâmetros ( e ), onde,  é um número inteiro que limita o número máximo de elementos permitidos em cada linha dos fatores triangulares L e U. Já  é um número real (tolerância) usado para descartar os elementos considerados pequenos. Basicamente,  é responsável pelo controle do uso da memória, enquanto que  é usado para reduzir o tempo computacional. O pré-condicionador associado a esta regra de eliminação é conhecido como ILUT(,ρ). A ideia central que permite definir a regra de preenchimento do pré-condicionador ILUT(,ρ) é descartar apenas os menores elementos (em valor absoluto) dos fatores triangulares. Neste caso, aparentemente, quando os menores elementos são descartados se produz um menor erro do que quando se descarta os maiores. A regra se fundamenta em duas heurísticas; uma usa o parâmetro  para decidir quais elementos pequenos devem ser descartados; e a outra usa  para decidir quantos elementos grandes serão preenchidos em cada linha dos fatores triangulares. Um elemento wk é substituído por zero se este for menor que uma tolerância relativa i obtida multiplicando-se a tolerância  pela norma da i-ésima linha da matriz de coeficientes original (norma-2). Entretanto, são mantidos apenas os  maiores elementos de cada linha no fator L e no fator U, em adição aos elementos diagonais, que são sempre mantidos. 3.11 Pré-condicionadores Aplicados a Problemas de Sistemas Elétricos de Potência Neste ponto, fica mais evidente a importância de se construir um pré-condicionador de qualidade para melhorar o desempenho do método iterativo e torná-lo competitivo em relação aos tradicionais métodos diretos, tanto em termos de eficiência computacional 69 70 como também de robustez. A seguir são analisados os principais avanços ocorridos durante os últimos anos sobre desenvolvimentos computacionais realizados em précondicionadores. 3.12 Características dos Principais Pré-condicionadores A Tabela 3.2 apresenta as características mais relevantes das estratégias de précondicionamento aplicadas a problemas de sistemas elétricos de potência. A primeira coluna identifica o pré-condicionador e a segunda indica a categoria a qual pertence cada estratégia de pré-condicionamento, ou seja, pré-condicionadores de aproximação esparsa da matriz inversa (PAEI) e pré-condicionadores de matriz descomposta (PMD). Na terceira coluna é informado o tipo de matriz de coeficientes que cada algoritmo de pré-condicionamento utiliza para calcular o pré-condicionador M. Na quinta coluna é mencionado o tipo de algoritmo usado para construir o pré-condicionador. Finalmente na última coluna apresenta-se o tipo de regra de preenchimento dos elementos do précondicionador. A primeira estratégia de pré-condicionamento, que se tem conhecimento, aplicada ao problema de fluxo de carga foi Cholesky incompleto utilizando-se o método iterativo Gradiente Conjugado. Para construir este tipo de pré-condicionador é necessário que a matriz de coeficientes seja numericamente simétrica e definida positiva. Este précondicionador é do tipo matriz descomposta devido à decomposição de Cholesky em fatores triangulares. A regra de preenchimento não permite a presença de novos elementos não-nulos nos fatores incompletos, isto é, apenas são preenchidas as posições antes ocupadas na matriz de coeficientes original. Por esta razão, esta regra não usa nenhum tipo de parâmetros. 70 71 Tabela 3.2 – Principais características dos pré-condicionadores propostos na área de sistemas de potência Nome do Pré-condicionador Tipo Tipo de Matriz de Coeficientes Cholesky Incompleto PMD Simétrica e definida positiva: B’ e B‖ Jacobiana da Primeira Iteração PMD Assimétrica: J Fatoração LU: 0 0 0 M=J =L ∙U --- Cholesky Completo da Primeira Iteração PMD Simétrica e definida positiva: B’ e B‖ Fatoração Cholesky: T M=L∙L --- Matriz Diagonal PAEI Assimétrica: J M=Diagonal(J) Preencher elementos da diagonal ILU(0) PMD Assimétrica: J Algoritmo KJI de Fatoração LU: M=L’∙U’ Sem novos elementos não-nulos Blocos Diagonais PMD Assimétrica: Blocos Pθ e QV de J FastD PMD Assimétrica: B’ --- ILU() Algoritmo KJI de Fatoração LU: M=B’=L∙U PMD Assimétrica: J Fatoração LU: M=L’∙U’ Com parâmetro limitante  Aproximação da PAEI Inversa por Chebyshev PAEI Assimétrica: J ILU(k) PMD Assimétrica: J ILUT() PMD Assimétrica: J PMD Assimétrica: J PMD Assimétrica: J Multifrontal LU Fatoração Cholesky: Sem novos elementos T M=L’∙L’ não-nulos Algoritmo KJI de Sem novos elementos Fatoração LU: não-nulos M=Blocos (Pθ e QV) Aproximação de B’-1 Simétrica e definida usando o Polinômio de positiva: B’ e B‖ Chebyshev JFNG ILUT(,p) Quase-Fixo Regra de Preenchimento Algoritmo Polinômio de Chebyshev Fórmulas de Broyden Fórmulas de Broyden Algoritmo IKJ de Fatoração LU: M=L’∙U’ Nível de novos elementos não-nulos Algoritmo IKJ de Fatoração LU: M=L’∙U’ Algoritmo IKJ de Fatoração LU: M=L’∙U’ Com parâmetro limitante  Com duplo parâmetro limitante  e p Algoritmo Multifrontal de Fatoração LU: M=J0=L0∙U0 --- Uma forma simples e prática de pré-condicionamento é usar os fatores triangulares de Cholesky completos, sem regra de preenchimento, consequentemente também sem parâmetros. Este tipo de pré-condicionador já foi para análise de contingências, usandose sempre o mesmo pré-condicionador durante todo o processo de solução, evitando-se uma nova construção por cada sistema linear. Uma versão deste tipo de précondicionador para matrizes assimétricas de fluxo de carga foi apresentada anteriormente, onde foi proposto usar a matriz Jacobiana da primeira iteração Newton71 72 Raphson para construir um pré-condicionador igual ao produto dos fatores triangulares completos, sem regra de preenchimento, sem parâmetros e fixo durante toda a simulação. No pré-condicionador conhecido como ILU(0) sua regra de preenchimento permite apenas usar as posições antes ocupadas na matriz de coeficientes original, não permitindo o preenchimento de novos elementos não-nulos. O pré-condicionador de Blocos Diagonais da matriz Jacobiana (J) também não oferece dificuldades na sua construção, pois considera apenas os blocos de Pθ e QV de J. O pré-condicionador é construído implicitamente na forma de fatores triangulares incompletos sem considerar novos elementos não-nulos. Já o pré-condicionador FastD é igual à matriz B’ do método desacoplado rápido, construído na forma de fatores triangulares completos. Outra estratégia conhecida como pré-condicionador de fatoração incompleta com parâmetro limitante ILU() é igual ao produto dos fatores triangulares incompletos, construído usando-se uma regra de preenchimento com parâmetro limitante . O pré-condicionador Newton-GMRES Livre de Jacobiana (JFNG) pode ser considerado do tipo aproximação esparsa da matriz inversa, já que em cada atualização preenche cada vez mais elementos para se aproximar mais da matriz inversa. A Tabela 3.3 apresenta as principais vantagens e desvantagens dos pré-condicionadores apresentados na Tabela 3.2. 72 Tabela 3.3 – Principais vantagens e desvantagens dos pré-condicionadores aplicados a problemas de sistemas de sistemas de potência. Précondicionador Cholesky Incompleto Jacobiana da Primeira Iteração Matriz Diagonal ILU(0) Blocos Diagonais Vantagens  Rápida construção devido a que apenas precisa calcular e armazenar elementos de L’ (U=L’T).  É construído só uma vez (pré-condicionador fixo).  Isento de erros introduzidos por regras de eliminação.  Custo computacional insignificante para sua construção.  Baixo custo computacional durante sua construção.  Baixo custo computacional durante sua construção. FastD  É construído só uma vez (pré-condicionador fixo). ILU()  Produz um desempenho aceitável do método iterativo. ILU(k)  Produz bom desempenho do método iterativo. Aproximação da Inversa por Chebyshev JFNG ILUT() ILUT(,p) Quase-fixo Multifrontal LU  Nenhuma.  Nenhuma.  Produz o melhor desempenho registrado em métodos iterativos aplicados ao problema de fluxo de carga.  Melhor robustez que os pré-condicionadores fixos.  Converge em fluxo de carga com controles e limites.  Produz o melhor desempenho registrado em métodos iterativos aplicados ao problema de simulação dinâmica. Desvantagens  Pode ser usado apenas com métodos desacoplados.  Não há convergência em fluxo de carga com controles e limites.  Produz péssimo desempenho do método iterativo em sistema de médio e maior-porte.  Produz baixo desempenho do método iterativo em sistema de médio e maior-porte.  Produz baixo desempenho do método iterativo em sistema de médio e maior-porte.  Erros introduzidos devido a que é usada B’ invés de J para costruir o pré-condicionador.  Não há convergência em fluxo de carga com controles e limites.  Péssimo critério para eleger o valor do parâmetro limitante .  A regra de preenchimento pode inserir erros quando a matriz J seja mal escalonada.  A regra de preenchimento não tem controle do erro inserido quando um elemento não é preenchido.  Apenas deixa de preencher novos elementos não-nulos.  Pode ser usado apenas com métodos desacoplados.  Custo computacional alto para a construção do pré-condicionador.  Custo computacional alto para a construção do pré-condicionador, precisa de várias atualizações para melhorar sua qualidade.  A regra de preenchimento não tem controle do erro inserido quando um elemento não é preenchido..  A regra de preenchimento não tem controle do erro inserido quando um elemento não é preenchido.  Não há convergência em fluxo de carga com controles e limites. 48 73 74 Pré-condicionadores baseados no cálculo da aproximação da inversa não devem ser aplicados no problema de fluxo de carga, já que o custo computacional associado à sua construção é elevado. Os pré-condicionadores baseados no cálculo dos fatores triangulares incompletos ILU(k) e ILUT são atualmente os mais eficientes e apropriados para o problema de fluxo de carga. Simulações apresentadas nesta apostila confirmam sua influência positiva no desempenho dos métodos iterativos, não apenas em termos de eficiência computacional, mas também em termos de robustez. 3.13 Avaliação do Desempenho do Método GMRES O objetivo principal é verificar qual a influência da estratégia de pré-condicionamento sobre o desempenho do método GMRES sem reordenamento. Para tal, são realizados experimentos numéricos nos sistemas-teste IEEE30, IEEE118, IEEE145, IEEE162, IEEE300 e o sistema Norte-Nordeste de 274 barras para comparar o desempenho de quatro configurações que agregam o método GMRES com quatro tipos de précondicionadores; ILU(k) e ILUT(,ρ) do tipo matriz descomposta, SPAI e AINV(,ρ) do tipo aproximação esparsa da inversa. As simulações são realizadas unicamente para o cenário do ponto mais Próximo do Máximo Carregamento (PMC) e os objetivos aqui são, comprovar quantitativamente os benefícios que o pré-condicionamento produz no desempenho do método GMRES, comparar os pré-condicionadores ILU(k) e ILUT(,ρ), e comprovar que a qualidade dos pré-condicionadores de aproximação da inversa é inferior a dos pré-condicionadores de fatoração incompleta. Na Figura 3.6 apresenta-se os resultados das simulações em 6 gráficos (para cada sistema elétrico). Em cada gráfico apresenta-se 4 grupos de diagramas de barras, cada um apresentando 3 barras onde a barra em preto informa o número de operações realizadas pelo GMRES, a cinza escura para a construção do pré-condicionador e a cinza clara é o número total de operações realizadas (soma das duas barras anteriores). Baseado no número de operações de ponto flutuante, o GMRES+ILUT(,ρ) apresentou o melhor desempenho em comparação aos demais, e por esta razão será utilizado nos próximos experimentos para comparações com o pré-condicionador proposto nesta apostila. O desempenho dos pré-condicionadores de aproximação esparsa da inversa é inferior já que realizam muitas operações de ponto flutuante durante o preenchimento de muitos elementos da matriz inversa, como visto principalmente em (b), (c) e (f). 74 75 Na Tabela 3.4 são apresentados os resultados dos desempenhos do GMRES+ILUT() aplicado ao problema de fluxo de carga sem o uso do parâmetro limitante ρ e o précondicionador ILUT(,ρ) aplicado para o problema de interesse, lembrando que a única diferença entre ambos os pré-condicionadores é que este último incorpora o parâmetro ρ na regra de preenchimento. Os resultados mostram que na maioria dos casos um valor pequeno para ρ é capaz de reduzir bastante o número de operações de ponto flutuante uma vez que se limita o número de elementos por linha nas matrizes L e U. 75 76 a) Sistema IEEE30-barras b) Sistema IEEE118-barras c) Sistema IEEE145-barras d) Sistema IEEE162-barras e) Sistema IEEE300 barras f) Sistema Norte-Nordeste de 274 barras Figura 3.6 – Desempenho computacional do método GMRES com pré-condicionadores de fatoração incompleta e de aproximação da inversa. 76 77 Por exemplo, para o ILUT() se uma linha tem 20 elementos cujos valores são maiores que , então todos são preenchidos. No entanto, para o ILUT(,ρ), se o valor de ρ for 10, mesmo que os 20 elementos sejam maiores que , apenas são preenchidos os 10 maiores. Neste caso, são evitadas todas as operações associadas com os outros 10 elementos não preenchidos. Um caso excepcional é apresentado no sistema IEEE300barras, já que o uso do parâmetro ρ não fez diferença. Tabela 3.4 – Comparação de desempenhos entre o GMRES com pré-condicionador ILUT(,ρ) e com ILUT() GMRES+ILUT(,ρ) Sistema Teste  GMRES+ILUT() ρ MFLOPS  -2 10 0,0141 10 -3 MFLOPS -2 0,0162 10 -3 0,3973 0,1464 10-4 1,2065 10 0,4432 10 -3 0,6976 -2 32 1,7386 10 -2 1,7386 -4 20 1,5843 10 -4 3,1925 IEEE30-barras 10 IEEE118-barras 10 10 0,2498 IEEE145-barras 10-4 9 IEEE162-barras 10 -3 IEEE300-barras 10 Norte-Nordeste de 274 barras 10 3.14 Reordenamento da Matriz Jacobiana São apresentados experimentos numéricos para avaliar o desempenho do método GMRES quando usadas conjuntamente estratégias de reordenamento e pré- condicionamento. O reordenamento melhora a qualidade do pré-condicionador e reduz o tempo necessário para sua construção, gerando subespaços de Krylov mais próximos da solução exata que quando não usado, chegando à convergência em poucas iterações. Portanto, a influência do reordenamento na convergência do método iterativo se dá de forma indireta, através do pré-condicionador. Sem o pré-condicionador o reordenamento não tem efeito nenhum sobre o subespaço de Krylov, como mostrado em alguns experimentos realizados. O objetivo principal é verificar qual a influência das estratégias de reordenamento MD e RCM sobre o desempenho do método GMRES com pré-condicionamento. Para tal, verifica-se a eficiência computacional e a robustez de quatro configurações, sendo estas: GMRES puro, GMRES+ILUT(,ρ), GMRES+ILUT(,ρ) e reordenamento RCM e 77 78 GMRES+ILUT(,ρ) e reordenamento MD. No teste de convergência do método GMRES considera-se que tanto a tolerância absoluta como a relativa são iguais a 10 -5. Os mesmos sistemas do experimento anterior são mantidos aqui, sendo que agora duas condições de operação são consideradas; uma para carga normal (base) e a outra para carga pesada (PMC). Na Tabela 3.5 apresenta-se o número de operações de ponto flutuante para ambos os casos referentes a cada técnica de reordenamento. Adicionalmente, para fins comparativos, na penúltima coluna apresenta-se o número de operações realizadas pelo solucionador direto MA28. Comparando a terceira e a quarta coluna desta tabela, observa-se que a estratégia de pré-condicionamento melhora consideravelmente a eficiência do método iterativo. No entanto, seu desempenho é inferior ao direto devido às operações realizadas pelo grande número de novos elementos não-nulos gerados no pré-condicionador. Portanto, já que o pré-condicionador é baseado em fatoração triangular é imprescindível considerar alguma estratégia de reordenamento que permita reduzir o número de novos elementos não-nulos nos fatores triangulares e, consequentemente reduzir o número de operações realizadas. O desempenho do método iterativo quando se usa reordenamento (RCM e MD) são apresentados na quinta e sexta coluna. Observa-se que para ambos os reordenamentos o método iterativo pré-condicionado foi superior ao direto, destacando-se o reordenamento MD. Para o sistema IEEE162, o MA28 apresentou uma interrupção durante a substituição direta/inversa. No entanto, foi solucionado normalmente pelo programa ANAREDE para ambos os casos (Base e PMC). 78 79 Tabela 3.5 – Número de operações de ponto flutuante para as 4 configurações em (MFLOPS), comparação de desempenho com o solucionador direto MA28. Sistema IEEE30barras IEEE118barras IEEE145barras IEEE162barras IEEE300barras Norte-Nordeste de 274 barras GMRES+ GMRES+ GMRES+ ILUT(,ρ)+ ILUT(,ρ)+ ILUT(,ρ) RCM MD MA28 EficiênciaSI EficiênciaSD 0,0039 0,0055 1,4 0,0100 0,0077 0,0111 1,4 0,1397 0,0583 0,0424 0,0894 2,1 575,7488 0,2498 0,0782 0,0565 0,1194 2,1 C. Base 227,1236 0,1464 0,0763 0,0751 0,1893 2,5 PMC 278,2075 0,1821 0,0868 0,0841 0,2012 2,4 C. Base 294,2203 0,2665 0,1951 0,1792 Interrup. --- PMC 410,3758 0,4432 0,3263 0,3198 Interrup. --- 0,9901 0,1962 0,1019 0,1763 1,7 5.663,5704 1,7386 0,3246 0,1753 0,2645 1,5 C. Base 4.741,2796 1,0788 0,1996 0,1062 0,1448 1,4 1,5843 0,2660 0,1416 0,1930 1,4 Cenário GMRES C. Base 1,4986 0,0141 0,0050 PMC 3,0990 0,0212 C. Base 534,1593 PMC C. Base 4.803,7987 PMC PMC 5.278,0255 Os melhores resultados são conseguidos pelo GMRES+ILUT(,ρ) com reordenamento MD. Portanto, esta é a melhor configuração. Na última coluna da Tabela 3.5 apresenta-se a relação de eficiências entre a melhor configuração e o solucionador direto. Esta coluna informa a superioridade do GMRES+ILUT(,ρ) com reordenamento MD em relação ao solucionador direto. A relação de eficiência é calculada como a divisão do número de operações realizadas pelo solucionador direto entre o número de operações realizadas pelo GMRES+ILUT(,ρ) com reordenamento MD. Por exemplo, esta configuração foi 2,1 vezes mais rápido que o solucionador direto na simulação de fluxo de carga do sistemateste IEEE118 para ambos cenários. Na Figura 3.7 apresentam-se os valores de taxa de convergência do método GMRES puro (barra em vermelho), GMRES+ILUT(,ρ) (em azul), GMRES+ILUT(,ρ) com reordenamento RCM (barra em verde) e GMRES+ILUT(,ρ) com reordenamento MD (barra em amarelo), para cada sistema elétrico no caso-base. As barras verdes e amarelas associadas com a taxa de convergência do método GMRES usando reordenamento RCM e MD, respectivamente, mostram que nestes casos o GMRES deve satisfazer o teste de convergência em menos iterações que os outros casos onde não é usado reordenamento (barras vermelhas e azuis). 79 80 Figura 3.7 – Taxa de convergência do método GMRES para as quatro configurações dos solucionadores iterativos no CASO BASE, na primeira iteração Newton-Raphson. Na Figura 3.8 apresentam-se os valores de taxa de convergência para as mesmas configurações anteriores do método GMRES, para a condição de PMC. De forma similar ao caso base, observa-se que para todos os sistemas os valores das taxas de convergência com reordenamento são indicadores de uma rápida convergência do método GMRES. Com estes resultados comprova-se que o reordenamento não serve apenas para reduzir o número de novos elementos não-nulos, como normalmente é observado nos métodos diretos e nos pré-condicionadores de fatoração incompleta, mas também influência na convergência do método iterativo GMRES com pré- condicionamento. Observa-se nas Figuras 3.7 e 3.8 que o efeito do reordenamento para o sistema-teste IEEE162 é desprezível, uma vez que não foi possível obter um pré-condicionador de boa qualidade, ou seja, próximo da matriz Jacobiana. Isto reduziria o erro em cada iteração linear do GMRES, o número de iterações lineares, aumentando a taxa de convergência. Os efeitos do reordenamento podem ser observados para outros valores de  e ρ. 80 81 Figura 3.8 – Taxa de convergência do método GMRES para as quatro configurações dos solucionadores iterativos no PMC, na primeira iteração Newton-Raphson. A ocorrência de zeros ou número próximos de zero na diagonal principal da matriz Jacobiana pode dificultar a construção do pré-condicionador nos seguintes aspectos: (a) incrementa o uso de pivotamento durante a fatoração para evitar interrupções, atrasando o processo; (b) diminui a precisão do cálculo dos elementos dos fatores incompletos; (c) torna os vetores coluna da matriz de coeficientes mais linearmente dependentes; (d) origina problemas de estabilidade numérica (aumenta o erro de arredondamento) durante a fatoração. Na Tabela 3.6 são apresentados dados estatísticos dos elementos da diagonal principal da matriz Jacobiana dos sistemas hipotéticos do IEEE, sistema Norte-Nordeste e configurações do sistema brasileiro, divididos em dois grupos: elementos pertencentes à submatriz H e elementos pertencentes à submatriz L. 81 82 Tabela 3.6 – Dados estatísticos dos elementos da diagonal principal das matrizes Jacobianas de típicos sistemas lineares de fluxo de carga. Informações na Diagonal Principal de J Sistema Elétrico Número de Zeros Menor elemento em: H L IEEE30-barras 0 1,7971 5,6438 IEEE118-barras 0 1,2926 14,9978 IEEE145-barras 0 1,5395 2,4066 IEEE162-barras 0 3,5996 8,8783 IEEE300-barras 0 0,1630 8,4786 Norte-Nordeste de 274 barras 0 0,2428 8,9685 Brasil de 2.256 barras 0 0,3300 3,1835 Brasil de 3.513 barras 0 0,1651 2,0241 Observa-se que em nenhum dos experimentos apareceu um valor zero para provocar interrupções durante a fatoração triangular, nem os menores elementos da diagonal principal. Portanto, a construção do pré-condicionador para todos os casos está livre de interrupções e nenhuma estratégia que evite problemas de estabilidade numérica causada por números próximos de zero é necessária para a solução do problema de fluxo de carga dos sistemas apresentados aqui. 82 83 3.15 Evolução do Desempenho do GMRES Para ilustrar resumidamente a evolução das melhorias em termos computacionais conseguidas no desempenho do GMRES, na Figura 3.9 mostra-se o gráfico com o número total de operações em MFlops associadas aos experimentos usando o sistema IEEE118 para os casos: GMRES puro (sem estratégias), GMRES(m), GMRES+ILUT(τ,ρ), GMRES+ILUT(τ,ρ) com reordenamento MD. Estes últimos estudos justificam a importância do uso das técnicas de pré-condicionamento e reordenamento em métodos iterativos do subespaço Krylov. A condição de operação considerada é para carga normal (base). Figura 3.9 – Evolução das melhorias conseguidas no desempenho do GMRES, comparação com o solucionado Direto A melhor alternativa conseguida até aqui é composta por GMRES+ILUT(τ,ρ) com reordenamento MD. Deve-se destacar que o pré-condicionador ILUT(τ,ρ) com duplo parâmetro limitante ainda não tinha sido aplicado ou avaliado no problema de fluxo de carga por outros pesquisadores, o quê leva a concluir que: • O reordenamento aplicado em conjunto com a estratégia de pré-condicionamento é determinante para que o método GMRES seja tão, ou mais eficiente que um direto. • A estratégia de pré-condicionamento substitui o sistema linear original malcondicionado e muitas vezes indefinido por outro mais fácil de resolver e com autovalores de módulos pequenos, com parte real positiva, agrupados próximos de um 83 84 e com um número de condicionamento menor que o do sistema original sem précondicionamento. • Uma vez que o pré-condicionamento reduz o número de condicionamento e agrupa os autovalores próximos de um, as propriedades espectrais do sistema linear original são melhoradas. Assim, o escalonamento é necessário apenas se o pré-condicionamento falhar. Neste capítulo foi mostrado que para o problema de fluxo de carga são pequenas as possibilidades de occorrerem interrupções resultantes de divisões por zero durante a fatoração triangular, uma vez que os elementos da diagonal da matriz Jacobiana normalmente apresentam valores distantes de zero. Portanto, não é necessária a estratégia de Máximo Produto Transversal com Escalonamento. Apenas uma estratégia de pivotação (permutações não-simétricas) é usada de forma preventiva. • Além de reduzir o número de operações realizadas por novos elementos não-nulos, o reordenamento usa a estratégia de pré-condicionamento para melhorar a qualidade dos subespaços de Krylov gerados em cada iteração do GMRES. Entenda-se que os subespaços de boa qualidade contêm sempre uma solução muito próxima da exata e o teste de convergência pode ser satisfeito em apenas uma (ou poucas) iteração (iterações). • Foram analisados e estudados grande parte dos trabalhos realizados na área de précondicionadores aplicados a problemas de sistemas elétricos de potência, enfatizando os pré-condicionadores ILUT(τ,ρ) e ILU(k). Ambos foram avaliados e comprovou-se que o ILUT(τ,ρ) é o mais eficiente. Este pré-condicionador na sua versão mais completa, como a apresentada neste capítulo, ainda não tinha sido avaliado no problema de fluxo de carga. Observou-se também que o parâmetro ρ ajuda a melhorar o desempenho e diminui o número de operações de ponto flutuante. • Foram avaliadas e estudadas as regras de preenchimento dos pré-condicionadores mais citados, e concluiu-se que: I. A regra de preenchimento no pré-condicionador ILU(k) não permite descartar elementos das posições ocupadas da matriz Jacobiana. Além disso, sua heurística não considera nem o tamanho dos elementos nem o erro associado a cada um. O uso da regra de preenchimento dentro do algoritmo de fatoração precisa de operações adicionais para calcular os níveis de preenchimento. Portanto, essa 84 85 regra torna o pré-condicionador menos eficiente que o ILUT(τ,ρ). Os resultados nas simulações realizadas neste capítulo confirmam que o ILUT(τ,ρ) é superior. II. A heurística da regra com duplo parâmetro limitante, em alguns casos, não distingue os menores dos maiores elementos (em valor absoluto). Nesses casos a regra falha ao descartar alguns dos maiores elementos ao invés dos menores. Comprovou-se também que, nem sempre a eliminação dos menores elementos produz um erro menor que o produzido quando eliminados alguns dos maiores elementos. Neste capítulo demonstrou-se que uma regra não pode ser baseada em simples heurísticas sem considerar fundamento matemático. Portanto, a regra com duplo parâmetro limitante deve ser substituída por outra baseada no cálculo do erro associado a cada elemento descartado ou preenchido. A maior dificuldade para propor uma regra baseada no erro associado a cada elemento dos fatores triangulares é saber como estimar o erro antes da posição ser preenchida. Outra dificuldade é saber se o algoritmo de fatoração vai facilitar ou não o cálculo deste erro. Uma terceira dificuldade está em propor uma forma de calcular o erro sem que isso resulte em operações de ponto flutuante adicionais, como nas regras dos pré-condicionadores ILU(k) e ILUT(τ,ρ). III. A escolha de valores apropriados para os parâmetros do pré-condicionador ILUT(,ρ) é uma tarefa difícil. Normalmente, a escolha desses parâmetros é realizada após de um processo de tentativa e erro, cujo resultado nem sempre garante eficiência e robustez do método iterativo. • Experimentos numéricos confirmaram que os pré-condicionadores de aproximação esparsa da inversa realizam muitas operações de ponto flutuante, não sendo indicado para o problema de interesse. Os resultados confirmaram que os pré-condicionadores de fatoração incompleta podem ser mais eficazes. • A configuração mais eficiente e robusta conseguida após os estudos realizados foi o GMRES+ILUT(τ,ρ) com reordenamento MD. 85