A benchmark test system for testing frequency dependent network equivalents for electromagnetic simulations

Electrical Power and Energy Systems 44 (2013) 364–374

Contents lists available at SciVerse ScienceDirect

Electrical Power and Energy Systems journal homepage: www.elsevier.com/locate/ijepes

A benchmark test system for testing frequency dependent network equivalents for electromagnetic simulations Y.P. Wang a,⇑, N.R. Watson b a b

Electrical Industry Research Institute of Korea, 533-2, Deungchon-2Dong, Gangseo-Gu, Seoul 157-718, Republic of Korea Department of Electrical and Computer Engineering, University of Canterbury, Private Bag 4800, Christchurch, New Zealand

a r t i c l e

i n f o

Article history: Received 13 December 2010 Received in revised form 4 April 2012 Accepted 9 April 2012 Available online 26 September 2012 Keywords: Electromagnetic transient Equivalents AC system representation

a b s t r a c t The complexity of modern power systems often makes it impractical to model it in its entirety for electromagnetic transient studies. Therefore areas outside the immediate area of interest must be represented by some form of frequency dependent network equivalent (FDNE). The equivalent must be ‘‘frequency dependent’’, that is match the frequency response of the system it represents, if it is to faithfully represent it. This paper presents a z-domain rational function formulation for developing a FDNE and demonstrates the use of it for the assessment of the steady-state and transient response of the Lower South Island of New Zealand. This paper also shows the accuracy that can be expected from using FDNEs as well as providing a test system for other researchers to benchmark their work against. This is achieved by using a well publicized test system (Lower South Island of New Zealand) the parameters of which are readily available as well as providing complete information of the developed FDNE. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction Although time domain techniques, such as used in EMTP and PSCAD/EMTDC packages, can accurately perform analysis of power system transients, detailed representation of a large complex power system will entail a prohibitive amount of computation. Therefore there is a need to have equivalents which adequately represent the areas outside the immediate area of interest. These equivalents must mimic the frequency response of the network it represents in order to correctly mimic the transient behaviour, hence these equivalents are called frequency dependent network equivalents (FDNEs). Conventional equivalents (Thevenin) based on the fundamental frequency short-circuit level are inadequate for representing the external network’s behaviour due to the presence of other frequency components. Early FDNEs modelled the external system by an appropriate network of lumped R, L and C components whose values are chosen so the equivalent network will have the same frequency response as the external network [1]. Using RLC components allows implementation in existing transient programs with minimal change, however it restricts the possible frequency response that can be represented. The use of rational functions in s (or z) domain is more general and this has been the dominant approach in recent years.

⇑ Corresponding author. Tel.: +82 2 3219 0692; fax: +82 2 3219 0559. E-mail addresses: [email protected] (Y.P. Wang), [email protected] (N.R. Watson). 0142-0615/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijepes.2012.04.065

Fitting in the s-domain always requires discretizing a continuous system and the inherent approximation when simulating the rational function. The work of Abur and Singh [2,3], Noda et al. [4,5], Gustavsen and Semlyen [6], Angelidis and Semlyen [7], and Morched et al. [8] are particularly noteworthy. The use of z-domain FDNEs was reported in [9]. One of the motivations for investigating z-domain fitting is that it can be directly implemented in a digital simulation program without any loss of accuracy as it is already a discrete formulation. Parameters such as the frequency range and weighting factors used for fitting do influence whether unstable poles are generated. Moreover a three-phase FDNE can be unstable even if all the self and mutual terms are fitted by stable rational functions, if the 3  3 matrices are not positive definite at all frequencies. This problem has been tackled by Gustavsen and Semlyen [10] for s-domain fitting. Most work on FDNE has used s-domain rational function and often in partial fraction form, whereas the present work uses a rational function in the z-domain. Ref. [9] introduced the concept of z-domain FDNE and demonstrated its use on simple single-phase systems (both single and two-port). A recent contribution outlines the least-squares fitting technique and applies it to a practical test system (Lower South Island of New Zealand) [11]. Although it used a full 3  3 representation (including mutual coupling), it is clear from the response that this FDNE is not optimal. A better FDNE has recently been developed and used specifically for harmonic analysis [12]. The paper compared various FDNEs (including the use of RLC networks) for harmonic assessment on a test system which includes an HVDC link.

365

Y.P. Wang, N.R. Watson / Electrical Power and Energy Systems 44 (2013) 364–374

