13pwrd02semlyen

Prévia do material em texto

596 IEEE. Transactions on Power Delivery, Vol. 13, No. 2, April 1998
COMBINED PHASE AND MODAL DOMAIN CALCULATION OF TRANSMISSION LINE TRANSIENTS
BASED ON VECTOR FITTING
Bjorn Gustavsen* (M) Adam Semiyen (LF)
Department of Electrical and Computer Engineering
University of Toronto
Toronto, Ontario, Canada M5S 3G4
* On leave from the Norwegian Electric Power Research Institute (EFI), Trondheim, Norway.
Abstract-This paper introduces two highly accurate transmission line
models. In the first one, particularly suitable for overhead lines, the
matrix transfer functions for both characteristic admittance and
propagation are fitted directly in the phase domain using the method of
vector fitting by optimal scaling. In the second model we use for
propagation a modal decomposition with a constant transformation
matrix and an optional phase domain correction term. Both models are
computationally highly efficient due to their time domain realizations
based on vector fitting.
1 INTRODUCTION
Transmission line models are needed in studies of
electromagnetic transients in power systems. Such models are
usually formulated in the time domain, as this facilitates easy
inclusion of the line model in general time domain simulation
programs, such as the EMTP. The model is then interfaced with
the host program by means of convolutions between impulse
responses characterizing the transmission line, and electrical
quantities at the line ends. The computational effort is greatly
reduced by the introduction of recursive convolutions [1], which
can be achieved by approximating the impulse responses in the
frequency domain by rational functions of low order.
Wedepohl [2] introduced modal decomposition to the theory of
wave propagation, thereby allowing the traveling waves to be
characterized by frequency dependent eigenvalues and
transformation matrices. The frequency dependence of the
transformation matrix T has traditionally been neglected [1,3]
because it is weak for many overhead lines of practical interest,
and because substantial savings in computational time is thus
achieved. However, this frequency dependence may be strong for
some line configurations (e.g. multi-circuit lines) and should
therefore be taken into account. This can be achieved by
introducing an additional convolution for T , as originally
suggested in [4] and later implemented in [5] and [6], Although
this approach seems to work fine with cable systems, recent
investigations [6] have shown that in the case of overhead lines it
may not always be possible to obtain an accurate rational function
approximation for T and for the characteristic admittance matrix
using stable poles only.
An alternative way of resolving the problem of frequency
dependent transformation matrices is to abandon the modal
approach and carry out the calculations directly in the phase
domain [7]. This method completely overcomes the problem of
PE-346-PWRD-0-01-1997 A paper recommended and approved
by the IEEE Transmission and Distribution Committee of the IEEE
Power Engineering Society for publication in the IEEE Transactions
on Power Delivery. Manuscript submitted July 26, 1996; made
available for printing January 8, 1997.
fitting T . An additional advantage is that the phase domain
transfer functions Yc and H are intrinsically stable, implying
that an accurate rational function approximation will automatically
be stable.
In the present paper we demonstrate that Yc and H can be
accurately fitted in the phase domain in the case of overhead lines.
Yc is in general very smooth and can easily be fitted with low
order rational functions. However, the elements of II have a
more complicated shape, which requires a higher order fit that
includes complex poles/zeros. This problem is overcome by the
application of a recently developed fitting methodology— vector
fitting by optimal scaling [6]. The new method produces a stable
and accurate fit for the elements of H , using both real and
complex poles/zeros. A procedure of optimal scaling is used to fit
H and Yc columnwise using the same set of poles, as this leads
to a more efficient time domain formulation [6] as compared to
element-by-element fitting. The method also works well for cable
systems as long as the modal velocities are not very different.
In the paper we also show that high accuracy can be achieved
for II by splitting it into the sum H =H0 +AH where H0 is the
propagation matrix calculated assuming a constant transformation
matrix. We then fit the modal propagation functions in H0 and
the remainder AH . In addition to giving increased accuracy for
overhead lines, this procedure makes the model applicable to
cable systems with widely different modal velocities.
2 METHODOLOGY
2.1 Model formulation
H ifar
compared to element-by-element fitting. In the case of Hm , the
D - matrix in (5) is 0. The time delays are removed prior to fitting
(unwinding) by multiplication with a factor exp(ycox) , where x
is calculated as shown in [6]. For Hm we use a separate x for
each element (mode). For H and AH we use a common x
which is chosen equal to the smallest of the modal x 's used in the
unwinding of Hm . The "starting poles" used in the fitting process
are in all calculated examples real and negative, logarithmically
distributed between 1 Hz and 10 MHz.
2.4 Time domain implementation
The combined modal and phase domain model defined by (1) and
(4) is transformed into the time domain by using convolutions :
Yc*v-i=2T10( Hm*i%r )+2AH*ifar (6)
where
ifar =T^ifar (7)
The convolution for the characteristic admittance in (6) can be
expressed as an instantaneous term plus a history term :
Yc*v=Yc0v+hl (8)
Because both Hm and AH contain a time delay greater than the
time step length used in the simulations, the right hand side in (6)
represents history terms only. Hence, insertion of (8) in (6) gives:
Yc0 v-i=h2 +h3 -hl (9)
This allows the transmission line model to be linked to the host
program by means of the Norton equivalent of figure 2.
h=h2+h3-hl Yco
•4-
+
V
Fig. 2 Network representation of line end
The instantaneous term Tc0 and the history terms hx , h2 and h3
are calculated from the respective SERs using trapezoidal
integration.
A similar representation is used for the full phase domain
model defined by (1).
3 CALCULATED RESULTS FOR OVERHEAD LINE
The transmission line models developed in the paper will in
the following be applied to the twin-circuit, untransposed
overhead line in figure 3.
A i 6.2m
6.7m Ground wires : d=1cm
Phase conductors : d=3.2cm
3 O 9-8m 0 6
6.1m
X 2 ° 14.6m 0 5
6.1m
y 1 O 9.8m O 4
14.5m
m m p=100Qm
Fig. 3 230 kV twin-circuit overhead line
3.1 Existing methods
The modal transformation matrix T can be strongly frequency
dependent in the case of untransposed overhead lines. This is
demonstrated in figure 4 for one eigenvector (ground mode) of the
twin-circuit line of figure 3.
598
1 , 4
2 , 5« 0 . 41
1
5 0.2 3 , 6
10z 104
Frequency [Hz]
Fig. 4 Magnitude Junctions of ground mode eigenvector
However, the currently used transmission line models neglect
the frequency dependence of T . It is assumed that the matrix
product YZ can be diagonalized using a real, constant
transformation matrix T10 , which is calculated at a user-specified
frequency. (Before discarding the imaginary parts of TI0 , its
columns are rotated so that the imaginary parts are minimized in
the least squares sense [9]).
This leads to the following relations for Yc and H :
Yc^ =TI0Y’ mTl0 (10)
H «H' =T10H' mTf0l (11)
Note that the modal propagation matrix H' m i n (l l) i s different
from the one in (4), as the diagonalization using TJ 0 leads to
inaccurate eigenvalues.
Although the assumption of a real, constant transformation
matrix often leads to satisfactory results, in other instances the
accuracy may be less good. As an example we show in figure 5 the
magnitude ftmctions of the first column of Yc as calculated by
(10), when using a T10 computed at 5 kHz. (5 kHz is the
fundamental frequency of a 15 km line.) Also is shown the
magnitude of the accurate Yc . It is seen that a noticeable error
occurs at low frequencies, including 60Hz.
xlO'3
2 rc5?
r;(5kHz)•8
M!i
IT
H
ir(60Hz)
 2
