**6. Application examples**

The simulation results presented as follows are obtained with the Matlab implementation of the model being described here. These results are compared against those from the phase domain line model in EMTP-RV. Two examples are presented next, first an aerial 9–conductor line and, finally, one for an underground cable. Also a basic m-code for the described phase domain line model is provided at the appendix. The code is given along with the companion routines to perform the first example presented in (Ramos- Leaños & Iracheta, 2010). The reader can readily modify the provided m-code for other applications of interest.

#### **6.1. Aerial line case**

The transversal geometry of this test case is shown in figure 8. Phase conductors are 1192.5 ASCR 54/19 and ground wires are 7 No 5 AWLD. This case consists of three coupled three– phase transmission lines. First line (or circuit 1) is composed of conductors 1 to 3, second line (or circuit 2) includes conductors 5 to 6 and the third line (or circuit 3) comprises conductors 7 to 9. The line length is 150 km. The test circuit is shown in figure 9 where the source is 169 kV, Y-grounded, source impedance is determined by its zero and positive sequence data in Ohms: *R0*=2, *R1*=1, *X0*=22, *X1*=15, and closing times are 0 s for phase a, 0.63 ms for phase b and 0.4 ms for phase c. The simulation time step is 5 s.

**Figure 8.** Aerial line transversal geometry.

**Figure 9.** Test circuit for the case of a nine–conductor line.

Simulation results are presented in figure 10 where the receiving end voltage waveforms of circuit 1 are shown, those for phase a are in blue, those for phase b are in green and those for phase c are in red. A dashed line is used for waveforms obtained with EMTP-RV, while a solid line is used for the results with the line model in Matlab. Notice that the two sets of results overlap and not difference can be seen. Figure 10 provides the differences between the two sets of results. Note that the largest difference is around 3e-9.

**Figure 10.** (a) Over voltages at receiving end for conductors 1, 2 and 3, (b) Differences between results with Matlab model and with EMTP–RV.

#### **6.2. Underground cable case**

The underground cable system used for this test consists of three single–phase coaxial cables, its transversal layout is shown in figure 11. The Corresponding connection diagram is provided in figure 12. Circuit parameters are given in table 1, the cable length is 6.67km and the time step used for the simulation is 1 s. The applied excitation is by a 3ph 169kV ideal source.

The simulation experiment consists in the simultaneous energizing of the three cable cores. The results presented in figure 13 correspond to the core voltages at the far end. Phase a voltages are in blue, phase b voltages are in green and those for phase c are in red. A dashed line is used for the results obtained with EMTP-RV, while a solid line is used for the Matlab model results. Notice that both sets of results overlap and that no difference can be seen by eye. Figure 13 also shows the difference between the two sets of results which is around 4e-9. Compared to the 1.69e+5 amplitude of the excitation source, this difference shows the outstanding accuracy of the Matlab model.

**Figure 11.** Cable layout.

284 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 1

+ SW1

Simulation results are presented in figure 10 where the receiving end voltage waveforms of circuit 1 are shown, those for phase a are in blue, those for phase b are in green and those for phase c are in red. A dashed line is used for waveforms obtained with EMTP-RV, while a solid line is used for the results with the line model in Matlab. Notice that the two sets of results overlap and not difference can be seen. Figure 10 provides the differences between

a b c

WB\_k

+ WB

c b a WB\_m2

WB\_m3

a b c WB\_m1

a b c

**0 0.005 0.01 0.015 0.02**

**Time (s)**

**Figure 10.** (a) Over voltages at receiving end for conductors 1, 2 and 3, (b) Differences between results

(a) (b)

**Absolute error**

**-4 -3 -2 -1 0 1 2 3 x 10-9**

The underground cable system used for this test consists of three single–phase coaxial cables, its transversal layout is shown in figure 11. The Corresponding connection diagram is provided in figure 12. Circuit parameters are given in table 1, the cable length is 6.67km and the time step

The simulation experiment consists in the simultaneous energizing of the three cable cores. The results presented in figure 13 correspond to the core voltages at the far end. Phase a voltages are in blue, phase b voltages are in green and those for phase c are in red. A dashed line is used for the results obtained with EMTP-RV, while a solid line is used for the Matlab