f

BUS 2 14 kV

BUS 1 14 kV

ROXBURGH 220 kV

MANAPOURI 220 kV 11 kV

220 kV INVERCARGILL

TIWAI

33 kV

220 kV Boundary Busbar for FDNE

Fig. 1. Lower South Island of New Zealand test system.

1

2

Voc1

3

Voc 3

Voc 2

− y12

− y23

− y13 yself 1

yself 2

yself 3

Fig. 2. Three-phase FDNE.

This paper demonstrates the accuracy achievable for transient analysis of a practical system by applying the enhanced FDNE to the Lower South Island of New Zealand. One of the objectives of this paper is to provide a benchmark for this type of work just as Ref. [13] is a benchmark for harmonic analysis. The Lower South Island was an actual system and all the input data, such as line geometry, and transformer parameters is available in published literature [14]. The inclusion of all the FDNE coefficients allows the transient response of this FDNE to be assessed easily and compared with those generated by other methods. The proposed method, although not as elaborate as [15,16], is straightforward to implement and has proven to give good and stable fits, and is a practical method for engineers to use.

Fig. 3. Norton representation of a rational function.

Table 1 FDNE voltage source parameters. Voltage (phase to neutral)

Phase A Phase B Phase C

Magnitude (kV)

Phase angle (°)

121.599 121.599 121.599

61.16755 59.15545 180.66100

2. Derivation of frequency response The development of a FDNE requires knowledge of the frequency response of the system to be represented. The modelling

366

Y.P. Wang, N.R. Watson / Electrical Power and Energy Systems 44 (2013) 364–374

300 Va

200

Vb

Terminal Voltage (kV)

100 Vc

0 −100 −200 −300 −400 0.8

0.81

0.82

0.83

0.84

0.85

0.86

0.87

0.88

0.89

0.9

Time (sec.) Fig. 4. Comparison of transient assessed using full system and FDNE representation, solid-full system, dotted-FDNE.

200 Va

Vb

Terminal Voltage (kV)

150 100

Vc

50 0 −50 −100 −150 −200 0.81

0.815

0.82

0.825

0.83

0.835

0.84

0.845

Time (sec.) Fig. 5. Fault inception, solid-full system, dotted-FDNE.

200

Va Vb

Vc

Terminal Voltage (kV)

150 100 50 0 −50 −100 −150 −200 0.86

0.865

0.87

0.875

Time (sec.) Fig. 6. Fault removal, solid-full system, dotted-FDNE.

0.88

0.885

367

Y.P. Wang, N.R. Watson / Electrical Power and Energy Systems 44 (2013) 364–374

0.04 Full 5−1250 5−2500 5−500

I (kA)

0.02

A

0 −0.02 −0.04 0.8

0.81

0.82

0.83

0.84

0.85

0.86

0.87

0.88

0.89

0.9

0.81

0.82

0.83

0.84

0.85

0.86

0.87

0.88

0.89

0.9

0.81

0.82

0.83

0.84

0.85

0.86

0.87

0.88

0.89

0.9

0.04

IB (kA)

0.02 0 −0.02 −0.04 0.8

0

C

I (kA)

5

−5 0.8

Time (sec.) Fig. 7. Comparison of current waveforms for both fault application and removal.

Table 2 RL values for 3  3 FDNE.

Table 5 Coefficients of rational function representing Yself1 term.

Term

Ohms

Henries

Order

a

b

Yself1 Yself2 Yself3 Y12 Y13 Y23

11.842 11.879 11.842 19.514 19.644 19.515

0.19647 0.20655 0.19647 0.34898 0.39479 0.34899

0 1 2 3 4 5 6

3.0078294923663079e03 1.6867157879437031e02 3.9394334773642417e02 4.9025038511999172e02 3.4265228069600508e02 1.2743659396358313e02 1.9684664251101679e03

1.0 5.8699942174618052e+00 1.4405484110295935e+01 1.8918195083417309e+01 1.4022321525318672e+01 5.5620099984815141e+00 9.2239367111043769e01

Table 3 Coefficients for rational function representing RL. Term

Yself1 Yself2 Yself3 Y12 Y13 Y23

Table 6 Coefficients of rational function representing Yself2 term.

Order a0 (=a1)

b1

Order

a

b

