Agnieszka Ligocka, Wojciech Bandurski ALGORYTMY PRZYSPIESZAJACE OBLICZENIA MACIERZY TRANZYCJNEJ DUŻYCH UKŁADÓW W ZASTOSOWANIU DO MODELOWANIU POŁĄCZEŃPolitechnika Poznańska, Instytut Elektroniki i Telekomunikacji

2003

Poznañskie Warsztaty Telekomunikacyjne

Poznañ 11-12 grudnia 2003

Wojciech Bandurski Agnieszka Ligocka Politechnika PoznaĔska

Instytut Elektroniki i Telekomunikacji ul. Piotrowo 3a, 60-965 PoznaĔ Wojciech.Bandurski@put.poznan.pl Agnieszka.ligocka@et.put.poznan.pl

ALGORYTMY PRZYSPIESZAJACE OBLICZENIA MACIERZY TRANZYCJNEJ DUĩYCH UKàADÓW W ZASTOSOWANIU DO MODELOWANIA POàĄCZEē

Streszczenie: Artykuá bazuje na wczeĞniejszych pracach autorów dotyczących modelowania poáączeĔ. Celem artykuáu jest wyraĪenie macierzy admitancji w dziedzinie czasu za pomocą zaleĪnoĞci analitycznej. Macierz tranzycjną wystĊpującą w tym wyraĪeniu obliczamy za pomocą algorytmu Schura–Frecheta. W artykule przedstawiamy wyniki obliczeĔ funkcji

wykáadniczej dla wersji rekursywnej i

semirekur-sywnej algorytmu Schura–Frecheta. Uzyskane wyniki zilustrowano przykáadem.

1.

WPROWADZENIE

W dziedzinie nauk elektronicznych i elektrycznych metody symulacji i modelowania mają ogromne znaczenie. PostĊp technologiczny powoduje jednak, Īe symulowane obwody są coraz wiĊksze i coraz bardziej záoĪone. W czasie ich analizy powstają równania liniowe, których macierze zawierają nawet dziesiątki tysiĊcy elementów. W symulacji ukáadów elektronicznych coraz wiĊkszą rolĊ odgrywa zagadnienie modelowania poáączeĔ w szybkich ukáadach cyfrowych, które stanowią podstawowe Ĩródáo opóĨnieĔ i znieksztaáceĔ sygnaáu. W związku z tym coraz wiĊksze znaczenie ma analiza i modelowanie linii transmisyjnych. Modelowanie i symulacja jest praktycznie jedynym sprawdzianem poprawnoĞci zaprojektowanego ukáadu przed jego fizyczną realizacją. Ze wzglĊdu jednak na duĪą záoĪonoĞü takich modeli, do ich analizy nie nadają siĊ tradycyjne metody stosowane do rozwiązywania ukáadów elektronicznych. Aktualnie najczĊĞciej [2] rozwaĪa siĊ modele poáączeĔ zbudowane z elementów skupionych RLC.

Analiza takiego obwodu elektrycznego opiera siĊ na rozwiązaniu ukáadu liniowych równaĔ róĪniczkowych postaci (1): x L i Bu Gx x C N N T dt d  (1) gdzie:

C, G, B, L – staáe macierze rzeczywiste

u – N–wymiarowy wektor reprezentujący m napiĊü na wrotach

i – N–wymiarowy wektor reprezentujący n prądów na wrotach

x – n–wymiarowy wektor zmiennych stanu

Rozwiązanie ukáadu (1) sprowadza siĊ do rozwiązania pierwszego równania przy pomocy metod numerycznych lub analitycznych, czyli przez poszukiwanie odpowiedzi w dziedzinie czasu. MoĪna takĪe rozwiązywaü problem prze poszukiwanie biegunów i zer transmitancji – sprowadzając to zagadnienie do poszukiwania wartoĞci wáasnych dwóch macierzy rzeczywistych, dla którego opracowano wiele metod numerycznych.

Dla duĪych ukáadów równaĔ (1) powstaje problem

záoĪonoĞci obliczeĔ związanych z wyznaczeniem

rozwiązania. Dlatego prace poĞwiĊcone temu zagadnieniu koncentrują siĊ gáównie na problemie redukcji rzĊdu ukáadu [6,7,8,9]. Redukcja rzĊdu ukáadu pozwala na stworzenie zastĊpczego modelu ukáadu o mniejszej liczbie

elementów. DziĊki temu ukáady takie moĪna

wykorzystując model parametrów zmniejszonej liczbie