used for the simulation is 1 s. The applied excitation is by a 3ph 169kV ideal source.

**Figure 9.** Test circuit for the case of a nine–conductor line.

<sup>+</sup> RL1

+ 1k R1

+

169kV /\_0

AC1

with Matlab model and with EMTP–RV.

**0 0.005 0.01 0.015 0.02**

**Time (s)**

**6.2. Underground cable case** 

**-3 -2 -1 0 1 2 3 x 105**

**Voltage (V)**

the two sets of results. Note that the largest difference is around 3e-9.


**Table 1.** Cable data.

**Figure 12.** Cable test circuit.

**Figure 13.** (a) Receivng end core voltages, (b) absolute error.

#### **7. Vector fitting**

The goal of VF is to approximate a complex function of frequency by means of a rational function; that is, a quotient of two polynomials of the frequency variable (Gustavsen & Semlyen, 1999). The function to be approximated could be trascendantal or could be specified by its values at a number of frequency points. The form of the approximation obtained with VF is that of a partial fraction expansion:

$$f(\mathbf{s}) \equiv \sum\_{n=1}^{N} \frac{r\_n}{\mathbf{s} + \overline{p}\_n} \tag{55}$$

VF estimates the system parameters by means of a two-stage linear least–squares procedure. First a set of initial poles for the partial fraction basis (55) is selected and relocated iteratively until a prescribed convergence criterion is attained. Then, convergence is tested by means of a second linear least–squares procedure in which the previously obtained poles are fixed and the corresponding residues are taken as the unknown parameters.

Consider the following relation (Gustavsen & Semlyen, 1999):

$$\sum\_{n=1}^{N} \frac{\hat{r}\_n}{s + \overline{p}\_n} \equiv f(s) \left( 1 + \sum\_{n=1}^{N} \frac{\tilde{r}\_n}{s + \overline{p}\_n} \right) \tag{56}$$

where, *N* is the order of approximation, *<sup>n</sup> p* represents the unknown poles and ˆ *n r* and *<sup>n</sup> r* are unknown residues. Poles are initialized by values distributed logarithmically over the frequency range of interest. Expression (56) is now rewritten as follows:

$$\sum\_{n=1}^{N} \frac{\hat{r}\_n}{s + \overline{p}\_n} - \left(\sum\_{n=1}^{N} \frac{\tilde{r}\_n}{s + \overline{p}\_n}\right) f(\mathbf{s}) \equiv f(\mathbf{s}).\tag{57}$$

An over–determined least squares equation–system is then obtained by evaluating (57) at a number *M* of specific frequencies, with *M>2N*:

$$\mathbf{Ax} = \mathbf{b}\_{\prime} \tag{58}$$

where *A* is the *M2N* matrix whose elements depend on the poles, *x* is the *2N*–dimension vector of unknown residues and *b* is the *M*–dimension vector with the values of the function to be approximated (Gustavsen & Semlyen, 1999). Special care is taken to accommodate next to each other those complex–conjugate pairs of pole–residues that can arise. Expression (58) is solved through an iterative process represented symbolically as follows:

$$\mathbf{A}^{(j-1)}\mathbf{x}^{(j)} = \mathbf{b},\tag{59}$$

were (j–1) and (j) represent super–indexes and j is the iteration index. *A*(0) is obtained from the initial poles with logarithmic distribution over the frequency range of interest (Gustavsen & Semlyen, 1999). As (59) is solved in the first iteration, a second step is to use the obtained residue values to recalculate new poles for the function to be fitted *f(s)*. This is accomplished by computing the eigenvalues of the following matrix *Q* (Gustavsen & Semlyen, 1999):

286 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 1

obtained with VF is that of a partial fraction expansion:

Consider the following relation (Gustavsen & Semlyen, 1999):

number *M* of specific frequencies, with *M>2N*:

The goal of VF is to approximate a complex function of frequency by means of a rational function; that is, a quotient of two polynomials of the frequency variable (Gustavsen & Semlyen, 1999). The function to be approximated could be trascendantal or could be specified by its values at a number of frequency points. The form of the approximation

1

