TELKOM NIKA , Vol. 11, No. 5, May 2013, pp. 2583 ~   2593   ISSN: 2302-4 046           2583      Re cei v ed  No vem ber 2 9 , 2012; Re vi sed  March 13, 20 13; Accepted  March 23, 20 13   The Research of Digital Algorithm Based on Frequency- Dependent Transmission Lines      Yongqing Liu 1,2 , Baina He* 3 , Yun w ei  Zhao 4 , Heng xu Ha 3 , Xinhui Zhang 3   1 Colle ge of el e c tric engi neer in g and a u tomati on, F u zhou U n i v ersit y , F u zh ou , China, 35 01 0 8   2 F u jian e l ectric al po w e r comp an y, F u zhou, C h in a, 350 003   3 Colle ge of Ele c trical an d Elec tronics Eng i ne erin g,  Shan don g Univ ersit y   of T e chnolog y, Z i Bo, Chin a,  255 04 9. T e l:  1368 53 313 24   4 Departme n t of Electric Engi n eeri ng, Sha ndo ng Indus tr y Po l y tec hnic C o ll eg e, Z i Bo, China,  2564 14.    *Corres p o ndi n g  author, e-ma i l : hbn7 70 425 @ 163.com       A b st r a ct   T he alg o rith m for obtaini ng  the discrete  re spo n se of p r opa gatio n fun c tion for frequ ency   dep en dent p a r a meter li ne  is prese n ted.  Co nsid er a  mi ni mum s a mpl i n g  p e rio d  T s m, tha t  is, the hig h e s t   freque ncy fH= 1 /(2T sm) i n  th e sign al is  tak en int o  acco un t. T he imp e d a n ce  z ( )   an d the a d m ittance  y ( )   are obta i ne d in  the frequency  rang e of  [0,fH]  by empl oyin g the Cars on s  fo rmu l a. T he pro pag atio n functi on   a t  e a c h  freq u e n cy p o i n t  i s  sub s e q u e n t l y  o b t a i ne d ,  th e  impu l s e  re sp on se   i n  d i scre t e  time  d o ma i n  i s  then  obtai ne d usin g  Poisio n Su m F o rmu l a. In or der to av oi d the lon g  len g th  of imp u ls e res pons e un der the   hig her s a mpl i n g  freq ue ncy, the  pol es  and   z e ros  of   z   tr ansfor m   of di screte pr opa g a tion  functi on  a r e   eval uate d  by t he Pro n y method. S ubse q u e n tly, the co effi cients of th di screte infi nite  i m p u ls e res pon se   of prop ag ation  function   are  obtai ne d. Usin g these  coeffic i ents  the  w a ve  transfer so urces can  be  eas i l y   compute d  by di screte convo l ut ion o perati on. T he  simul a tion tests  show   that the results usi ng the pro pos e d   meth od is acc u rate, the error i s  not mor e  tha n   1% in co ntra st of the  results gener ated by  EMT P   Ke y w ords : Propa gati on F u n c tion Co efficie n t; F r equency  Dep end ent Mo del; Prop ag atin g T r aveli ng W a ves         Copy right  ©  2013 Un ive r sita s Ah mad  Dah l an . All rig h t s r ese rved .       1. Introduc tion  In the ea rly 80’s  of the 2 0 th ce ntury,  J.  Marti  pro p o se d a tran si ent mod e l fo r the  freque ncy - de pend ent pa rameter lin es,  in whi c h t he dist ribute d  paramete r s vary with  the   freque ncy [1 -5]. In Marti’s mod e l, the surge impe dan ce  Z c ( )  is synthe sized with an  RC  netwo rk, in  which the  con s tants R’ s and   C’s a r e dete r mined by the  pole s  and  ze ros of   Z c ( ) , the  equivalent m odel is  sho w n  in Figure 1.           Figure 1. J M a rti’s Equival ent Model       In the equiva lent ci rcuit  of Figure 1, th wave t r an sfer source s a r e repe ctively the   traveling waves propa gate  from the othe r terminal.     d f t a t f t a t e d b t a t b t a t e M L M L N N L N L M ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) (   (1)    Evaluation Warning : The document was created with Spire.PDF for Python.
                               ISSN: 23 02-4 046   TELKOM NIKA  Vol. 11, No . 5, May 2013 : 2583 – 259 3   2584 Whe r e:  L a  is propag ation fun c tion,   N b  is backward traveling  wave   M f   is  forward traveling wave.   The propa ga tion function  is obtain ed b y  taking inve rse  Lapl ce transfo rm of  A L ( ) w h er   ) ) ( ) ( exp( ) ( l y z A L   (2)     The z and  y  are re sp ectiv e ly  the seri es  imped an ce and sh unt  a d m ittance   pe r length  of the transmissi on line,  which  are frequ en cy-d ep enda nt [6-8]. In orde r to get the inverse   Lapla c e t r an sform  of the  prop agatio n f unctio n , Ma rt i pro p o s ed th at the fun c tio n  in fre que ncy  domain  can b e  approximat ed to be the synthesi s  by a sum of first o r der terms:     ) ( ) ( 1 min P m m m T j L p s K e A   (3)     Subse que ntly, the propag ation functio n   in time domain ca n be  obtained by  the  following formula:    P m T t p m L m e K t a 1 ) ( min ) (   (4)     Obviou sly, the di screte v a lue  of tran sfer  so urse  i n  eq uivalent  circuit  (Fig u r e 1 )   requi re the discrete co m putation  of convolution int egral  with  forward an d ba ckwa rd traveling   wave s by em ploying form u l a (1). Sup p o s e that if it  directly get the  discrete im pu lse respon se  o f   the prop agati on functio n   a L ( k ), the tra n s fer  sou r ces  can b e  di re ct ly obtained b y  the discrete   convol ution al gorithm, that is,    n m M L N n m N L M m k f m a k e m k b m a k e 1 1 ) ( ) ( ) ( ) ( ) ( ) (  (5)     This could b e  simplify the comp utatio n of  transie n t  analysis. In this pape r,  the   Poisson’ s Su mmation fo rmula i s  e m p l oyed to  obt ain the  digit a l coefficie n ts of  propag a t ion   function,  a L ( k ) .  Ho wever,  the discrete prop agatio n functio n  ca n be re garded  as an Infini te  Impulse  Re spon se (IIR) d i gital filt er [9-10], that is the length of  a L ( k ) is infinit e . The Prony ’s   method is e m ployed to cal c ulate the ze ro s and p o le s o f  the z transfo rm of  a L (k ) Therefore the  traveling wa ves pro pag ating to  anothe r terminal (tra nsfer  sou r ce s) can   be cal c ul ated  in discrete t i me domai with le ss am ount of co m putati on. Co mpared with  the  results by the  ne w te ch niq ue  with th at b y  PSCAD, th e results  sh o w  that  the  propo sed  alg o ri thm  is efficien cy a nd accu rate e v en at a con s idera b le lareg e  step s.       2. Basic Prin ciple  2.1. Frequen c y  Depende nt Parame ter  Line  For a  singl e  pha se tra n smissi on lin e MN which le ngth is  l , associate d  with  the  freque ncy d e pend ent dist ri buted pa ram e ters, the lin e  equation s  in  freque ncy d o m ain at location   x is shown as follows:   Evaluation Warning : The document was created with Spire.PDF for Python.
TELKOM NIKA   ISSN:  2302-4 046     The Re se arch of Digital Algorithm  Base on Fre que n c y-Depe nde n t  ... (Yongqin g  Liu)  2585 ) , ( ) ( ) , ( ) , ( ) ( ) , ( j x U j y dx j x dI j x I j z dx j x dU  (6)     Whe r e,  and  y are  re spe c t i vely se ries i m peda nce a n d  shunt  admit tance  pe r l e n g th of  line. Th e val ue of  freq uen cy-de pen dent  impe dan ce   and  admittan c at ea ch  freque ncy  ca n  be  obtaine d by Carso n ’s fo rmula [11-13].  One  can  get  the pro pag a t ion relatio n s of forwa r a nd ba ckward  traveling  wa ves in  freque ncy do main:    ) ( ) ( ) ( ) ( ) ( ) ( j B j A j B j F j A j F N L M M L N  (7)     Whe r e  F=U+ZcI i s  fo rwa r d travelin wave  an d B = U-ZcI i s  ba ckward  travelin wave.  Zc is  surge i m peda nce, A L  is prop agati on fucntio n   ) ( ) ( ) ( y z z c  (8)     ) ) ( ) ( exp( ) ( l y z A L  (9)     2.2. Poisson’s Summation Formula  Acco rdi ng to the theories of signal an d syst em s, the discrete Fouri e r tran sf orm of a  sign al x(n )  i s  the p e rio d ic extensi on of  the Fo urie r transfo rm of its  corre s po n d ing contin uo us  sign al x(t). That is, if x(n) is the discre te si gnal by  sampli ng the  signal x(t)  with the sampl i ng   freque ncy  N, the Fou r ier transfo rm of x(t) is   X ( ) , then the Fou r ie r transfo rm of  discrete  sig n a x( n )  is   s h ow n a s   k kN X X ) 2 ( ) ( ~  (10 )     Usi ng Poi s son’s  sum fo rmula, one  can get  the  discrete  sig n a ls from th e  Fourie r   transfo rm of continuo us fu cntion:     kk N jk e N k X N kN X ) ( ˆ 2 1 ) 2 (  (11 )     Whe r e:     N N a k j d e X N k X / ) ( ) ( ˆ  (12 )     Then the  z transfo rm of di screte  signal  c an b e  obtain ed as the foll owin g formul a.    k k z N k X N z X ) ( ˆ 2 1 ) (  (13 )     Therefore the  discrete  sign al is obtain e d :   Evaluation Warning : The document was created with Spire.PDF for Python.
                               ISSN: 23 02-4 046   TELKOM NIKA  Vol. 11, No . 5, May 2013 : 2583 – 259 3   2586 ) ( ˆ 2 1 ) ( N k X N n x  (14 )     Certai nly usi ng Poisson’ s Sum Formul a, one  ca n e a sily get the discrete p r o p agation  function  a nd discrete su rg impeda nce.    2.3. The Disc rete Prop aga t ion Func tio n  and Surge  Impedance   The discrete  Fouri e r tra n sf orm of propa gation fun c tio n  can b e  writt en as:     k L L kN A A ) 2 ( ) ( ~  (15 )     By employing  Poisso n’s  su m formula,  propag ation fun c tion will b e  written a s :     k k L k a k j L L z k a e k a z A ) ( ) ( ) ( /  (16 )   Similar proce dure fo r su rg e impeda nce:    k k c k a k j c c z k Z e k Z z Z ) ( ) ( ) ( /  (17 )     In ord e r to  si mplify the co mputation, th e minim u m traveling time   N mi n  is se pa rated fro m   A L   ) ( 0 min z A Z A L N L  (18 )     whe r Nmin  is the minim u m delaye d  time of pro p a gation traveli ng wave s. It can b e   determi ned b y  the light velocity:    ) / ( floor min s cT L N  (19 )     By employin g Poi s son  sum formula,  re cu rsive  d i screte  sequ ence of  pro pagatio function a nd  surge imp eda nce  will be ob tained.         ) 2 exp( ) ( 2 1 ) ( ) 2 exp( ) ( 2 1 ) ( 0 0 0 0 f N n f k j f k Z N n z f N n f k j f k A N n a N k s c s c N k s L s L  (20 )     Ho wever, the  propa gation  function is  an Infini te Impulse Lo w Pass filter, that is, its  discrete  coefficient i s  infinit e  long. In  ord e r fo r sim p lified computati on, the z tran sform  of AL a n d   Zc re quire to be expre s sed  as re cu rsive  form:    n N m M N L z z z z Z z A ... 1 ... ) ( 1 1 1 1 0 min  (21 )     So doe s the surge imp eda nce. It need t o  deter mi ne the co efficient s of a and b f r om the   discrete  pro p agation fun c t i on and  su rge imped an ce. The Pron y’s method i s  empl oyed  to   es timate these c oeffic i ents .     2.4. The Rec u rsiv e Coeffi cients by  Pr on y s Metho d   The c o effic i ents   k  can be  estimated by  solving h o mo gene ou s differen c e eq uatio n.  Evaluation Warning : The document was created with Spire.PDF for Python.
TELKOM NIKA   ISSN:  2302-4 046     The Re se arch of Digital Algorithm  Base on Fre que n c y-Depe nde n t  ... (Yongqin g  Liu)  2587 N L L L L L L L L L L L L N p n a n a p n a N n a n a n a N n a n a n a p n a n a n a ... ) ( ... ) 2 P ( ) 1 ( ... ... ... ... ) 1 ( ... ) 3 ( ) 2 ( ) ( ... ) 2 ( ) 1 ( ) ( ... ) 1 ( ) ( 2 1 0 0 0 0 0 0 0 0 0 0 0 0  (22 )     By means of Lea st Square  Method, coef ficient    can  be estimate d as:     L0 T 1 T ) ( A B B B  (23 )       The z tra n sfo r m of AL 0  can  be synthe sized by a sum  of first orde r term s.    k k k n N m M L z z b z z z z z A 1 1 1 1 1 0 0 1 ... 1 ... ) (  (24 )     Whe r z k  is the pole s  dete r mine d by the following formula.    ]) ,..., , 1 [ roots( 1 N k z  (25 )       B k  can be cal c ulate d  by usi ng the followi ng formul a.    N n N n n n N n n n N n n L L L b b b z z z z z z z z z p n a n a n a ... ... ... ... ... ... ... ... ) ( ... ) 1 ( ) ( 2 1 P P 2 P 1 1 1 2 1 1 2 1 0 0 0  (26 )     By employing  Least Squa re Method, co efficient b ca n be solve d  a s    L0 1 A H H Z Z Z b  (27 )     At finally,  k z k b have been o b tained, then  we can g e t pro pagatio n function  ) ( z A L   N 1 1 k min 1 ) ( k k N L z -z b Z z A N N M M N z z z z Z 1 1 1 1 0 min 1  (28 )     The p r o c edu re to obtainin g  the re cursive  coe ffici ents  of prop agatio n functio n  an d su rge   impeda nce is sho w n in Fig u re 2.     2.5. Propaga tion Func tio n  Verificatio On the  ba sis of the frequ ency  depe nd ent mod e l, propag ation fu nction  will  be  verified.   whe n  Sampli ng interval  ms T s 4 10 5 , prop agatio n functio n  and  surge impe d ance coeffici en are  sho w n i n   Figure 3. The  comp ari s o n  result b e twe e n  the a c tual  simulation  re sults an d Fou r i e r   algorith m , sh own in Fig u re  4.    Evaluation Warning : The document was created with Spire.PDF for Python.
                               ISSN: 23 02-4 046   TELKOM NIKA  Vol. 11, No . 5, May 2013 : 2583 – 259 3   2588 mi n s T mi n mi n mi n N T T s N k a a T j L p j k e j A 1 mi n ) ( mi n s f s T K T min / 1 s s T N N N f s / 5 . 0 0 L a ) / ) 1 ( e xp( ) ( 2 1 ) ( 0 0 0 s f N k L s L N K n j f k A N n a ) ( 0 mi n j A e L T j 0 0 L L a a N N M M N k a a L z z z z p z k z A 1 1 1 1 0 1 0 1 ) (     Figure 2. Pro c ed ure of Obt a ining Th e Recu rsive  Coef ficients          Figure 3. Pro pagatio n Fun c tion an d Surge Impeda nce Coeffici ent   Evaluation Warning : The document was created with Spire.PDF for Python.
TELKOM NIKA   ISSN:  2302-4 046     The Re se arch of Digital Algorithm  Base on Fre que n c y-Depe nde n t  ... (Yongqin g  Liu)  2589     Figure 4. The  compa r i s on result s betwee n  discrete  se quen ce p r op a gation fun c tio n  with the  simulat i o n  re sult s       From  Figu re  4, on can  se e that  propag at ion  co efficient h a ve 40  poi nts  non-ze ro   element s, if li ne i s  lo nge r,  sampli ng f r eq uen cy is  hig h e r, the n  p r op agation  fun c tion  will b e  lo n ger  delay time, L ead to in crea sed  com putat ion. By  empl oying Pro n y algorith m  to reasona bly limit  the len g th of  se que nce, compa r ed  wit h  propa gatio n fun c tion  gi ven by PSCAD, are  sho w n i n   Figure 5.        Figure 5. The  compa r i s on result s Prony  algorith m       From figu re 5, it is concl uded that Prony  algorithm  is corre c t an d accurate, and have   less com puta t ion quantitie s.      3. The Propa gating Trav eling Wav e s   Suppo se that , at the location  x =0, the traveling wave  can be e s tabl ishe d by mea s ured  voltage and  current at  origin al termi nal M, forward a nd ba ckward t r aveli ng wave  dig i tal  expre ssi on a r e written a s :     ) ( * ) ( ) ( ) ( * ) ( ) ( k u Y k i k b k u Y k i k f m C m m m C m m  (28 )   Evaluation Warning : The document was created with Spire.PDF for Python.
                               ISSN: 23 02-4 046   TELKOM NIKA  Vol. 11, No . 5, May 2013 : 2583 – 259 3   2590   Whe r e,  ) ( k Y C is a d mittance  co efficient, Symbol “*” m e ans  revol u tio n  ope rato r.  Therefore the  solution of di fferential trav e ling wave eq uation is  sho w n in the follo ws.     ) ( * ) ( ) ( ) ( * ) ( ) ( k a k b k e k a k f k e L n m L m n  (29 )     The p r op agat ing relation s i n  discrete tim e  dom ain b e twee n the t w o  terminal s M   and  are sho w n in  the followin g   N l m l n M l l m N l n l M l m l n l k e N l k b k e l k e N l k f k e 1 min 0 1 0 min ) ( ) ( ) ( ) ( ) ( ) (  (30 )       4. Simulation Resul t s   This sectio n pre s ent s sim u lation re sult of the new algorithm compa r ed wit h  the   results of EMTP, asso ciate d  with real 5 0 0 kV  tran smi s sion n e two r in ShanDong  powe r  syste m sho w n in Fi g u re 6. Th e int e re sted transmissi on lin is that of ZiBo  Station to ZouXian Plant, the  line length i s  328 km, the line with freq u ency de pen d ent para m ete r s.   Suppo se that  the voltage s and  cu rre nts are  m e a s u r e d  at Zibo Sta t ion, the sa m  ling  perio d is 0.1 m s, that is, 200 sam p le s p e r cycl e.  The  forwa r d travel ing wave s at ZouXian term inal  and p ba ckward traveli ng  wave s at Zibo Station are  requi re d to ca lculate.           Figure 6. The  associ ated 5 00kV po we r n e twork      4.1. Compar ed  w i th th e Calcula t ed  Results   Suppo se th at a  singl e p h a se  to e a rth  fault is take n pla c e  in th e line  of Z o uxian  Station, the  forwa r d  trav eling  wave s at Zo uX ian  termin al a n d  ba ckward   traveling  wa ves,  cal c ulate d  by  the data  of Zibo Station  usin g t he n e w  alg o rithm,  comp ared  with the  wavefo rms  given by PSCAD, are sho w n in Figure 7 and 8.   From  Figu re   8 an d 9,  one   can  see th at  t he  re sult s ca lculate d   by n e w algo rithm have   relatively a c cura cy, compa r ed  with  the  result s by PS CAD the m a ximum erro r i s  n o  m o re  th an  1.04% u nde r the  sam p lin g pe rio d  0.1 m s, at  the  sa me time,  wit h  le ss am ou nt of  comp utation,  only have 10  times sum an d multiply.    Evaluation Warning : The document was created with Spire.PDF for Python.
TELKOM NIKA   ISSN:  2302-4 046     The Re se arch of Digital Algorithm  Base on Fre que n c y-Depe nde n t  ... (Yongqin g  Liu)  2591     Figure 7. The  waveform of  the cal c ulate d  forw a r d trav eling wave at terminal Zo u X ian Plant  comp ared wit h  the re sult           Figure 8. The  waveform of  the cal c ulate d  ba ckward traveling wave at terminal Zi bo Plant  comp ared wit h  the re sult       4.2. The Res u lts under V a rious Samp ling Periods   The wavefo rms of traveli ng wave s co m puted by the pro p o s ed  method, co mpared  with wavefo rms given by PSCAD, und er the sampli ng peri od of 0 . 01ms, are sh own in Fig u re  9.  Evaluation Warning : The document was created with Spire.PDF for Python.
                               ISSN: 23 02-4 046   TELKOM NIKA  Vol. 11, No . 5, May 2013 : 2583 – 259 3   2592     Figure 9. The  forwa r d trave ling wave  co mpari s o n  re sults at termin al ZouXian  station wh en  ms T S 01 . 0           Figure 10. Th e backward traveling wave comp ar i s io n result s at term inal Zibo Station wh en  ms T S 01 . 0       The m a ximu m erro rs com pare d   with th e results  give n by PSCA D   at variou s sa mpling  perio ds a r shown in Tabl e 1.      Table 1. Maxi mum errors u nder va riou sampli ng pe ri ods  Sampling  period (ms)   The for w ard t r aveling  w a v e  E rro r ( % )   The back w ard t r av eling  w a v e  E rro r ( % )   0.20 3.14  2.85  0.10 0.90  1.14  0.05 0.74  0.84  0.01 0.73  0.89  Evaluation Warning : The document was created with Spire.PDF for Python.