parametrów symulowaü za pomocą narzĊdzi takich jak SPICE. Takie podejĞcie wykorzystują miĊdzy innymi

algorytmy Arnoldiego i Lanczosa [2,3,4], bazujące na

metodzie podprzestrzeni Kryáowa. Wadą metod opartych na podprzestrzeniach Kryáowa jest jednak brak moĪliwoĞci oszacowania rzĊdu ukáadu zredukowanego a priori. Dlatego istotne są dalsze badania nad innymi bardziej precyzyjnymi metodami modelowania.

W artykule przedstawiono metodĊ rozwiązania duĪego ukáadu równaĔ (1). PodejĞcie to wykorzystuje analityczne rozwiązanie ukáadu równaĔ (1), rekurencyjne algorytmy splotowe, a do wyznaczenia macierzy

tranzycyjnej ostatnio zaproponowany algorytm Schura-Frecheta [1] w wersji rekursywnej oraz semirekursywnej.

2.

TEORIA

Zagadnienie analitycznego rozwiązania równania (1) sprowadza siĊ do wyznaczenia równania postaci:

    » » ¼ º « « ¬ ª 

_³

    t t N t t t T N e t e d 0 0 0 0 _{( ) ( )} 0 W W W Wu x L i A A , (2) gdzie: G C A0 1  B C

W 1 x

 

t0 – warunki początkowe.

W analitycznym rozwiązaniu ukáadu (1) pojawia siĊ

macierz tranzycyjna exp(-A0t). Funkcja wykáadnicza

exp(-A0t), gdzie A0 jest macierzą kwadratową rzĊdu n,

definiowana jest nastĊpującym szeregiem:









, ! 1 ! 2 1 1 0 2 0 0 0t _{   t}n_ n t t eA A A A (3)

Funkcja ta jest jak wynika z (3) macierzą kwadratową stopnia n, której elementy są funkcjami czasu. Jej obliczenie stanowi gáówną trudnoĞü rozwiązania

analitycznego. Dla ukáadów duĪych rzĊdów macierz A0

jest bardzo duĪa, co uniemoĪliwia wykorzystanie standardowych metod numerycznych obliczania macierzy tranzycyjnej. Dla duĪych macierzy wiĊkszoĞü algorytmów nie daje oczekiwanych rezultatów. Jest to spowodowane zarówno duĪymi nakáadami obliczeniowymi potrzebnymi do uzyskania wyniku, jak i záym uwarunkowaniem numerycznym macierzy. Macierze okreĞlające parametry ukáadu są zazwyczaj macierzami rzadkimi. Do wyznaczania funkcji wykáadniczej duĪych macierzy moĪna zastosowaü dwie techniki. Jedną z nich jest aproksymacja Padego. MoĪe ona byü jednak wykorzystywana tylko wtedy, gdy norma ||A0t||2 jest maáa. Gdy jednak wartoĞü

normy jest duĪa metoda aproksymacji Padego okazuje siĊ

maáo dokáadna. Drugie podejĞcie bazuje na metodzie „potĊgowania i pierwiastkowania” (ang. scaling and squaring) wykorzystującą nastĊpującą zaleĪnoĞü:

k k e e 2 2 _¸ ¹ · ¨ © § X X _{. (4)}

Dobór wartoĞci parametru k uzaleĪniony jest od dwóch czynników. Jego wartoĞü musi byü na tyle duĪa by

exp(X/2k) moĪna byáo áatwo aproksymowaü przez

rozwiniĊcie wielomianowe lub wymierne np. przy pomocy aproksymacji Padego. Jednak gdy k jest zbyt duĪe to

exp(X/2k)~I. WáaĞciwy dobór wartoĞci k powinien

zapewniü dobre uwarunkowanie numeryczne obliczeĔ oraz osiągnąü wymaganą dokáadnoĞü. Nie jest to jednak proste. Dodatkowo utrudnieniem jest fakt, Īe w róĪnych fragmentach macierzy wartoĞci mogą znacznie siĊ od siebie róĪniü, wiĊc korzystne byáoby skalowanie róĪnych fragmentów macierzy z róĪnymi wartoĞciami wspóáczynnika k.

Korzystając z (2) po przyjĊciu zerowych warunków początkowych oraz podstawieniu t0=0, otrzymamy macierz

admitancji dla N-tego wyjĞcia:

W LTe At t

y₍₎  0 _(5a)

oraz w dziedzinie zmiennej s:

>

@

L

>

A

@

W Y( ) () 1 0 1  s t y s L T (5b)