2 deviation
5
0
102 io4
Frequency [Hz]
Fig. 8 First column of H (15th order approximation)
599
2
'a
3ma
£
-200
-400
3,4
5Accurate
Fitted 6
10s102 104
Frequency [Hz]
Fig. 12 Phase angles of Hm (15th order approximation)
The correction term AH was calculated from (4) with 7}0
evaluated at 10 MHz. Figure 13 shows the magnitude functions of
the first column of All , and the resulting fit when using vector
fitting with 15 poles.
Accurate0.04
Fitted
£3
•3
5 3
1,0.02 4
deviation 62
2,5
10° 102 104
Frequency [Hz]
Fig. 13 First column of AH (15th order approximation)
Figure 13 shows that the calculated SER gives a good
approximation of the correction term. Including the fitted
correction term in the line model as defined by (4) will reduce the
magnitude of the deviation in the approximation for 77(1,1) from
0.04 to about 0.0005. This is about 20 times smaller than when
H was fitted directly (figures 8-10).
The phase angles corresponding to figure 13 are shown in
figure 14.
600
0 5
100
102 104
Frequency [Hz]
Fig.14 Phase angles of AH (15th order approximation)
3.4 Time domain results
In the following we show calculated results in the time domain for
the 100 km transmission line of figure 3. Using vector fitting the
frequency domain transfer responses were approximated by SERs
in the interval 1 Hz-10 MHz. We used 9 poles for each column of
Yc , 15 poles for each column of H , 15 poles for each column of
AH and 15 poles for each diagonal element of Hm . TI 0 was
calculated at 10 MHz.
3.4.1 Oven circuit response
In this test a 1 p.u. step voltage is applied to conductor #1, as
shown in figure 15. The other conductors are grounded at the
sending end, and all conductors are open circuited at the receiving
end. (The numbering of conductors is defined in fig. 3).
t — 0 1
2
3
4
5
lp.u. y ^ ) 6
100 km
v6
mmmmmmmmmmmmmmmmmmmimm.
Fig.15 Open circuit test
Figure 16 shows the open circuit voltage on conductors #1 and
#6, as calculated by the hybrid model (model #2 ). Also is shown
the theoretically accurate voltages by a Fourier transform. The
calculated responses are seen to match the Fourier solution very
accurately when the correction term AH is used. Neglecting AH
leads to a significant deviation for response #1.
2
1.5
I .
to,
£
o
-0.5
vi
vs
0 1
Fourier method
Hybrid with AH
Hybrid without AH
,/
2 3
Time [ms]
4 5
Fig.16 Open circuit responses
Very high accuracy was also obtained when using the full
phase domain method (model #1). Figure 17 shows the deviation
from the Fourier solution for response #1 for the following cases :
a) Full phase domain model
b) Hybrid model including AH
c) Hybrid model neglecting AH
0.1
a)
b)
0.05 c)
& 0
£ -0.05
-0.1
0 41 2 3
Time [ms]
Fig. 17 Deviation from Fourier solution (response #1)
5
It is seen that the deviation is mostly limited to about 0.01 for
models a) and b). The deviation is significantly higher near the
steep fronts of the response (figure 16). This is due to the fact that
the upper frequency limit used in the Fourier solution was limited
to 500 kHz (in order to limit the number of frequency samples to
be calculated), while the responses used in the simulations were
fitted up to 10 MHz. Thus, the simulated responses have steeper
fronts than those of the Fourier solution.
3.4.2 Short circuit response
In this test a 1 p.u. step voltage is applied to conductor #1, as
shown in figure 18. The other conductors are grounded at the
sending end, and all conductors are short circuited at the receiving
end.
f =0
1P -u.
100 km
>16
1
2
3
4
5
6
mmmmmmmmmmmmmmsmmmmmmmii
Fig. 18 Short circuit test
Figure 19 shows the short circuit current flowing in conductor
#6 (sending end), as calculated by the
• Fourier method
• Full phase domain model
• Hybrid model including AH
• Hybrid model neglecting AH
It is seen that the simulated responses match the response by the
Fourier method quite accurately. However, when the correction
term is neglected in the hybrid model, then the deviation becomes
very evident.
601
x 10 3
OrV
r— n
P “1
\r~
&
1, ^Fourier method l A
•— - Full phase domain model 1 ,
'l
-- -- Hybrid with AH ~‘‘ 1
Hybrid without AH vg
.3 I . •
0 1 2 3 4 5
Time [ms]
Fig.19 Short circuit response
4 DISCUSSION
The two simulation models developed in the paper have been
demonstrated to produce accurate results for overhead lines. This
result is due to the fact that the modal time delays are not very
different in the case of overhead lines. Accordingly, it is possible
to fit the propagation matrix H (model #1) or a remainder
AH (model #2) directly in the phase domain, using a fairly small
number of poles.
Results not shown in the paper have demonstrated that the two
models can also be applied to underground cables provided that
the modal velocities are not very different. For instance, the
sheaths of coaxial cables are normally grounded at both ends. One
may therefore eliminate the sheath conductors from the system of
conductors, which leads to modes having equal time delays, and
thus smooth functions for H and AH . However, in the case of
crossbonded cables the velocities will be widely different as the
sheaths cannot be neglected. Therefore, an accurate fitting of H
or AH would require a very high order fit. But model #2 could
still be used if one accepted to neglect AH , which can be justified
since AH is usually quite small. In view of these results it
appears that model #2 is somewhat more general than model #1.
The two models introduced in this paper are expected to be
significantly faster in a time domain simulation than the one
introduced in [7]. This is because the usage of complex poles leads
to a lower order fit when approximating H , and because the
usage of columnwise realization for H and Yc gives
approximately a 3-fold increase in efficiency as compared to an
element-by-element realization as used in [7],
5 CONCLUSIONS
In this paper we have developed two models for simulation of
electromagnetic transients on overhead lines which produce
accurate results in the presence of a frequency dependent
transformation matrix.
Model #1 is a fully phase domain model in which the matrix
transfer functions for the characteristic admittance Yc and the
propagation H are fitted with State Equation Realizations (SER)
directly in the phase domain. The columns of Yc can be accurately
fitted with as little as 5-10 poles per column. A higher order fit is
needed for H since the modal contributions have in general
different time delays. An order of 15 will normally produce an
accurate fit.
Model #2 splits the propagation matrix H into two
terms : Tj0HmT[0 and AH , where 7}0 is a constant
transformation matrix evaluateU at High frequencies, H" is a
diagonal matrix, and AH is a correction term in the phase
domain. This model appears very attractive because AH is
usually quite small and could therefore be considered an optional
correction term when high accuracy is desired.
The matrix transfer functions H and AH will in general
contain oscillating components due to differences in the modal
time delays. The success achieved in fitting these functions is
mainly due to the adopted fitting methodology— vector fitting
using optimal scaling. This method automatically produces the
combination of real and complex poles/zeros needed to fit the
complicated shapes of H and AH .
The application of vector fitting has the additional advantage
that the resulting time domain implementation is about 3 times
more efficient than an implementation based on element-by¬
element fitting, as used in [7],
Further investigations, not shown in the paper, have
demonstrated that both models #1 and #2 can be applied to cable
systems, provided the modal velocities are not widely different.
Model #2 can even be used in such cases if AH is neglected.
6 ACKNOWLEDGEMENTS
Financial assistance by the Natural Sciences and Engineering
Research Council of Canada is gratefully acknowledged. The first
author wishes to express his gratitude to the Norwegian Electric
Power Research Institute (EFI), Trondheim, Norway, for granting
and financing his leave at the University of Toronto.
7 REFERENCES
[1] A. Semiyen and A. Dabuleanu, "Fast and Accurate Switching Transient
Calculations on Transmission Lines With Ground Return Using
Recursive Convolutions", IEEE Trans. PAS, vol. 94, March/April 1975,
pp. 561-571.
[2] L. M. Wedepohl, "Application of matrix methods to the solution of
travelling-wave phenomena in polyphase systems", Proc. IEE, vol. 110,
no. 12, December 1963, pp. 2200-2212.
[3] J. R. Marti, "Accurate Modelling of Frequency-Dependent Transmission
Lines in Electromagnetic Transient Simulations", IEEE Trans. PAS, vol.
101, no. 1, January 1982, pp. 147-157.
[4] A Ametani, "Refraction Coefficient Method for Switching-Surge
Calculations on Untransposed Transmission Lines (Accurate and
Approximate Inclusion of Frequency Dependency)", IEEE PES Summer
Meeting, C 73-444-7, 1973.
[5] L. Marti, "Simulation of Transients in Underground Cables With
Frequency-DependentTransformation Matrices", IEEE Trans. PWRD,
vol. 3, no. 3, July 1988, pp. 1099-1110.
[6] B. Gustavsen and A Semiyen, "Simulation of Transmission Line
Transients Using Vector Fitting and Modal Decomposition", accepted for
publication in IEEE Transactions and presentation at the IEEE PES
Winter Meeting, New York, 1997.
[7] H. V. Nguyen, H. W. Dommel, and J. R. Marti, "Direct Phase-Domain
Modelling of Frequency-Dependent Overhead Transmission Lines", paper
no. 96 SM 458-0-PWRD.
[8] L. M. Wedepohl. H. V. Nguyen and G. D. Irwin, "Frequency-Dependent
Transformation Matrices for Untransposed Transmission Lines Using
Newton-Raphson Method", paper 95 SM 602-3 PWRS, presented at
IEEE Summer Meeting, Portland, 1995.
[9] V. Brandwajn, "Modifications of user's instructions for "MARTI
SETUP"", EMTP Newslwetter, vol. 3, no.l, August 1982, pp. 76-80.
8 BIOGRAPHIES
Bjorn Gustavsen was bom in 1965 in Harstad, Norway. He received the
M.Sc. degree in Electrical Engineering from The Norwegian Institute of
Technology, Trondheim, in 1989, and the Dr. Ing. Degree in 1993. Since
1994 he has been working at the Norwegian Electric Power Research Institute,
mainly in the field of transient studies. He is currently on leave at the
Department of Electrical and Computer Engineering, University of Toronto.
Ajam Semiyen was bom in 1923 in Rumania where he obtained 9 Dipl Ing.
degree and his Ph.D. He started his career there with an electric power utility
and held academic positions at the Polytechnic Institute of Timisoara. In 1969
he joined the University of Toronto where he is a professor in the Department
of Electrical and computer engineering, emeritus since 1988. His research
interests include steady state and dynamic analysis as well as computation of
electromagnetic transients in power systems.
602
Discussion
Brandão Faria, Senior Member, IEEE (CETME, Instituto Superior
Técnico, Av. Rovisco Pais, 1096 Lisboa Codex, Portugal):
This paper, together with a companion one [1], presents an impor¬
tant and lasting value contribution to the field of transmission line
transients computation. The authors are, thus, to be commended by
the excellent work they have done. The few comments and ques¬
tions that next follow have the sole purpose of permitting the clarifi¬
cation of some minor details.
1) The so-called full phase domain model (model #1) referred to by
the authors is based on the fit of both the propagation matrix func¬
tion for currents H = exp(-/YZ 1) and the characteristic admittance
matrix Yc = Z-1 JZY . Do the authors mean that the matrix func¬
tions H and Yc can be evaluated without resorting previously to the
modal domain? To determine JZY (OTJYZ ) one must go through
the computation [2] of
JzV = Tv Y T[ and JYZ = Ti y T^where y is the diagonal matrix of the modal propagation constants
and Ti and Tv are the right eigenvector matrices of YZ and ZY re¬
spectively. How can one then avoid the modal domain?
2) In Fig.10 several curves are presented with the purpose of show¬
ing the increasing of fit accuracy with the increasing of the number
of poles. Rather strangely the plots reveal that, for low frequencies,
a fit with 10 poles yields better results than a fit with 15 or even 20
poles. Can it be found a justification for this?
3) The frequency domain analysis carried out by the authors extends
from 1Hz to 10MHz; this upper frequency corresponding to propa¬
gation modes with wavelengths around A » 30m . The overhead line
configuration used for simulation purposes depicted in Fig.3 shows
conductor heights in the range 14.5mof these models. As the
authors pointed out, the idea of fitting Yc and H functions directly
in the phase domain was also shared by the discusser and his col¬
leagues [7], It is clear that the advantage of fitting these functions
in the phase domain is the reduction in the convolution operations
required in time domain simulation and that vector fitting by opti¬
mal scaling using both real and complex poles/zeros with pole
sharing would lead to more efficient in time-domain formulation.
The fitting methodology utilized in[7], however, has the advantage
that it uses strictly the minimum-phase rational functions and thus
guarantees the stability of the solution. The solution stability as¬
pect of the phase-domain transmission line model in the paper was
addressed and confirmed, as the authors indicated, ". . .the phase do¬
main transfer functions Yc and H are intrinsically stable, implying
that an accurate rational function approximation will automatically
be stable." However, the discusser feels that additional time-
domain testing for different circuit configurations, especially for
those having asymmetrical geometries, with long-duration simula¬
tion (3 to 4 seconds), may be needed in order to strengthen this as¬
pect of the model. The significant improvement may even be
highlighted further if such a model can be mathematically proven
for its stability.
It would be helpful for the discusser and the readers if the authors
could describe in more detail the time-domain implementation pro¬
cedure which utilizes recursive convolutions for functions having
complex poles/zeros.
The phase angles of the off-diagonal elements of 1st column of H
showed in Figure 9 seem to start out at a zero angle value. This
characteristic is opposite of that obtained in [7], The different con¬
ductance values used in the line admittance matrix for the models
may attribute to this difference. A value of 3.0 x 10'9 mho/km was
normally used in [7], Did the author assume a zero conductance
value? If not, what was the typical value used for it?
The authors' comments on the above points are highly appreciated.
Manuscript received March 5, 1997.
Bjorn Gustavsen and Adam Semiyen: We wish to
thank the discussers for their valuable comments and
useful contributions. Our answers are as follows.
Professor Brandão Faria :
Ad 1 The designation as full phase domain refers only
to the time domain simulations. As properly emphasized
by Professor Brandão Faria, in the preparatory phase,
the frequency domain evaluation of matrix functions for
phase variables is performed in the conventional way by
use of similarity transformations. These phase domain
matrices are the object of the state equation
approximation obtained by vector fitting. Their ultimate
use, in the context of this paper, is for direct phase
domain simulation of transmission line transients.
Ad 2 The very high accuracy at low frequencies
achieved with 10 poles is very local. The overall
impression from figure 10 is that increasing the number
of poles increases the accuracy of the fitting.
Ad 3 The discusser is absolutely right that at higher
frequencies the validity of a methodology based on pure
longitudinal wave propagation is increasingly reduced. A
remark to this effect should have been made in the paper.
Ad 4 We also agree with Professor Brandão Faria that
at higher frequencies the double circuit transmission line
should be represented as a system of eight conductors.
Then, discrete two-port results could be obtained for the
individual spans so that an equivalent six-conductor line
may be conceived with the characteristics, including
resonances, of the actual line. We did not however
examine the applicability of vector fitting for such
problems.
Dr. Taku Noda and Professor Akihiro Ametani :
The switching-back procedure can only be used for
eliminating artificial eigenvector switchovers. By this we
mean the abrupt switchovers that occur when the
sequence of eigenvalues returned by the diagonalization
routine changes between two successive frequency
points. On the other hand, the intrinsic swithcovers
referred to by the discussers take place as smooth
variations as function of frequency and will have to be
fitted.
For each line end, the full phase domain model requires
In2 convolutions. The combined model requires
n + 2«2 or n+n2 convolutions when, respectively, the
correction term AH is included or neglected (n2
convolutions come from the characteristic admittance).
Thus, the combined model is slightly less efficient than
the full phase domain model when including AH , and
more efficient when neglecting AH .
Dr. Huven V. Nauven :
We agree with the discussor that fitting with stable poles
will in general not guarantee a stable simulation. This
property is, however, assured when fitting with stable,
minimum phase-shift functions, as in the methodology
adopted by the discussor in [7], So far we have not
observed any stability problems. The full phase domain
model has also been tested with trapped charge
604
simulation, without any sign of instability. When fitting
with non-minimum phase-shift functions (as was done in
this paper), the only way of guaranteeing a stable model
is to check that the line admittance matrix has
eigenvalues with positive real parts for all frequencies. It
is also possible to check the eigenvalues of Yc and H
separately, as shown in [A] If, however, the eigenvalues
suggest the possibility of an unstable simulation, it is not
clear how the model should be changed ’ so as to
guarantee stability.
Complex poles can be handled in the time domain
implementation in the same way as the real poles. The
only difference is that some poles and residues (elements
in the C-matrix) are now complex numbers, which has
the implication that the state variables x become complex
numbers. For further details regarding the
implementation, we refer to the closure of [6].
When extending the frequency range in figure 9 to lower
frequencies, the phase angles of the diagonal elements
approached zero while those of the off-diagonal elements
approached +45 degrees. We used a zero value for the
shunt conductance G.
[A] T. Noda, N. Nagaoka, and A. Ametani, “ Further
Improvements to a Phase-Domain ARMA Line
Model in Terms of Convolution, Steady-State
Initialization, and Stability” , IEEE Trans, on Power
Deliveiy, Vol. 12, No. 3, July 1997, pp. 1327-1334.
Manuscript received September 9, 1997.