VF estimates the system parameters by means of a two-stage linear least–squares procedure. First a set of initial poles for the partial fraction basis (55) is selected and relocated iteratively until a prescribed convergence criterion is attained. Then, convergence is tested by means of a second linear least–squares procedure in which the previously obtained poles are fixed and the corresponding residues are taken as the

*<sup>r</sup> f s*

1 1

where, *N* is the order of approximation, *<sup>n</sup> p* represents the unknown poles and ˆ

frequency range of interest. Expression (56) is now rewritten as follows:

1 1

is solved through an iterative process represented symbolically as follows:

*N N n n n n n n*

*sp sp*

*N N n n n n n n r r f s s p s p*

unknown residues. Poles are initialized by values distributed logarithmically over the

*r r <sup>f</sup> s fs*

 

An over–determined least squares equation–system is then obtained by evaluating (57) at a

vector of unknown residues and *b* is the *M*–dimension vector with the values of the function to be approximated (Gustavsen & Semlyen, 1999). Special care is taken to accommodate next to each other those complex–conjugate pairs of pole–residues that can arise. Expression (58)

were (j–1) and (j) represent super–indexes and j is the iteration index. *A*(0) is obtained from the initial poles with logarithmic distribution over the frequency range of interest

()1 ,

 

( ) ( ).

*2N* matrix whose elements depend on the poles, *x* is the *2N*–dimension

*<sup>N</sup> <sup>n</sup> n n*

*s p*

(55)

(56)

(57)

**Ax b** , (58)

( 1) ( ) , *j j* **Ax b** (59)

*n r* and *<sup>n</sup>*

*r* are

( )

**7. Vector fitting** 

unknown parameters.

where *A* is the *M*

$$\mathbf{Q} = \mathbf{W} - \mathbf{g}\mathbf{\tilde{x}}^T,\tag{60}$$

where *W* is a diagonal matrix containing previously calculated poles *<sup>n</sup> p* , *g* is a vector of ones and **x** is a vector containing the *r* terms only. The reason for using (60) is explained next. Let (56) be rewritten as follows:

$$\frac{\sum\_{n=1}^{N} \frac{\tilde{r}\_n}{s + \overline{p}\_n}}{\sum\_{n=1}^{N} \frac{\tilde{r}\_n}{s + \overline{p}\_n} + 1} = \frac{\prod\_{n=1}^{N} (s + \hat{z}\_n)}{\prod\_{n=1}^{N} (s + \tilde{z}\_n)} \equiv f(s). \tag{61}$$

It is clear in (61) that the two polynomials containing the poles *<sup>n</sup> p* cancel each other, and that the zeros *nz* become the poles of *f(s)*. Notice further that the denominator on the left– hand–side of (61) can be written as follows:

$$\sum\_{n=1}^{N} \frac{\tilde{r}\_n}{s + \overline{p}\_n} + 1 = \frac{\prod\_{n=1}^{N} (s + \tilde{z}\_n)}{\prod\_{n=1}^{N} (s + \overline{p}\_n)}.\tag{62}$$

Zeros *nz* are then obtained by finding the roots of (Gustavsen & Semlyen, 1999)

$$\sum\_{i=1}^{N} \left( \tilde{r}\_i \prod\_{n=1, n \neq i}^{N} (\mathbf{s} + \overline{p}\_n) \right) + \prod\_{n=1}^{N} (\mathbf{s} + \overline{p}\_n) = \mathbf{0},\tag{63}$$

which is equivalent to finding the eigenvalues of *Q* in (60) (Gustavsen & Semlyen, 1999).

The newly found set of poles is replaced in (55) to determine a new set of residues *rn*. This is again an over–determined linear system. The fitting error is tested at this stage for each available sample of *f(s)*. If the error level is not acceptable, the new poles are used to restart the procedure as with (56). If the desired error limit is not met after a pre-specified number of iterations, then, the order of approximation *N* is increased and the iterative procedure is restarted (Gustavsen & Semlyen, 1999).

Even in most cases where initial poles are not chosen adequately, VF is capable of finding a solution at the expense of more iterations. In some cases an iteration can produce unstable poles; these poles simply are flipped into the left–hand–side part of the complex plane (i.e., the stable part) and a new solution is searched (Gustavsen & Semlyen, 1999).