1.2705353612953927e04 1.2086071925372443e04 1.2705353612953927e04 7.1537366252693527e05 6.3246286117103626e05 7.1534755297287628e05

9.9699092078614049e01 9.9712853546340041e01 9.9699092078614049e01 9.9720809601966065e01 9.9751517353535690e01 9.9720806584012944e01

0 1 2 3 4 5 6

2.4466368842016076e03 1.3673558091512085e02 3.1809746270842362e02 3.9405381653974635e02 2.7395429619140748e02 1.0125536046686263e02 1.5526659700345685e03

1.0 5.8774107818227614e+00 1.4441679174219205e+01 1.8988979106062335e+01 1.4091654909120685e+01 5.5960216493182253e+00 9.2907746088487098e01

Table 4 Rational function order for various fitting frequency ranges.

Table 7 Coefficients of rational function representing Yself3 term.

Term

Order (5–500)

Order (5–1250)

Order (5–2500)

Order

a

b

Yself1 Yself2 Yself3 Y12 Y13 Y23

6 6 6 8 7 8

7 12 7 11 11 11

15 15 15 17 17 17

0 1 2 3 4 5 6

3.0078294923663079e03 1.6867157879437031e02 3.9394334773642417e02 4.9025038511999172e02 3.4265228069600508e02 1.2743659396358313e02 1.9684664251101679e03

1.0 5.8699942174618052e+00 1.4405484110295935e+01 1.8918195083417309e+01 1.4022321525318672e+01 5.5620099984815141e+00 9.2239367111043769e01

of the frequency dependence of overhead lines and cables is well advanced in electromagnetic transient packages, but this is not the case with other system components; for instance the standard

models of generators, transformers and loads do not represent the increase in resistance (and slight reduction in inductance) associated with skin effect. Therefore the use of a frequency domain pro-

368

Y.P. Wang, N.R. Watson / Electrical Power and Energy Systems 44 (2013) 364–374

gram will give a frequency response closer to reality. However, if the FDNE is developed from the frequency response obtained from a frequency domain program its accuracy can only be assessed by corresponding measurements in the real system, which is not a practical proposition. Hence for verification purposes the approach taken in this paper is the derivation of the frequency response from the less accurate time domain simulation, as this allows the complete time domain model to be used as the benchmark. The New Zealand Lower South Island, shown in Fig. 1, is used as a test system. As the main interest is the operation of an aluminium smelter connected to Tiwai busbar, this is the terminal busbar for the FDNE and FDNE represents the AC system as seen from this busbar. To derive the self and mutual impedances at each frequency three tests were performed by injecting current in one phase at a time and with the fundamental frequency voltage sources removed. The resulting 3  3 impedance matrices are then inverted to provide the admittance matrices to be fitted. The off-diagonal elements of this matrix are the negative of the admittances of the interconnecting branches to be fitted. The diagonal elements (which are the sum of all the branch admittances connected to the node) allow calculation of the shunt terms (Yself) to be fitted. As the admittance matrix is symmetrical (i.e. Y12 is equal to Y21 and Y13 is equal to Y31) only six terms need to be fitted.

Table 8 Coefficients of rational function representing Y12 term. Order

a

b

0 1 2 3 4 5 6 7 8

6.6876629787707654e03 5.1491335355811199e02 1.7424585227718672e01 3.3855097683174190e01 4.1315377108245688e01 3.2433392314869508e01 1.5996468074640230e01 4.5324861745415362e02 5.6491299976239055e03

1.0 7.1955469088034087e+00 2.2489101993989813e+01 3.9786978675039855e+01 4.3445034023350161e+01 2.9847753735867922e+01 1.2512802920585504e+01 2.8933709408255108e+00 2.7671132262133130e01

Table 9 Coefficients of rational function representing Y13 term. Order

a

b

0 1 2 3 4 5 6 7

1.7399872694598111e03 1.1263856073053968e02 3.1262831838239538e02 4.8223950008405743e02 4.4649362195819012e02 2.4814410002616423e02 7.6656127060372083e03 1.0155778304361271e03

1.0 6.7774986009704321e+00 1.9725734762103109e+01 3.1956619989073147e+01 3.1120004223321381e+01 1.8215063623245992e+01 5.9329131941526958e+00 8.2946996474582613e01

Table 10 Coefficients of rational function representing Y23 term.

3. Structure of the frequency dependent network equivalents

Order

a

b

0 1 2 3 4 5 6 7 8

