CILAMCE
Ouro Preto/MG - Brazil
2003
XXIV IBERIAN LATIN-AMERICAN CONGRESS ON COMPUTATIONAL METHODS IN ENGINEERING
˜ DOS METODOS ´ APLICAC ¸ AO DE PONTOS INTERIORES AO PROBLEMA DE ´ ˆ FLUXO DE POTENCIA OTIMO AC UTILIZANDO COORDENADAS CARTESIANAS
Aurelio R. L. Oliveira Departamento de Matem´atica Aplicada – IMECC Universidade Estadual de Campinas, 13083-859, Campinas - SP, Brazil
[email protected] Adriano Thomaz Secundino Soares Departamento de Engenharia de Sistemas – FEEC Universidade Estadual de Campinas, 13083-970, Campinas - SP, Brazil
[email protected] [email protected] Abstract. The primal-dual interior point method are developed to the AC optimal power flow problem and the resulting matricial structure is studied. A representation of the voltage through cartesian coordinates is adopted, once that Hessian of the problem is constant and the expansion in Taylor is accurate for the second order term. The advantage of working with polar coordinates, that easily shape the voltage magnitudes, lose importance due to the efficient treatment of inequalities proportionated by the interior point methods. These methods is developed applying Newton’s methods to the optimality conditions of the problem. Before the application of the method, the number of variables of the problem is reduced through the elimination of free dual variables. This reduction does not modify the sparse structure of the problem. The linear system obtained can be reduced the dimension of twice number of buses. Moreover, such matrix is symmetric in structure. This feature can be explored reducing the computational effort per iteration. Keywords: Optimal power flow, Interior point methods, Nonlinear programming
1.
˜ INTRODUC ¸ AO
O problema de fluxo de potˆencia o´ timo e´ um problema de programac¸a˜ o n˜ao-linear de grande porte n˜ao convexo. Torna-se complicado na pr´atica pela presenc¸a de vari´aveis discretas. Dada sua importˆancia no planejamento e operac¸a˜ o de sistemas de potˆencia, os problemas de fluxo de potˆencia o´ timo tˆem sido assunto de intensa pesquisa. As ferramentas de planejamento e operac¸a˜ o de redes el´etricas precisam trabalhar com problemas com alto grau de n˜ao-linearidade no comportamento do sistema. Uma t´ecnica mais recentemente utilizada para a resoluc¸a˜ o de problemas de fluxo de potˆencia o´ timo de grande porte AC e´ a dos m´etodos de pontos interiores (Momoh, El-Hawary & Adapa, 1999, Quintana, Torres & Medina-Palomo, 2000). Problemas de programac¸a˜ o n˜ao-linear tˆem sido resolvidos eficientemente pelos m´etodos de pontos interiores para programac¸a˜ o n˜ao-linear derivados da abordagem de func¸a˜ o barreira logar´ıtmica. Embora tenha sido desenvolvida para resolver problemas de programac¸a˜ o n˜aolinear em geral, foi no campo da programac¸a˜ o linear que sua excelente eficiˆencia computacional foi primeiramente demonstrada e amplamente aceita (Adler et al, 1989, Mehrotra, 1992). 2.
´ ˆ PROBLEMA DE FLUXO DE POTENCIA OTIMO AC
O problema de fluxo de potˆencia o´ timo AC e´ um dos mais importantes na a´ rea de sistemas de potˆencia, servindo como base para diversas outras aplicac¸o˜ es. Uma dificuldade que este problema apresenta e´ instabilidade num´erica proveniente dos m´etodos de soluc¸a˜ o tradicionais. Os m´etodos de pontos interiores trouxeram a` tona uma nova linha de pesquisa na a´ rea de sistemas de potˆencia (Granville, 1994, Oliveira, Nepomuceno & Soares, 2003, Oliveira, Nepomuceno & Soares, 2001, Oliveira & Soares, 2000). Esses m´etodos s˜ao reconhecidos atualmente por sua robustez (Momoh, El-Hawary & Adapa, 1999, Quintana, Torres & Medina-Palomo, 2000). Al´em disso, o tratamento eficiente de desigualdades permite uma revis˜ao dos procedimentos geralmente adotados. 2.1 Motivac¸a˜ o A tens˜ao (complexa) da barra (E) e´ definida em coordenadas cartesianas como E = r + js, onde r e s s˜ao os componentes real e imagin´ario de (E) respectivamente. Optou-se pela utilizac¸a˜ o de coordenadas cartesianas para as tens˜oes pois desta forma, tanto as restric¸o˜ es do problema como as func¸o˜ es objetivos porventura adotadas s˜ao quadr´aticas. Conseq¨uentemente, as matrizes do problema s˜ao mais f´aceis de trabalhar e o c´alculo do termo de correc¸a˜ o do m´etodo preditor-corretor pode ser feito de forma menos custosa do ponto de vista computacional. Outra vantagem e´ que a Hessiana do problema e´ constante e a expans˜ao em Taylor e´ exata para o termo de ordem dois. Finalmente, a vantagem em se trabalhar com coordenadas polares, que modelam mais facilmente os limites de magnitude de tens˜ao, perde importˆancia devido ao tratamento de desigualdades eficiente proporcionado pelos m´etodos de pontos interiores (Momoh, El-Hawary & Adapa, 1999, Quintana, Torres & Medina-Palomo, 2000). 2.2 Formulac¸a˜ o do Problema O problema de fluxo de potˆencia o´ timo com coordenadas cartesianas pode ent˜ao ser escrito da seguinte forma, onde as equac¸o˜ es s˜ao dadas pela formulac¸a˜ o do fluxo de potˆencia apresentado anteriormente e pela representac¸a˜ o do fluxo entre as linhas de transmiss˜ao e as inequac¸o˜ es
representam a canalizac¸a˜ o de vari´aveis. minimizar 21 pt Hp + ct p sujeito a RGr + SGs + SBr − RBs = p − lp SGr − RGs − RBr − SBs = q − lq v min ≤ r2 + s2 ≤ v max fkm = gkm (rk2 + s2k − rk rm − sk sm ) + bkm (rk sm − rm sk ) pmin ≤ p ≤ pmax q min ≤ q ≤ q max f min ≤ f ≤ f max ,
(1)
onde: G representa a matriz de condutˆancia; B representa a matriz de susceptˆancia; p representa a gerac¸a˜ o de potˆencia ativa; q representa a gerac¸a˜ o de potˆencia reativa; fkm representa o fluxo de potˆencia ativa da linha k para a linha m; gkm representa a parte real da admitˆancia entre as linhas k e m; bkm representa a parte imagin´aria da admitˆancia entre as linhas k e m; H matriz diagonal representando o termo quadr´atico do custo de gerac¸a˜ o; c representa a componente linear do custo de gerac¸a˜ o; lp representa as demandas de potˆencia ativa; lq representa as demandas de potˆencia reativa; v min e v max s˜ao os limites de tens˜ao ao quadrado; pmin e pmax s˜ao os limites de gerac¸a˜ o de potˆencia ativa; q min e q max s˜ao os limites de gerac¸a˜ o de potˆencia reativa; f min e f max s˜ao os limites de fluxo de potˆencia ativa. A minimizac¸a˜ o das perdas na gerac¸a˜ o e´ utilizada como crit´erio de otimizac¸a˜ o. Essas perdas podem ser modeladas como uma func¸a˜ o quadr´atica separ´avel tanto para geradores t´ermicos representando os custos, como hidrel´etricos representando as perdas (Soares & Salmazo, 1997). Vale ressaltar que outras func¸o˜ es objetivo podem ser adotadas sem muitas alterac¸o˜ es no desenvolvimento apresentado a seguir. 2.3 Simplificac¸a˜ o da Matriz Hessiana No mesmo esp´ırito de obter Hessianas mais f´aceis de trabalhar utilizando as coordenadas cartesianas, optou-se por acrescentar a restric¸a˜ o r2 + s2 = v, onde v representa o quadrado da magnitude da tens˜ao. Assim, as equac¸o˜ es: v min ≤ r2 + s2 ≤ v max e fkm = gkm (rk2 + s2k − rk rm − sk sm ) + bkm (rk sm − rm sk ), s˜ao substitu´ıdas por: v min ≤ v ≤ v max e fkm = gkm (v 2 − rk rm − sk sm ) + bkm (rk sm − rm sk ), transformando um conjunto de restric¸o˜ es em canalizac¸a˜ o de vari´aveis e simultaneamente, simplificando outro conjunto de restric¸o˜ es. Apesar deste modelo conter um conjunto adicional de vari´aveis e restric¸o˜ es, as derivadas s˜ao mais simples e a Hessiana obtida no desenvolvimento dos m´etodos de pontos interiores mais esparsa.
3.
´ DESENVOLVIMENTO DO METODO
A aplicac¸a˜ o do m´etodo de Newton a` s condic¸o˜ es de otimalidade leva a um m´etodo de pontos interiores primal-dual espec´ıfico para este modelo. As condic¸o˜ es de otimalidade por sua vez podem ser obtidas atrav´es da func¸a˜ o lagrangiana do problema onde as restric¸o˜ es de desigualdade s˜ao representadas por func¸o˜ es de barreira logar´ıtmicas das vari´aveis de folga. Para simplificar o desenvolvimento do m´etodo, vamos considerar um problema onde cada barra pode gerar potˆencia ativa e reativa e est´a diretamente conectada a todas as outras barras. Desta forma, a representac¸a˜ o do fluxo entre as linhas pode ser escrita da seguinte forma: f = V G − RGR − SGS + RBS − SBR. Uma vez desenvolvido o m´etodo, basta considerar os elementos de f que representam linhas de transmiss˜ao existentes no sistema em estudo preservando assim a esparsidade do problema. Na pr´atica e´ muito comum resolver problemas desconsiderando os fluxos nas linhas de transmiss˜ao, pois as restric¸o˜ es tˆem utilidade somente para verificar as capacidades das linhas. Esta representac¸a˜ o do modelo pode ser facilmente alterada caso n˜ao existam muitas linhas carregadas. Com o objetivo de reduzir o n´umero de vari´aveis do problema, antes de construir a func¸a˜ o lagrangiana, vamos fazer uma mudanc¸a de vari´aveis de tal forma que todos os limites inferiores das vari´aveis canalizadas sejam anulados. O modelo adotado pode ser ent˜ao escrito da seguinte forma, onde as vari´aveis de folga para os limites superiores tamb´em s˜ao acrescentadas: minimizar φ(p) sujeito a RGr + SGs + SBr − RBs − p = −lp SGr − RGs − RBr − SBs − q = −lq V G − RGR − SGS + RBS − SBR − f = lf r2 + s2 − v = lv p + sp = pmax q + sq = q max f + sf = f max v + sv = v max (p, q, f, v, sp , sq , sf , sv ) ≥ 0,
(2)
onde φ(p) = 12 pt Hp + ct p e, por abuso de notac¸a˜ o, utilizamos os mesmos s´ımbolos para representar os vetores antes e depois das mudanc¸as de vari´aveis. 3.1 Func¸a˜ o Barreira Logar´ıtmica Para tratar as restric¸o˜ es de desigualdade das vari´aveis no m´etodo de pontos interiores, utiliza-se func¸o˜ es de barreira logar´ıtmica (Wright, 1996) que s˜ao adicionadas a` func¸a˜ o objetivo: φ(p) − µ
n X
ln(xi ),
(3)
i=1
onde x = (p, q, f, v, sp , sq , sf , sv ), n e´ a dimens˜ao do vetor x e µ > 0 representa o parˆametro de barreira.
3.2 A Func¸a˜ o Lagrangiana A func¸a˜ o lagrangiana da Eq. (2) com a barreira logar´ıtmica Eq. (3) e´ dada por (Luenberger, 1984): X L(r, s, x, l) = φ(p) − µ ln(xi ) + lt Ld (x), (4) onde lt = (yp , yq , yf , yv , wp , wq , wf , wv ) representa os multiplicadores de Lagrange e, RGr + SGs + SBr − RBs − p + lp SGr − RGs − RBr − SBs − q + lq V G − RGR − SGS + RBS − SBR − f − lf r2 + s2 − v − lv . Ld (x) = max p + s − p p max q + sq − q f + sf − f max max v + sv − v Uma soluc¸a˜ o de Eq. (2) deve satisfazer ∇(r,s,x,l) L = 0, ou seja, ∇l L = Ld (x) = 0 e: ∇p L ∇q L ∇f L ∇v L ∇r L ∇s L
= = = = = =
Hp + c − µP −1 e − yp + wp , −µQ−1 e − yq + wq , −µF −1 e − yf + wf , Gyf − µV −1 e − yv + wv , ∇r Ltp yp + ∇r Ltq yq + ∇r Ltf yf + 2Ryv , ∇s Ltp yp + ∇s Ltq yq + ∇s Ltf yf + 2Syv ,
∇sp L ∇sq L ∇sf L ∇sv L
= = = =
−µSp−1 e + wp , −µSq−1 e + wq , −µSf−1 e + wf , −µSv−1 e + wv ,
onde: Lp = RGr + SGs + SBr − RBs − p + lp , Lf = V G − RGR − SGS + RBS − SBR − f − lf , Lq = SGr − RGs − RBr − SBs − q + lq , assim: ∇r Ltp = GR + diag(Gr) + BS − diag(Bs), ∇s Ltp = GS + diag(Gs) − BR + diag(Br), ∇r Ltq = GS − diag(Gs) − BR − diag(Br), ∇s Ltq = −GR + diag(Gr) − BS − diag(Bs). ˜ = RCS onde C e´ uma matriz constante. Ent˜ao: E para ∇r Ltf e ∇s Ltf considere o termo L t (CE ) R diag (Ce S ) i−1 i−1 i−1 ˜ ˜ ∂L ∂L , = diag (Cei Si ) = (CEi )t R , ∂ri ∂si diag (Cei+1 Si+1 ) (CEi+1 )t R onde Ei = diag(ei ), i = 1, . . . , N B. Com isso podemos construir as express˜oes de ∇r Ltf e ∇s Ltf , pois Lf e´ composto de combi˜ se ignorarmos aqueles independentes de r e s. Al´em dessas nac¸o˜ es de termos similares a L, relac¸o˜ es, devemos ter x ≥ 0 o que implica (wp , wq , wf , wv ) ≥ 0. Para obter um m´etodo estritamente primal-dual resta ainda definir as vari´aveis de folga duais: zp = µP −1 e,
zq = µQ−1 e,
zf = µF −1 e,
Estas vari´aveis tamb´em s˜ao n˜ao-negativas por definic¸a˜ o.
zv = µV −1 e.
3.3 Eliminac¸a˜ o de Vari´aveis Livres y E´ poss´ıvel aplicar o m´etodo de pontos interiores primal-dual diretamente ao sistema n˜aolinear proveniente do reescalamento dos conjuntos de equac¸o˜ es referentes a` s condic¸o˜ es de complementaridade, que corresponde a` s condic¸o˜ es de otimalidade de primeira ordem do problema Eq. (2). No entanto, parece ser mais aconselh´avel eliminar as vari´aveis y do sistema antes da aplicac¸a˜ o do m´etodo. As vari´aveis y podem ser eliminadas trivialmente atrav´es das equac¸o˜ es: yp = Hp + c − zp + wp , yf = wf − zf ,
yq = wq − zq , yv = Gyf − zv + wv = Gyf + y˜v .
Obtendo o seguinte sistema:
RGr + SGs + SBr − RBs − p + lp SGr − RGs − RBr − SBs − q + lq V G − RGR − SGS + RBS − SBR − f − lf r2 + s2 − v − lv p + sp − pmax q + sq − q max f + sf − f max v + sv − v max −µe + P zp −µe + Qzq −µe + F zf −µe + V zv t t ˜ r Lt (wf − zf ) + 2R(wv − zv ) ∇r Lp (Hp + c − zp + wp ) + ∇r Lq (wq − zq ) + ∇ f ˜ s Lt (wf − zf ) + 2S(wv − zv ) ∇s Ltp (Hp + c − zp + wp ) + ∇s Ltq (wq − zq ) + ∇ f −µe + Sp wp −µe + Sq wq −µe + Sf wf −µe + Sv wv
= 0 (5)
˜ r Lt = ∇r Lt + 2RG e ∇ ˜ s Lt = ∇s Lt + 2SG. com (x, t) ≥ 0, onde ∇ f f f f 3.4 M´etodo de Pontos Interiores Primal-Dual Dada uma classe de problemas, a forma padr˜ao para desenvolver um m´etodo de pontos interiores consiste na aplicac¸a˜ o do m´etodo de Newton (Dennis & Schnabel, 1996) a` s condic¸o˜ es de otimalidade, desconsiderando as restric¸o˜ es de capacidade. A convergˆencia do m´etodo a uma soluc¸a˜ o e´ obtida partindo-se de um ponto estritamente positivo e nunca permitindo que estas vari´aveis se tornem negativas. Este controle e´ realizado atrav´es do tamanho do passo. O m´etodo resultante e´ essencialmente um m´etodo primal-dual espec´ıfico para esta classe de problemas (El-Bakry et al, 1996). O m´etodo de pontos interiores primal-dual para o problema Eq. (2) consiste portanto, na aplicac¸a˜ o do m´etodo de Newton a` Eq. (5) desconsiderando as restric¸o˜ es de n˜ao-negatividade (x, t) ≥ 0. Esta aplicac¸a˜ o resulta no seguinte m´etodo: M´etodo 1 (M´etodo de Pontos Interiores) Dados (x0 , t0 ) > 0 e (r0 , s0 ) livres. Para k = 0, 1, 2, . . ., fac¸a: ³ k´ (1) Escolha: β k ∈ (0, 1) e fac¸a µk = β k γn , onde γ k = (xk )t tk e n e´ a dimens˜ao do vetor x.
(2) Calcule as direc¸o˜ es de Newton: ∆xk e ∆tk . (3) Calcule o tamanho do passo para permanecer em um ponto interior: © ª −τ k −τ k n k o e αdk = n k o para τ k ∈ (0, 1) e αk = min 1, αpk , αdk . αpk = ∆x ∆t mini xki mini tki i
(4) Calcule o novo ponto: (x
i
k+1
,t
k+1
) = (xk , tk ) + αk (∆xk , ∆tk ).
Direc¸o˜ es de Newton As direc¸o˜ es de Newton s˜ao definidas pelo seguinte sistema linear: −∆p + ∇r Lp ∆r + ∇s Lp ∆s = r1 −∆q + ∇r Lq ∆r + ∇s Lq ∆s = r2 −∆f + G∆v + ∇r Lf ∆r + ∇s Lf ∆s = r3 −∆v + 2R∆r + 2S∆s = r 4 ∆p + ∆s = r p 5 ∆q + ∆s = r q 6 ∆f + ∆sf = r7 ∆v + ∆sv = r8 Zp ∆p + P ∆zp = r9 Z ∆q + Q∆z = r q q 10 Z ∆f + F ∆z = r f f 11 Zv ∆v + V ∆zv = r12 ˜ r Lt (∆wf − ∆zf )+ ∇r Ltp (H∆p − ∆zp + ∆wp ) + ∇r Ltq (∆wq − ∆zq ) + ∇ f 2R(∆w − ∆z ) + M ∆r + N ∆s = r13 v v t t t ˜ ∇s Lp (H∆p − ∆zp + ∆wp ) + ∇s Lq (∆wq − ∆zq ) + ∇s Lf (∆wf − ∆zf )+ 2S(∆wv − ∆zv ) − N ∆r + M ∆s = r14 W ∆s + S ∆w = r p p p p 15 W q ∆sq + Sq ∆wq = r16 Wf ∆sf + Sf ∆wf = r17 Wv ∆sv + Sv ∆wv = r18
(6)
onde: M = GYp + Yp G − (BYq + Yq B) + GYf + 2Y˜v , N = BYp − Yp B + GYq − Yq G + BYf , e os res´ıduos de r1 a r18 s˜ao dados pela aplicac¸a˜ o do ponto corrente (x, t) ao lado esquerdo do sistema de equac¸o˜ es Eq. (5) com o sinal trocado. 4.
˜ COMPUTACIONAL IMPLEMENTAC ¸ AO
Para a implementac¸a˜ o foram utilizados alguns parˆametros e vari´aveis, como a precis˜ao, representada por ², o n´umero de barras de gerac¸a˜ o N G, n´umero de barras com limites de gerac¸a˜ o de potˆencia reativa N H e tamb´em representamos o n´umero total de barras por N B. 4.1 Atualizac¸a˜ o das Vari´aveis As novas vari´aveis primais e duais s˜ao calculadas da seguinte maneira: xk+1 = xk + αk ∆x
e
tk+1 = tk + αk ∆t,
onde o escalar αk ∈ (0, 1] e´ o parˆametro de comprimento do passo.
(7)
4.2 C´alculo do Comprimento do Passo O comprimento m´aximo do passo αk e´ determinado por: ½ k¾ ½ k¾ ª © kmax −xi −ti kmax kmax kmax k αp = min , α = min , α = min τ α ; τ α ; 1, 0 . d p d ∆xki ∆tki ∆xki