3.

ALGORYTM SCHURA FRECHETA

Przy zaáoĪeniu, Īe norma macierzy ||A0t||2 oraz

rozmiary macierzy są duĪe oraz rozmiary macierzy są duĪe obliczenie funkcji wykáadniczej duĪej macierzy exp(-A0ǻt)

jest trudne. Do jej obliczenia zastosowano w przedstawianym podejĞciu metodĊ Schura-Frecheta. Jest to ostatnio opracowany efektywny algorytm obliczania duĪych macierzy wykáadniczych w oparciu o metodĊ „scaling and squaring” (4). Metoda wykorzystuje

dekompozycjĊ macierzy A0ǻt na macierz trójkątną przy

pomocy dekompozycji Schura oraz oblicza pochodną kierunkową Frecheta z macierzy R powstaáej w wyniku dekompozycji.

Kolejne kroki algorytmu obejmują:

1.

DekompozycjĊ Schura macierzy A= A0ǻt

R = QTAQ .

2.

DekompozycjĊ macierzy R » ¼ º « ¬ ª  » ¼ º « ¬ ª  0 0 R 0 R 0 0 R Z D R 12 22 11

macierz Z jest macierzą trójkątną górną i nilpotentną.

3.

RozwiniĊcie exp(R)

exp(R) = exp(D) + Lexp(Z,R) gdzie

Lexp (Z,R) – pochodna kierunkowa Frecheta z exp(R)

w kierunku macierzy Z .   _{e ds} e  e _ds e L t s s t s s

³

³