6.6879759169563204e03 5.1491829884128787e02 1.7424233072110848e01 3.3853666697027623e01 4.1313051174614734e01 3.2431375927022815e01 1.5995509304323929e01 4.5322602381628350e02 5.6489470797505929e03

1.0 7.1943188253680557e+00 2.2480775327153474e+01 3.9762719625266101e+01 4.3405670023731673e+01 2.9809336292797784e+01 1.2490254511167217e+01 2.8860022832352925e+00 2.7567716462752406e01

Current Error (kA)

x 10

Current Error (kA)

Pn k IðzÞ a0 þ a1 z1 þ a2 z2 þ    an zn k¼0 ðak  z Þ P ¼ n 1 2 n VðzÞ 1 þ b1 z þ b2 z þ    bn z 1 þ k¼1 ðbk  zk Þ

5 0 −5

0.8

0.82

0.84

0.86

10

0.88

0.9

0.92

0.94

0.96

0.98

1

Phase B

x 10

5−1250 5−2500 5−500 RL

5 0 −5

0.8

0.82

0.84

0.86

0.88

0.9

0.92

0.94

0.96

0.98

1

Phase C

2

5−1250 5−2500 5−500 RL

RL 5−2500 5−1250

1 0

5−500

−1

0.8

0.82

0.84

0.86

ð1Þ

The a and b coefficients are calculated using the algorithm detailed in Appendix A so as to give a rational function that matches the given frequency response.

5−1250 5−2500 5−500 RL

−3

Current Error (kA)

HðzÞ ¼

Phase A

−3

10

The three-phase system equivalent of the test system will be of the form shown in Fig. 2, where each frequency dependent block is represented by a rational function of the form:

0.88

0.9

0.92

Time (sec.) Fig. 8. Comparison of current error.

0.94

0.96

0.98

1

Y.P. Wang, N.R. Watson / Electrical Power and Energy Systems 44 (2013) 364–374

369

Fig. 9. Expanded view of current in unfaulted phases (phase A and B) for fault application.

The rational function (1) can be implemented exactly as a Norton equivalent by identifying the instantaneous and history terms. From Eq. (1)

IðzÞ ¼

a0 VðzÞ |fflfflffl{zfflfflffl}

Instantaneous Term

þ ða1 z1 þ a2 z2 þ    an zn ÞVðzÞ  ðb1 z1 þ b2 z2 þ    bn zn ÞiðzÞ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} Histrory Term

¼ Gequiv VðzÞ þ IHistrory

ð2Þ

Hence in discrete time:

ik ¼

a0 mk |ffl{zffl}

Instantaneous

þ a1 mk1 þ a2 mk2 þ    an mkn  ½b1 ik1 þ b2 ik2 þ    bn ikn  |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} History

¼ Gequiv mðnDtÞ þ IHistory

ð3Þ

where mkm ¼ mðkDt  mDtÞ. Hence the rational function is implemented in PSCAD/EMTDC as a Norton equivalent (a current source with a parallel resistance as shown in Fig. 3). The resistance represents the instantaneous term (1/a0) while the current source represents all the past history terms

in current (a1, a2. . .an) and voltage (b1, b2. . .bn). The complete formulation and implementation for this approach is given in Ref. [13]. Although increasing the order improves the fit, some poles are sometimes unstable. Normally the highest order fit which is still stable is used. The least-squares approach to solving (2) is biased to give a better fit at the higher frequencies compared to the lower due to multiplying both sides of (A-3) by the denominator. In Refs. [14,15] an excellent iterative method of using adaptive weighting and QR decomposition is described to counter this bias. QR decomposition was selected due to its computational efficiency, because it must be used repetitively. An initial singular value decomposition (SVD) solution is used to determine the model order and hence remove redundant unknowns. In the present work a simpler approach of inverse frequency weighting is used so that the method is direct and not iterative. SVD is used due to its computational accuracy and robustness. The common practice of preconditioning the matrix before solving has been used. This technique has proven to give accurate and stable fits. The voltage sources (magnitude and phase angle) for all the FDNEs are set to be that of the full system under open circuit. These values are given in Table 1.

Fig. 10. Comparison of terminal voltage for fault application.

370

Y.P. Wang, N.R. Watson / Electrical Power and Energy Systems 44 (2013) 364–374

4. Simulation

4.1. Transient tests

The Lower South Island of New Zealand was chosen as the test system for the following reasons:

4.1.1. Base-case (5–1250 Hz fitting) The first test was to apply a single-phase fault to phase C at Tiwai. This is a more demanding fault as any error in modelling the mutual terms will be more evident in the unfaulted phase voltages and currents than if a three-phase fault was applied. The fault impedance was 0.001 X, and was applied to Tiwai phase C busbar at 0.815 s and removed after 0.05 s (arbitrary fault clearance time to allow evaluation transient due to fault removal). The comparison of the Full system and FDNE (using a frequency range of 5–1250 Hz for fitting) solutions is displayed in Fig. 4 and the coefficients used are given in [12]. The period immediately after fault inception (expanded in view given in Fig. 5) is practically identical in both cases, while differences are noticeable after fault removal (expanded view in Fig. 6). The phase current is shown in Fig. 7. Close inspection of these shows that the FDNE equivalent mimics the full system well; however, there is a low frequency component that causes a growing phase discrepancy in the terminal voltages (for the two unfaulted phases) during the fault period. This means that when the fault is removed there

 It is because it is a practical system. It was the system in use in New Zealand for this region for many years.  It is a demanding test system as it has untransposed lines and hence unbalanced. It can show the inaccuracies that can be expected in using a positive sequence (balanced) representation compared to full three-phase for a practical system.  The parameters, including transmission line geometry, are all available in previous publications. This allows electrical parameters of overhead transmission lines to be calculated, including all the frequency dependent self and mutual couplings.  The system has a large aluminium smelter fed from Tiwai 220 kV busbar and hence for studies of this system an equivalent of the rest of the system is desirable. All simulations were performed using PSCAD/EMTDC using a simulation time-step of 50 ls.

Fig. 11. Comparison of terminal voltage for fault removal.

Full 5−1250 1−1250

0.02 0.015

Current (kA)

0.01 0.005 0 −0.005 −0.01 −0.015 −0.02 0.814

0.816

0.818

0.82

0.822

0.824

0.826

Time (sec.) Fig. 12. Expanded view current in unfaulted phases (phase A and B).

0.828

0.83

371

Y.P. Wang, N.R. Watson / Electrical Power and Energy Systems 44 (2013) 364–374

Current (kA)

Phase A 0.02

Full 5−1250 1−1250

0.01 0 −0.01 −0.02 0.8

0.81

0.82

0.83

0.84

0.85

0.86

0.87

0.88

0.89

0.9

Current (kA)

Phase B 0.02

Full 5−1250 1−1250

0.01 0 −0.01 −0.02 0.8

0.81

0.82

0.83

0.84

0.85

0.86

0.87

0.88

0.89

0.9

Current (kA)

Phase C 6

Full 5−1250 1−1250

4 2 0 −2 −4 0.8

0.81

0.82

0.83

0.84

0.85

0.86

0.87

0.88

0.89

0.9

Time (sec.) Fig. 13. Comparison of the current waveforms.

Table 11 Injected harmonic current. Order

Phase

Magnitude (kV)

Phase angle (°)

11

A B C

0.07620 0.07620 0.07620

45.43 74.63 165.42

13

A B C

0.04790 0.04790 0.04790

40.60 160.62 79.44

23

A B C

0.01230 0.01230 0.01230

124.07 4.18 115.90

25

A B C

0.01297 0.01297 0.01297

123.60 116.40 3.52

is a slight discrepancy in the voltage between the Full and FDNE solutions which accounts for the slight differences in the transient on fault removal. This phase difference vanishes within approximately 0.25 s of fault removal. 4.1.2. Influence of fitting frequency range To investigate the influence of the fitted frequency range on the transient performance, FDNEs were also developed using fitting over the 5–2500 Hz and 5–500 Hz frequency ranges. Also for comparison a Thevenin equivalent (where each block in Fig. 2 is represented by an RL branch) based on the 50 Hz admittance was simulated. The RL is implemented using the appropriate rational 1 Dt function where a0 ¼ a1 ¼ Rþ2L= and b1 ¼ R2L . Table 2 gives the Rþ2LDt Dt R and L values used to model each of the admittance terms while Table 3 displays the coefficients used to implement this RL branch. In these comparisons the fault duration was not specified as removal occurs on the first current zero after 0.025 s following fault application. Table 4 shows the order of the rational function that was used for each term, for the various frequency ranges. Tables 5–10 give