» ¼ º « ¬ ª   0 12 1 0 1 exp 22 11 ) , ( 0 0 R 0 Z R Z R R R R _.

4.

Obliczenie Lexp(Z,R) poprzez zastosowanie metody

“scaling and squaring” (2) oraz metody áaĔcuchowej dla pochodnej Frecheta [1]

k R₁₂/2 0 12 R   _{e ds} e k s k t s 0 /2 12 0 2 / 1 1 12 11 22 R R R R

_³

 k k e e 1 /2 12 1 12 2 / 2 12 11 22 R R R R R  (7) … i k i k e e i i i 11 1 _ 22/21 12 12 2 / 1 12 R R R R R … » » ¼ º « « ¬ ª _ 0 0 R R 0 R Z R R k k e e L k k /2 12 12 2 / exp 22 11 ) , ( . W celu obliczenia 1 12 R , w artykule wykorzystano

metodĊ trapezów. W kolejnych krokach obliczano wartoĞci

i

12

R aĪ do uzyskania Lexp(Z,R)MoĪna wyodrĊbniü dwie

zasadnicze metody pracy algorytmu Schura Frecheta:

kaĪdym kroku obliczania pochodnej Frecheta wykorzystuje siĊ kaĪdorazowo rekursywne obliczanie macierzy wykáadniczej (mniejszego rzĊdu niĪ macierz wyjĞciowa), w kaĪdym kroku obliczania pochodnej Frecheta konieczne jest wiĊc ponowne obliczenie wartoĞci macierzy wykáadniczej obliczonej w poprzednim kroku, jednak odpowiednio skalowanej, w zaleĪnoĞci od wartoĞci normy macierzy obliczanej w tym kroku. Wersja semirekursywna, zamiast obliczania wprost wartoĞci macierzy wykáadniczej, wykorzystuje poprzednio obliczone wyniki skalując je przy pomocy pierwiastkowania numerycznego macierzy. MoĪe to prowadziü do powstania báĊdów numerycznych, dlatego wersja semirekursywna zalecana jest jedynie gdy norma macierzy znajdującej siĊ w wykáadniku nie jest duĪa[1]. Jej zaletą jest jednak krótszy czas obliczeĔ ze wzglĊdu na mniejszą iloĞü mnoĪeĔ.

Zasada pracy algorytmu jest nastĊpująca. Algorytm

rozpoczyna siĊ wyznaczeniem wartoĞci

 

a

ii R exp oraz

 

b ii R

exp na diagonalii macierzy R (Rys.1). Po

dekompozycji Schura wartoĞci na diagonalii są bliskie wartoĞciom wáasnym macierzy R. Pozwala to na obliczenie wartoĞci na dialgonalii wprost jako funkcji wykáadniczej wartoĞci skalarnej.

R12 R22 a 11 Ra R11b R22 b R12b R12a 0

Rys.1 Schemat algorytmu obliczania pochodnej Frecheta NastĊpnie stosując pochodną Frecheta algorytm oblicza exp R

 

a_ii,b . W ten sposób uĪywając do obliczeĔ coraz wiĊkszych fragmentów macierzy R algorytm oblicza caáe wyraĪenie exp(R) poprzez podwajanie rozmiaru macierzy w kaĪdym kroku. JeĪeli rozmiar macierzy R nie jest dokáadnie równy potĊdze 2 wtedy macierz musi zostaü uzupeániona przez dodanie zerowych wierszy i kolumn.

4.

PRZYKàAD

Dziaáanie algorytmu zaprezentowano na przykáadzie analizy stratnej linii transmisyjnej trójprzewodowej sprzĊĪonej (Rys.2). Dla celów analizy zostaáa ona dyskretyzowana na 40 segmentów, kaĪdy záoĪony z elementów RLC.

Rys.2 Model segmentu stratnej linii transmisyjnej prezentowany w przykáadzie

Parametry linii wynosiáy odpowiednio:

>

@

>

/

@

, 4.98 0.76 0.15 0.76 4.98 0.76 0.15 0.76 4.98 , / 1.08 0.197 -0.006 -0.197 -1.124 0.197 -0.006 -0.197 -1.08 cm nH cm pF » » » ¼ º « « « ¬ ª » » » ¼ º « « « ¬ ª L C R=diag(3.5)[:/cm], dáugoĞü linii d=5 cm.

Równania obwodu dla N sekcji, które powstaáy po dyskretyzacji macierzowego równania telegraficznego [5] moĪna przedstawiü w postaci:

N m dt d n n n n n n G v i i 1 1,..., v C 1 ,..., 1 1     m N dt d n n n n n n R i v v i L

W przedstawionym przykáadzie po dyskretyzacji przy pomocy algorytmu Schura-Frecheta obliczono macierz

exp(-A0't). ObliczeĔ dokonano przy pomocy wersji

rekursywnej, semirekursywnej algorytmu oraz dla porównania funkcji expm, która jest standardową funkcją pakietu MATLAB. Porównanie otrzymanych wyników oraz czasu obliczeĔ przedstawiono w tabeli 1. W celu obliczenia pierwiastka macierzy w wersji semirekursywnej korzystano z funkcji pakietu MATLAB – sqrtm. Dla

maáych wartoĞci parametru k dodatkowe przyspieszenie

pracy algorytmu uzyskiwano przez zastąpienie kolejnych pierwiastkowaĔ potĊgą uáamkową, jednak obliczanie potĊgi uáamkowej macierzy obarczone jest dodatkowymi báĊdami, dlatego dla wiĊkszych macierzy, metoda ta okazaáa siĊ nieefektywna.

Tabela 1 Porównanie wyników dziaáania algorytmów obliczania macierzy wykáadniczej

Maksymalna wzglĊdna róĪnica

algorytm

Czas obliczeĔ

[s] Do expm _rekursywnej Do wersji

Expm 0.313 - 7.0074 e-007 SF w. rekursywna 13.93 7.0074 e-007 -SF w. semirekur-sywna (sqrtm) 4.96 0.7032 0.7032 SF w. semi-rekursywna (potĊga uáamkowa) 1.40 2.8958 e+094 2.8050 e+094

Badania przeprowadzone przez Kenney i Laub dowodzą, Īe dla macierzy Ĩle uwarunkowanych metoda expm daje wyniki mniej dokáadne niĪ obliczenia przy pomocy algorytmu Schura Frecheta [1]. Na podstawionego powyĪszego przykáadu moĪna zauwaĪyü, Īe w tym przypadku wyniki te są zbliĪone, dlatego korzystniejsze jest stosowanie szybszych algorytmów. W przypadku

omawianej macierzy jest nią metoda expm.

Przeprowadzone badania dowodzą jednak, Īe metoda ta nie zawsze daje sobie radĊ z obliczeniami (dla macierzy Ĩle uwarunkowanych), zwiĊksza siĊ teĪ znacznie czas obliczeĔ. (Dla macierzy o tym samym rozmiarze co w

przykáadzie - 243x243 nawet do 450s) podczas gdy czas

obliczeĔ metodą Schura Frecheta dla tego samego przykáadu wyniósá 69,7s.

Kolejnym etapem badaĔ byáo porównanie wyników

dokáadnych z wynikami obliczonymi przy pomocy

algorytmu Schura Frecheta w wersji semirekursywnej. W tym celu wykorzystano wyraĪenie:





W

L

T A t m m

a

t

e

o

y

(

)

 ' (tm = m 't, m=0,...,M)

i otrzymano macierz admitancji w (M+1) punktach czasowych.

5.

OTRZYMANE WYNIKI

Wykorzystanie przeksztaácenia FFT wyników

algorytmu Schura Frecheta Ya(jZ) pozwoliáo na

porównanie otrzymanego wyniku z metodą dokáadną (równanie 5b). Porównanie to jest przedstawione na wykresach (Rys 4, 5, 7, 8). MoĪna zauwaĪyü, Īe metoda Schura Frecheta daje wyniki pokrywające siĊ z metodą dokáadną. Lepsze rezultaty uzyskuje siĊ dla czĊĞci urojonej wyniku. Wykresy przygotowano w oparciu o metodĊ

semirekursywną algorytmu Schura Frecheta obliczania

macierzy wykáadniczej ze wzglĊdu na krótszy niĪ w przypadku metody rekursywnej czas obliczeĔ, przy zachowanej wymaganej dokáadnoĞci.

Na wykresach przedstawiono nastĊpujące zaleĪnoĞci:

Wyniki bezpoĞredniego obliczenia macierzy admitancji

(4b):





>

ij

@

ij T j i s L y t Y, ( ) , () _, dla

i

,

j

0

,

1

...

5

Wyniki analizy czasowej macierzy admitancji otrzymane przy pomocy algorytmu Schura Frecheta:





>

@

ij m W t A T m j i SF

t

L

e

y

,

(

)

, ' ˜  dla

i

,

j

0

,

1

...

5

Przedstawione poniĪej wykresy obrazują odpowiedĨ czĊstotliwoĞciową ukáadu na napiĊciowe wymuszenie impulsowe. Przedstawiono wyniki dla elementu (0,0) i (4,5) macierzy admitancji.

Rys. 3 Element (0,0) – odpowiedĨ czĊstotliwoĞciowa (amplituda)

Rys. 4 Element (0,0) – odpowiedĨ czĊstotliwoĞciowa (faza)

Rys 5 Element (4,5) – odpowiedĨ czĊstotliwoĞciowa (amplituda)

Rys 6 Element (4,5) – odpowiedĨ czĊstotliwoĞciowa (faza)

6.

PODSUMOWANIE

W artykule zaprezentowano algorytm Schura

Frecheta obliczania macierzy wykáadniczej oraz

moĪliwoĞci jego zastosowania w obliczaniu macierzy admitancji w dziedzinie czasu. Przedstawiono równieĪ prosty przykáad zastosowania powyĪszej metody obliczania oraz porównanie róĪnych metod algorytmu Schura Frecheta. Planowane są dalsze prace nad oceną przydatnoĞci prezentowanego algorytmu w analizie duĪych obwodów RLC.

SPIS LITERATURY

[1] C. S. Kenny, A. J. Laub, A Schur-Frechet algorithm for computing the logarithm and exponential of a matrix, SIAM J. Matrix Anal. Appl., vol. 19, 1998. [2] A. Odabasioglu, M. Celik, L. T. Pileggi, PRIMA:

Passive reduced-order interconnect macromodeling algorithm, IEEE Tran CAD, vol. 17, August, No. 8, 1998, 645-654.

[3] R. W. Freud, Krylov-subspace methods for reduced-order modeling i circuit simulation, J. of Comput. Appl. Math., (123): 395-421,2000.

[4] C. S. Kenny, A. J. Laub, Condition estimates for matrix functions, SIAM J. Matrix Anal. Appl., vol. 10, 191-209, 1989.

[5] L. Knockaert, D. De Zutter, Laguere-SVD reduced-order modeling, IEEE Tran. MTT, vol.48, No. 9, Sept., 2000, 1469- 1475.

[6] A. Ligocka, W. Bandurski, Schur-Frechet Algorithm For Simulation Large RLC Circuit, Proc of 7th IEEE Workshop on SPI May 2003, 155-156

[7] B. Gustavsen, A. Semylyen, Enforcing Passivity for Admittance Matrices Approximated by Rational Functions. IEEE Trans. Power Systems Vol 16 No 1 February 2001, 99-104

[8] B. Gustavsen, A. Semylyen, Rational approximation of Frequency Domain Responses by Vector Fitting. IEEE Trans. Power Delivery Vol 14 No 3 July 1999, 1052-1061.

[9] R. Neumayer, A. Stelzer, W. Enrsens, R. Weigel „Package Modeling Using Subspace Identification Algorithms” Proc of 7th IEEE Workshop on SPI May 2003, 29-32.