the coefficients used in PSCAD/EMTDC simulation for the 5– 500 Hz frequency range. Fig. 7 shows a comparison of the phase currents (without the RL equivalent) for both fault inception and removal. Although some differences exist they all give a reasonable representation of the transient response of the system. Careful inspection of the current error, displayed in Fig. 8, shows that the greater the frequency range used for the fitting, the smaller the initial error; however, after 0.02 s all equivalents, even the 50 Hz based RL equivalent, give the same error. This error is due to the low frequency component causing the phase shift. Inspection of the error in the terminal voltage shows a very similar trend to the current error. To illustrate this Fig. 9 shows an expanded view of the current in the unfaulted phases. Close inspection of Fig. 9 shows that the RL equivalent gives the worst representation, followed by 5–500 Hz matched equivalent then 5–1250 Hz and 5–2500 Hz being the best. This demonstrates that the greater the frequency range of the matching process, the better the FDNE mimics the system response. Figs. 10 and 11 show an expanded view of the terminal voltage on fault application and removal. All except the 50 Hz based RL equivalent give a good representation of the full system, with the main error at fault removal due to the voltage phase shift. 4.1.3. Importance of fitting at low frequencies The comparison of the full system, FDNE (using a frequency range of 5–1250 Hz with 5 Hz spacing) and FDNE (using a frequency range of 1–1250 Hz with 1 Hz spacing) solutions is displayed in Figs. 12 and 13. These clearly show a greatly improved accuracy and the slow low frequency phase-shift at fault inception has also disappeared. This has demonstrated the importance of accurately fitting and low frequency, even below the fundamental. The simulation time for the full LSI test system was approximately 10 s (16 s if the calculation of transmission line parameters is considered) compared to 4 s for the FDNE representation, however, the speed advantage will increase as the size of the system increases. As the system size increases, the full system will take

372

Y.P. Wang, N.R. Watson / Electrical Power and Energy Systems 44 (2013) 364–374

Full 5−1250 5−2500 5−500 RL

200

Terminal Voltage (kV)

150 100 50 0 −50 −100 −150 −200 0.814

0.816

0.818

0.82

0.822

0.824

0.826

0.828

0.83

Time (sec.) Fig. 14. Comparison of the terminal voltage (Phase A).

Full 5−1250 5−2500 5−500 RL

160

Voltage Magnitude (kV)

140 120 100 80 60 40 20 0

0

5

10

15

20

25

Harmonic Order Fig. 15. Spectrum of the terminal voltage (Phase A).

Table 12 Harmonic voltage level (kV).

Full system FDNE (5–1250 Hz) FDNE (5–2500 Hz) FDNE (5–500 Hz) RL

Table 13 Harmonic voltage maximum error (kV).

11

13

23

25

3.3996 3.8321 3.6929 4.0353 53.9312

3.1579 3.1621 2.6539 2.4272 20.9641

1.5032 1.4770 1.5417 3.3695 17.1984

2.2906 2.3912 2.3355 2.5694 11.0162

longer to simulate while the time required for the FDNE will essentially stay the same. The FDNE fitting software takes less than a 1s and only has to be run once to obtain the rational functions. 4.2. Steady-state tests The next set of tests were performed to look at the steady-state response when a nonlinear (distorting load) is present. A threephase nonlinear load was represented by a harmonic current injection at Tiwai busbar. This source injected the typical harmonic currents of a power electronic converter (given in Table 11). Again the response of the full system is compared with the responses from the FDNEs (using a frequency range of 5–1250 Hz,

Maximum error FDNE (5–1250 Hz) FDNE (5–2500 Hz) FDNE (5–500 Hz) RL

43.253530 50.394044 2.3403355 50.530920

5–2500 Hz and 5–500 Hz) as well as RL representation. The comparison of the phase A terminal voltage time waveform is displayed in Fig. 14. The harmonic components of these waveforms are displayed in Fig. 15. The harmonic levels predicted by the FDNEs match the full system well. Tables 12 and 13 give the assessed harmonic voltage level and error. All the FDNEs give acceptable steady-state performance while as expected the RL equivalent based on 50 Hz does not (see Table 13). 4.3. Converter tests The purpose of a FDNE is to allow test without representing the AC system fully. The next simulation demonstrates this by model-

373

Y.P. Wang, N.R. Watson / Electrical Power and Energy Systems 44 (2013) 364–374

Fig. 16. Expanded view of DC current (fault recovery).

Fig. 17. Comparison of DC current waveform.

ling in detail a converter load connected to Tiwai busbar and simulating how the converter behaves when an AC fault occurs. Figs. 16 and 17 show the relevant comparison. It is clear that any FDNE is better that a representation based on 50 Hz. 5. Conclusions A simple, practical weighted least-squares fitting of a z-domain rational function has been presented and its application to the development of a frequency dependent network equivalent of an actual power system has been demonstrated. Full information has been given on the derived FDNE, allowing the transient response of this FDNE to be easily assessed and compared to other FDNEs developed for the same test system. The Lower South Island of New Zealand test system was chosen as it is a practical system for which all the input parameters are readily available from other publications. The effect of frequency range used to develop the FDNE has been shown and compared to a full representation, as well as an RL equivalent based on fundamental frequency information. Although the frequency range determines the error in the first cycle after a disturbance, after this period the phase shift from a low frequency component dominates. This phase shift can be reduced by weighting the lower frequencies. The importance of fitting at

low frequencies, even those below fundamental frequency has been shown. Fitting down to 1 Hz will greatly remove the low frequency phase-shift that can occur after fault application. For steady-state studies and particularly harmonic assessment, the frequency range was not critical. This work has shown that selection of the upper frequency limit (1250 or 2500 Hz) is far less critical than the weighting of the lower frequencies, particularly for the transient following fault recovery. Acknowledgments This Work was support by the Power Generation & Electricity Delivery of the Korea Institute of Energy Technology Evaluation and Planning (KETEP) grant funded by the Korea Government Ministry of Knowledge Economy (No. 20101020300430). Appendix A. Rational function fitting to frequency response Each term is represented by a rational function of the form.

HðzÞ ¼ ¼

IðzÞ VðzÞ

Pn k a0 þ a1 z1 þ a2 z2 þ    an zn k¼0 ðak  z Þ P 1 þ b1 z1 þ b2 z2 þ    bn zn 1 þ nk¼1 ðbk  zk Þ

ðA1Þ

374

Y.P. Wang, N.R. Watson / Electrical Power and Energy Systems 44 (2013) 364–374

Evaluating the frequency response of the rational function and equating it to the required values (sample data) gives:

HðjxÞ ¼

Pn  kjxDt  k¼0 ak e P 1 þ nk¼1 ðbk ekjxDt Þ

C T ¼ ðW jx1 cðjx1 Þ; W jx2 cðjx2 Þ; . . . ; W jx1N cðjxN ÞÞ DT ¼ ðW jx1 dðjx1 Þ; W jx2 dðjx2 Þ; . . . ; W jx1N dðjxN ÞÞ

ðA2Þ Rik ¼ W jxi ½cðkxi Þ cosðkxi DtÞ þ dðkxi Þ sinðkxi DtÞ

Adding frequency selective weighting yields:

W jx HðjxÞ ¼ W jx a0 

Sik ¼ W jxi ½dðkxi Þ cosðkxi DtÞ  cðkxi Þ sinðkxi DtÞ

n X ½W jx ðbk  HðjxÞ  ak ÞekjxDt 

ðA3Þ

k¼1

Splitting into real and imaginary ekjxDt ¼ cosðkxDtÞ  j sin ðkxDtÞ) gives:

components

(using

n X  W jx cðjxÞ ¼  bk W jx ðcðjxÞ cosðkxDtÞ þ dðjxÞ sinðkxDtÞÞ k¼1

  ak W jx cosðkxDtÞ þ a0 W jx

ðA4Þ

n X W jx dðjxÞ ¼  ½bk W jx ðdðjxÞ cosðkxDtÞ  cðjxÞ

In Eq. (6) C and D are vectors of the negative of the real and imaginary parts of the sample data, respectively, at each frequency. The system of linear equations is derived by evaluating the frequency response of H(z) and equating to the required response at each sample point. This is solved via weighted least-squares. Two equations result from each sample point, one for the real component and the other for the imaginary component. A weighting factor of 100 is applied to equations representing the fundamental frequency so as to ensure minimal steady-state error. The choice of 100 was made on the basis of many simulations and inspecting the steady-state error that occurred.

k¼1

 sinðkxDtÞÞ  ak W jx sinðkxDtÞ

ðA5Þ

The values of a and b coefficients are determined by setting up an over-determined system of linear equations of the form:



A11

A12



a

A21

A22

b

¼

C

ðA6Þ

D

where

2

3

W jx1

W jx1 cosðx1 DtÞ



W jx1 cosðnx1 DtÞ

6 W jx2 6 A11 ¼ 6 6 .. 4 .

W jx2 cosðx2 DtÞ .. .

 .. .

W jx2 cosðnx2 DtÞ 7 7 7; 7 ... 5

W jxN W jxN cosðxN DtÞ    W jxN cosðnxN DtÞ 3 R11 R12    R1n 7 6R 6 21 R22    R2n 7 7 ¼6 . . . . 6 . .. .. .. 7 5 4 . 2

A12

RN1 2

RN2

   RNn 3

0

W jx1 sinðx1 DtÞ



W jx1 sinðnx1 DtÞ

60 6 A21 ¼ 6 6 .. 4.

W jx2 sinðx2 DtÞ .. .

 .. .

W jx2 sinðnx2 DtÞ 7 7 7; .. 7 5 .

0 W jxN sinðxN DtÞ    W jxN sinðnxN DtÞ 3 S11 S12    S1n 7 6 S 6 21 S22    S2n 7 ¼6 .. .. .. 7 7 6 .. 4 . . . . 5 2

A22

SN1 0 SN2

   SNn

aT ¼ ða0 ; a1 ; a2; . . . ; an Þ bT ¼ ðb1 ; a2; . . . ; bn Þ

References [1] Hingorani NG, Burbury M. Simulation of AC system impedance in HVDC system studies. IEEE Trans Power Syst 1970;89(5/6):820–8. [2] Abur A, Singh H. Time domain modeling of external systems for electromagnetic transient programs. IEEE Trans Power Syst 1993;8(2):671–9. [3] Singh H, Abur A. Multi-port equivalencing of external systems for simulation of switching transients. IEEE Trans Power Deliv 1995;10(1):374–82. [4] Noda T, Nagaoka N, Ametani A. Phase domain modeling of frequencydependent transmission lines by means of an ARMA model. IEEE Trans Power Deliv 1996;11(1):401–11. [5] Noda T, Nagaoka N, Ametani A. Further improvements to a phase-domain ARMA line model in terms of convolution, steady-state initialization, and stability. IEEE Trans Power Deliv 1997;12(3):1327–34. [6] Gustavsen B, Semlyen A. Rational approximation of frequency domain responses by vector fitting. In: IEEE/PES winter meeting 1997. Paper no. PE194-PWRD-0-11-1997. [7] Angelidis G, Semlyen A. Direct phase-domain calculation of transmission line transients using two-sided recursions. IEEE Trans Power Deliv 1995;10(2):941–7. [8] Morched AS, Ottevangers JH, Marti L. Multi port frequency dependent network equivalents for the EMTP. IEEE Trans Power Deliv 1993;8(3):1402–12. [9] Watson NR, Gole AM, Irwin GD. Z-domain frequency-dependent network equivalents for electromagnetic transient studies. In: Proceedings of international conference on power system transients (IPST’99); 1999. p. 37–42. [10] Gustavsen B, Semlyen A. Enforcing passivity for admittance matrices approximated by rational functions. IEEE Trans Power Syst 2001;16(1):97–104. [11] Wang YP, Watson NR. Z-domain frequency-dependent A.C. system equivalent for electromagnetic transient simulation. IEE Proc Genr Transm Distrib 2003;150(2):141–6. [12] Watson NR, Arrillaga J. Harmonic assessment using electromagnetic transient simulation and frequency-dependent network equivalents. IEE Proc Genr Transm Distrib 2003;150(6):641–50. [13] Abu-Hashim R, Burch R, Chang G, Grady M, Gunther E, Halpin M, et al. Test systems for harmonics modeling and simulation. IEEE Trans Power Deliv 1999;14(2):579–87. [14] Watson NR, Arrillaga J. Power systems electromagnetic transients simulation. UK: IEE Books; 2002. [15] Noda T. Reduced-order realization of a nonlinear power network using companion-form state equations with periodic coefficients. IEEE Trans Power Deliv 2003;18(2):1478–88. [16] Noda T. Identification of a multiphase network equivalent for electromagnetic transient calculations using partitioned frequency response. IEEE Trans Power Deliv 2005;20(2):1134–42.