**4. Phase domain line model**

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

**Figure 2.** Frequency domain circuit representation of a multi-conductor line.

instance, the transformation of (21) to the time domain yields:

represented as follows:

auxiliary sources *Iaux,0* and *Iaux,L*.

currents can be represented as

with

and

with, *Ish,0* =*YCV0* , *Iaux,0* =*HIrfl,L* ,and *Irfl,L*= *IL+ YCVL* .

In the same way as it has been previously done for (19), expression (20) is conveniently

0 ,0 ,0 *sh aux* **II I** (22)

Expressions (21) and (22) constitute a traveling wave line model for the segment of length *L* depicted in figure 1. The former set of expressions represents end *L* of the line segment, while the latter set represents end *0*. A schematic representation for the whole model is provided by figure 2. Note that the coupling between the two line ends is through the

The line model defined by expressions (21) and (22) is in the frequency domain. Power system transient simulations require this model to be transformed to the time domain. For

Note that at (23), (24), (25) the lowercase variables represent the time domain images of their uppercase counterparts at (22) and that the symbol \* represents convolution. Reflected

Expressions (23)-(26) constitute a general traveling–wave based time–domain model for line end *0*. The model corresponding to the other end is obtained by interchanging sub–indexes "*0*" and "*L*" at (23)-(26). Equation (23) essentially provides the interface of the line–end *0* model to the nodal network solver that, for power system transient analysis, usually is the EMTP (Dommel, 1996). Expressions (24) and (25) require the performing of matrix–to–vector

0 ,0 ,0 *sh aux* **ii i** (23)

,0 \* *sh C L* **i yv** (24)

,0 ,0 \* *aux rfl* **i hi** (25)

, , , 2 *rfl L sh L aux L* **i ii** (26)

Since rational fitting and model solutions are carried out directly in the phase domain, the model described here is said to be a phase domain line model. Rational fitting for this model is carried out using the Vector fitting (VF) tool (Gustavsen, 2008). In the case of *YC,* the whole fitting process is done in the phase domain, whereas for *H* initial poles and time delays are first calculated in the modal domain.

#### **4.1. Rational approximation of Yc**

The following rational representation has been proposed for *YC* in (Marti, 1982) and (Morched et al., 1999):

$$\mathbf{Y}\_{\mathcal{C}} = \mathbf{G}\_0 + \sum\_{i=1}^{Ny} \frac{\mathbf{G}\_i}{\mathbf{s} - \boldsymbol{\eta}\_i} \tag{27}$$

where *Ny* is the fitting order, *qi* represents the *i–th* fitting pole, *Gi* is the corresponding matrix of residues and *G0* is a constant matrix obtained at the limit of *YC* when *s=j*. Note in (27) that common poles are used for the fitting of all the elements of *YC* obtained by fitting the matrix trace and finally the fitting of the matrices of residues *Gi* and of proportional terms *G0* is done in the phase domain. Section 7 provides an overview of the VF procedure and further information is to be found in (Gustavsen & Semlyen, 1999) and (Gustavsen, 2008).

#### **4.2. Rational approximation of H**

Rational fitting of transfer matrix *H* is substantially more involving than the one of *YC* above. To attain an accurate and compact (low order) rational representation for *H* it is essential to factor out all terms involving time delays (Marti, 1982). The major difficulty here is that its elements could involve a mix of up to *N* different delay terms due to the multimode propagation on an *N*–conductor line (Wedepohl, 1965). Separation of matrix *H* into singledelay terms is obtained from the following modal factorization (Wedepohl, 1965):

$$\mathbf{H} = \mathbf{M} \mathbf{H}\_m \mathbf{M}^{-1} \tag{28}$$

where *Hm* is a diagonal matrix of the form

$$\mathbf{H}\_m = \text{diag}\{e^{-\gamma\_1 L}, e^{-\gamma\_2 L}, \dots, e^{-\gamma\_N L}\} \tag{29}$$

and =(*YZ*) is the propagation constant of a conductor line (Wedepohl, 1965). As the triple product in (28) is performed by partitioning *M* in columns and *M–1* in rows, the following expression is derived:

$$\mathbf{H} = \sum\_{i=1}^{N} \mathbf{D}\_i e^{-\boldsymbol{\gamma}\_i L} \tag{30}$$

where *Di* is the rank–1 matrix obtained from pre–multiplying the *i–th* column of *M* by the *i– th* row of *M–1*. Matrix *Di* in fact is an idempotent (Marcano & Marti, 1997). The exponential factors at (30) can be further decomposed as follows:

$$e^{-\gamma\_i L} = e^{-\gamma\_i L} e^{-s\tau\_i}; i = 1, 2, \dots, N \tag{31}$$

where exp *<sup>i</sup> L* is a minimum phase–shift function (Bode, 1945) and *<sup>i</sup>* is the time delay associated to the velocity of the *i–th* mode. Hence:

$$\mathbf{H} = \sum\_{i=1}^{N} \mathbf{D}\_i e^{-\mathcal{F}\_i L} e^{-s\mathbf{r}\_i} \tag{32}$$

The time delays in (32) can be initially estimated by applying Bode's relation for minimum phase complex functions (Bode, 1945) to the magnitudes of *exp(-iL)* in (30). Although (32) provides the desired separation of *H* as a sum of terms, each one involving a single delay factor, the following consideration is brought in for computational efficiency (Morched et al., 1999). Modal delays often occur in groups with almost identical values. Suppose that a number *Ng* of these groups can be formed and that (32) can be modified as follows:

$$\mathbf{H} = \sum\_{k=1}^{N\_{\mathcal{S}}} \tilde{\mathbf{H}}\_k e^{-s\tau\_k} \tag{33}$$

where *Ng* is less or equal to *N*, and *<sup>k</sup>* is the representative delay for the *k–th* group. Clearly, by comparing (33) and (32):

$$\tilde{\mathbf{H}}\_k = \sum\_{i=1}^{lk} \mathbf{D}\_i e^{-\tilde{\gamma}\_i L} e^{-s\tau\_i} k = 1, 2, \dots, \text{Ng} \tag{34}$$

with *Ik* being the number of modal terms in the *k–th* group. To determine whether a set of exponential factors can be grouped or not, the maximum phase shifts associated to their time delays are compared. The set is grouped into a single delay group if the phase shift differences are below a pre-established value typically chosen at 10 (Morched et al., 1999). Each term **<sup>H</sup>***<sup>k</sup>* at (34) can now be considered free of delay factors and can be fitted as follows:

$$\tilde{\mathbf{H}}\_k = \sum\_{i=1}^{Nh(k)} \frac{\mathbf{R}\_{k,i}}{s - p\_{k,i}} k = 1, 2, \dots, Ng \tag{35}$$

where *Nh(k)* is the fitting order for the *k–th* term **<sup>H</sup>***<sup>k</sup>* , *pk,i* represents its *i–th* fitting pole and *Rk,i* is the corresponding matrix of residues. Note in (35) that common poles are being used to fit all elements at each matrix **<sup>H</sup>***<sup>k</sup>* . As (35) is introduced in (33), the following rational form is obtained for *H*:

$$\mathbf{H} = \sum\_{k=1}^{N\_{\mathcal{S}}} e^{-s\tau\_k} \sum\_{i=1}^{Nh(k)} \frac{\mathbf{R}\_{k,i}}{s - p\_{k,i}} \tag{36}$$

Initial estimates for the poles as well as time delays are obtained in the modal domain. The poles result from applying VF to the sum of the modal exponential factors conforming each delay group. The time delays proceed from Bode's formula as it has been said before. With all the poles *pk,i* and group delays *<sup>k</sup>* of (36) being estimated in the modal domain, the overall fitting of *H* is completed in the phase domain by obtaining the matrices of residues *Rk,i* and recalculating the poles (Gustavsen & Nordstrom, 2008). The fitting parameters so obtained can be taken as final or can be further refined by an iterative process. VF is applied throughout all these fitting tasks and detailed descriptions of these processes can be found in (Gustavsen & Nordstrom, 2008; Gustavsen & Semlyen, 1999).

#### **4.3. State-space analysis**

With the rational representation of *YC* and *H* state–space forms to evaluate *ish,0* and *iaux,0* arise naturally. Consider first the case of *ish,0*. Taking (22) and introducing (27) yields

$$\mathbf{I}\_{sh,0} = \mathbf{G}\_0 \mathbf{V}\_0 + \sum\_{i=1}^{Ny} \mathbf{W}\_i \tag{37}$$

where

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

 

 

*L* is a minimum phase–shift function (Bode, 1945) and

1

The time delays in (32) can be initially estimated by applying Bode's relation for minimum

provides the desired separation of *H* as a sum of terms, each one involving a single delay factor, the following consideration is brought in for computational efficiency (Morched et al., 1999). Modal delays often occur in groups with almost identical values. Suppose that a

1

*k*

*Ik L s*

( ) , 1 ,

*i k i*

*e*

 

 

where *Nh(k)* is the fitting order for the *k–th* term **<sup>H</sup>***<sup>k</sup>* , *pk,i* represents its *i–th* fitting pole and *Rk,i* is the corresponding matrix of residues. Note in (35) that common poles are being used to fit all elements at each matrix **<sup>H</sup>***<sup>k</sup>* . As (35) is introduced in (33), the following rational

 *s p* **R**

*Nh k k i*

 

with *Ik* being the number of modal terms in the *k–th* group. To determine whether a set of exponential factors can be grouped or not, the maximum phase shifts associated to their time delays are compared. The set is grouped into a single delay group if the phase shift differences are below a pre-established value typically chosen at 10 (Morched et al., 1999). Each term **<sup>H</sup>***<sup>k</sup>* at (34) can now be considered free of delay factors and can be fitted as follows:

*Ng <sup>s</sup> k*

*e* 

1,2,..., *i i*

*e e k Ng*

1,2,...,

*k Ng*

( ) ,

*s p*

**R**

1 1 , *k Ng Nh k <sup>s</sup> k i k i k i*

number *Ng* of these groups can be formed and that (32) can be modified as follows:

1

*k i i*

*k*

phase complex functions (Bode, 1945) to the magnitudes of *exp(-*

factors at (30) can be further decomposed as follows:

associated to the velocity of the *i–th* mode. Hence:

where exp *<sup>i</sup>*

where *Ng* is less or equal to *N*, and

by comparing (33) and (32):

form is obtained for *H*:

1

where *Di* is the rank–1 matrix obtained from pre–multiplying the *i–th* column of *M* by the *i– th* row of *M–1*. Matrix *Di* in fact is an idempotent (Marcano & Marti, 1997). The exponential

; 1,2, , *i ii L Ls e ee i N*

*i <sup>N</sup> <sup>L</sup> i i*

*i i <sup>N</sup> L s i i*

*k*

*e e* 

**H D** (30)

(31)

**H D** (32)

**H H** (33)

*<sup>k</sup>* is the representative delay for the *k–th* group. Clearly,

**H D** (34)

**H**(35)

**H** (36)

*<sup>i</sup>* is the time delay

*iL)* in (30). Although (32)

*e* 

$$\mathbf{W}\_{i} = \frac{\mathbf{G}\_{i}}{\mathbf{s} - \mathbf{q}\_{i}} \mathbf{V}\_{0} \text{ } i = \mathbf{1} \text{ } \mathbf{2} \text{ } \dots \text{N} \text{y} \tag{38}$$

Application of the Inverse Laplace Transform to (37) and (38) gives the following continuous-time-state-space (CTSS) form for *ish,0*:

$$\mathbf{i}\_{sh,0} = \mathbf{G}\_0 \mathbf{v}\_0 + \sum\_{i=1}^{Ny} \mathbf{w}\_i \tag{39}$$

and

$$\frac{d\mathbf{w}\_i}{dt} = q\_i \mathbf{w}\_i + \mathbf{G}\_i \mathbf{v}\_0 \text{ \(i = 1, 2, \dots, Ny\)}\tag{40}$$

The CTSS form to evaluate *iaux,0*, is derived now. On replacing the fitted form of *H* given by (36) into (22):

$$\mathbf{I}\_{\text{aux},0} = \sum\_{k=1}^{\text{Ng}} \sum\_{i=1}^{\text{Nh}(k)} \mathbf{X}\_{k,i} \tag{41}$$

with

$$\mathbf{X}\_{k,i} = \frac{R\_{k,i}}{s - p\_{k,i}} I\_{r\emptyset,\iota} e^{-s\tau\_k} ; \quad \begin{array}{l} i = 1, 2, \ldots, N\text{y} \\ k = 1, 2, \ldots, N\text{g} \end{array} \tag{42}$$

Application of the Inverse Laplace Transform to (41) and (42) renders the following CTSS form:

$$\mathbf{i}\_{aux,0} = \sum\_{k=1}^{\text{Ng}} \sum\_{i=1}^{\text{Nh}(k)} \mathbf{x}\_{k,i} \tag{43}$$

$$\frac{d\mathbf{x}\_{k,i}}{dt} = p\_{k,i}\mathbf{x}\_{k,i} + \mathbf{R}\_{k,i}\mathbf{i}\_{rfl,L}(t-\tau\_k);\ \begin{array}{c} i = 1,2,...,Ny\\k = 1,2,...,Ny \end{array} \tag{44}$$

CTSS forms (39), (40), (43) and (44) provide the basis for a phase domain line model (Morched et al., 1999). Nevertheless, their solution by a digital processor requires the conversion to discrete–time state–space (DTSS). This is accomplished by applying a numerical differentiation rule to the CTSS forms. The one adopted here is the mid–point rule of differentiation, which is equivalent to the trapezoidal integration rule extensively used in EMTP (Dommel, 1969, 1992). Application of this rule to (44) with *t* as the solution time step results in:

$$\mathbf{x}\_{k,i} = a\_{k,i}\mathbf{x}\_{k,i}^{\prime} + \tilde{\mathbf{R}}\_{k,i}[\mathbf{i}\_{r\emptyset,L}(t-\tau\_k) + \mathbf{i}\_{r\emptyset,L}^{\prime}(t-\tau\_k)]; \quad \begin{array}{l} i = 1,2,...,Ng \\ k = 1,2,...,Ng \end{array} \tag{45}$$

where *ak,i*=(*2+tpk,i*)/ (*2-tpk,i*) and **R***k i*, =(*tRk,i*)/ (*2-tpk,i*). *xk,i* are discrete-state variables and primed variables denote their value at one previous time step *x*'*k,i= xk,i*(*t-t*). The discrete– time version of (43) maintains its original form:

$$\mathbf{i}\_{\text{aux},0} = \sum\_{k=1}^{N\_{\text{S}}} \sum\_{i=1}^{Nh(k)} \mathbf{x}\_{k,i}$$

Transmission line simulation of EMTs requires the use of time steps *t* smaller than any of the travel times *<sup>k</sup>* in the line. Hence, (45) provides the update of state vectors *xk,i* using only past values of variables already available, either from initial conditions or from previous simulation time steps.

The differentiation mid–point rule is now applied to (40):

$$\mathbf{w}\_{i} = a\_{i}\mathbf{w}\_{\ i}^{\prime} + \tilde{\mathbf{G}}\_{i}(\upsilon\_{0} + \upsilon\_{0}^{\prime}); i = 1, 2, \dots, N y \tag{46}$$

where *ai*=(*2+tqi*)/ (*2-tqi*) and **G***<sup>i</sup>* =(*tGi*)/ (*2-tqi*)

Expression (46) is not a proper DTSS form, as *wi* depends on the present–time value of *v0* which still is to be determined (Gustavsen & Mahseredjian, 2007). This problem is fixed here with the following redefinition of the state variable vector:

$$\mathbf{y}\_{i} = (\tilde{\mathbf{G}}\_{i}^{-1} \mathbf{w}\_{i} - \upsilon\_{0}) / (a\_{i} + 1); i = 1, 2, \dots, Ny \tag{47}$$

Introducing (47) in (46) and (39) the following expressions are obtained:

$$\mathbf{y}\_{i} = a\_{i}\mathbf{y}\_{i}^{\prime} + \boldsymbol{\upsilon}\_{\;\;0}^{\prime}; i = 1, 2, \ldots, N \\ y \tag{48}$$

$$\mathbf{\dot{i}}\_{sh,0} = \mathbf{G} \mathbf{v}\_0 + \mathbf{\dot{i}}\_{y-aux,0} \tag{49}$$

where

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

,, , ,

, ,, , , ,

primed variables denote their value at one previous time step *x*'*k,i= xk,i*(*t-*

Transmission line simulation of EMTs requires the use of time steps

The differentiation mid–point rule is now applied to (40):

*tqi*) and **G***<sup>i</sup>*

with the following redefinition of the state variable vector:

 =(*tGi*)/ (*2-*

1

Introducing (47) in (46) and (39) the following expressions are obtained:

*tpk,i*) and **R***k i*, =(

*k i k i k i rfl L k <sup>d</sup> i Ny p t dt k Ng*

CTSS forms (39), (40), (43) and (44) provide the basis for a phase domain line model (Morched et al., 1999). Nevertheless, their solution by a digital processor requires the conversion to discrete–time state–space (DTSS). This is accomplished by applying a numerical differentiation rule to the CTSS forms. The one adopted here is the mid–point rule of differentiation, which is equivalent to the trapezoidal integration rule extensively used in EMTP (Dommel, 1969, 1992).

> 1,2,..., ' [ ( ) ' ( )]; 1,2,..., *k i k i k i k i rfl L k rfl L k i Ny a tt k Ng*

> > *tRk,i*)/ (*2-*

past values of variables already available, either from initial conditions or from previous

0 0 ' ( ' ); 1,2,..., *i ii i* **w wG** *a v v i Ny* (46)

*tqi*)

Expression (46) is not a proper DTSS form, as *wi* depends on the present–time value of *v0* which still is to be determined (Gustavsen & Mahseredjian, 2007). This problem is fixed here

,

*k i*

**x**

Application of this rule to (44) with

*tpk,i*)/ (*2-*

time version of (43) maintains its original form:

where *ak,i*=(*2+*

the travel times

where *ai*=(*2+*

simulation time steps.

*tqi*)/ (*2-*

( ) ,0 , 1 1 *Ng Nh k aux k i k i*

1,2,..., ( ); 1,2,...,

*t* as the solution time step results in:

( ) ,0 , 1 1 *Ng Nh k aux k i k i* **i x**

*<sup>k</sup>* in the line. Hence, (45) provides the update of state vectors *xk,i* using only

<sup>0</sup> ( ) / ( 1); 1,2,..., *i ii i v a i Ny* **y Gw** (47)

<sup>0</sup> ' ' ; 1,2,..., *i ii* **y y** *a v i Ny* (48)

*sh*,0 0 ,0 *y aux* **i Gv i** (49)

 **x x Ri i** (45)

**i x** (43)

**x Ri** (44)

*tpk,i*). *xk,i* are discrete-state variables and

*t*). The discrete–

*t* smaller than any of

$$\mathbf{i}\_{y-aux,0} = \sum\_{i=1}^{Ny} \hat{\mathbf{G}}\_i \mathbf{y}\_i \; ; \; \hat{\mathbf{G}}\_i = (a\_i + 1)\tilde{\mathbf{G}}\_i \; ; \; \mathbf{G} = \mathbf{G}\_0 + \sum\_{i=1}^{Ny} \tilde{\mathbf{G}}\_i \mathbf{e}\_i$$

Expression (48) is a proper DTSS form for the sequential evaluation of *ish,0* at the phase domain line model.

Finally, the introduction of (43) and (49) in (23) results in:

$$
\dot{\mathbf{i}}\_0 = \mathbf{G} \mathbf{v}\_0 + \dot{\mathbf{i}}\_{hist, 0} \tag{50}
$$

with

$$\mathbf{i}\_{hist,0} = \mathbf{i}\_{y-aux,0} - \mathbf{i}\_{aux,0} = \sum\_{i=1}^{Ny} \hat{\mathbf{G}}\_i \mathbf{y}\_i - \sum\_{k=1}^{Ng} \sum\_{i=1}^{Nh(k)} \mathbf{x}\_{k,i}$$

Expression (50), along with (48) and (49), provides a discrete time–domain model for end *0* of the line segment at figure 1. The expressions for the model at end *L* are simply obtained by exchanging sub–indexes *0* and *L* at (48), (49) and (50). Obviously, state variables "*yi*" and "*xk,i*" of end *L* model are different from those of end *0*. Figure 3 provides a discrete–time circuit–model for the line segment of length *x=L*. This model is based on expression (50) and its companion for line end *L*. Note that the model consists of parallel arrangements of shunt conductances and auxiliary sources of currents comprising historic terms of ends (or nodes) *0* and *L*. Figure 3 model is thus in an appropriate form for computer code implementation. In this chapter, the Matlab environment has been chosen for this end.

**Figure 3.** Discrete time domain circuit representation of a multi-conductor line.

#### **5. Line model implementation in Matlab**

The discrete–time line model depicted in figure 3 and defined by (50) has been programmed by these authors in Matlab as an M–code function (see Appendix). This function consists of two sub–blocks, one for each multi-conductor line end. This model is to be used with a nodal network solver, a complete explanation on the nodal solver can be found in (Dommel, 1969 & 1992). Expression (50) constitutes essentially the interface between the line model and the nodal solver. Each one of the two sub–blocks in the line model performs iteratively the six tasks that are described next for line–end *0* sub–block. Figure 4 provides the block diagram of the complete line/cable model, along with its interfacing with the nodal solver.


$$\mathbf{i}\_{sh,0} = \mathbf{G}\mathbf{v}\_0 + \mathbf{i}\_{y-au\alpha,0}$$

**Step 3.** Auxiliary current source value, due to the reflected traveling waves at the remote line end, is updated by (43):

$$\mathbf{i}\_{\text{aux},0} = \sum\_{k=1}^{N\_{\text{S}}} \sum\_{i=1}^{Nh(k)} \mathbf{x}\_{k,i}$$

**Step 4.** Vector of reflected currents at the local line–end (node) "*irfl,0*" is calculated for the present time by means of (26) being modified to suit line–end *0*:

$$\mathbf{i}\_{r\emptyset,0} = \mathbf{2i}\_{sh,0} - \mathbf{i}\_{au\alpha,0}$$

This vector is delivered to end *L* sub–block through a delay buffer. Although branch current vector *i0* usually is not explicitly required, it is conveniently evaluated here by (50):

$$\mathbf{i}\_0 = \mathbf{G}\mathbf{v}\_0 + \mathbf{i}\_{hist, 0}$$

**Step 5.** Internal states inside the line model are updated by (48) and (45):

$$\mathbf{y}\_i = a\_i \mathbf{y}'\_i + \boldsymbol{\upsilon}'\_0$$

$$\mathbf{x}\_{k,i} = a\_{k,i} \mathbf{x}'\_{k,i} + \tilde{\mathbf{R}}\_{k,i} [\mathbf{i}\_{r\emptyset,L}(t - \boldsymbol{\tau}\_k) + \mathbf{i}'\_{r\emptyset,L}(t - \boldsymbol{\tau}\_k)]$$

**Step 6.** The vector of history currents for end (node) *0* is updated by means of (50) and the update is delivered to the nodal-network solver.

Steps *1* to *6* are iterated *Nt* times until *Ntt* spans the total simulation time of interest.

#### **5.1. Handling of line-travel delays**

It follows from expressions (43) and (45) that the calculation of *iaux,0* requires the reflected currents vector *irfl,L* being evaluated with various time delays , …, *Ng*. Recall that the delays are due to the travel time needed by a wave to travel from one line end to the other. Past values of *irfl,L* can be obtained either from line initial conditions or from previous simulation steps; nevertheless, these values are given by samples regularly distributed *t* seconds apart. Since the involved travel times (or line delays) usually are not integer multiples of *t*, the required values of *irfl,L* must be obtained by means of interpolations. The standard procedure for this is to interpolate linearly (Dommel, 1992) and this is adopted here.

Evaluation of the delayed values require a memory buffer spanning at least the largest travel time

$$
\pi\_{\text{max}} = \max \{ \tau\_{1'} \tau\_{2'} ... \tau\_{\text{Ng}} \} \, \, \, \tag{51}
$$

and buffer length *Nb* is calculated as follows:

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

to determine line end (nodal) voltages *v0* and *vL*.

repeated here for convenience:

line end, is updated by (43):

**Step 1.** State–variable and history–current values are assumed known, either from initial conditions or from previous simulation steps. These values are used by the nodal solver

**Step 2.** Shunt current due to the characteristic admittance of the line is calculated by (49)

*sh*,0 0 ,0 *y aux* **i Gv i**

**Step 3.** Auxiliary current source value, due to the reflected traveling waves at the remote

**Step 4.** Vector of reflected currents at the local line–end (node) "*irfl,0*" is calculated for the

,0 ,0 ,0 2 *rfl sh aux* **i ii**

This vector is delivered to end *L* sub–block through a delay buffer. Although branch current

0 0 ,0 *hist* **i Gv i**

<sup>0</sup> ' ' *i ii* **y y** *a v*

, ,, , , , ' [ ( ) ' ( )] *k i k i k i k i rfl L k rfl L k* **x x Ri i** *a tt*

**Step 6.** The vector of history currents for end (node) *0* is updated by means of (50) and the

It follows from expressions (43) and (45) that the calculation of *iaux,0* requires the reflected

delays are due to the travel time needed by a wave to travel from one line end to the other. Past values of *irfl,L* can be obtained either from line initial conditions or from previous simulation steps; nevertheless, these values are given by samples regularly distributed

seconds apart. Since the involved travel times (or line delays) usually are not integer

standard procedure for this is to interpolate linearly (Dommel, 1992) and this is adopted

*t*, the required values of *irfl,L* must be obtained by means of interpolations. The

*t* spans the total simulation time of interest.

, …, 

*Ng*. Recall that the

*t*

vector *i0* usually is not explicitly required, it is conveniently evaluated here by (50):

present time by means of (26) being modified to suit line–end *0*:

**Step 5.** Internal states inside the line model are updated by (48) and (45):

update is delivered to the nodal-network solver.

currents vector *irfl,L* being evaluated with various time delays

Steps *1* to *6* are iterated *Nt* times until *Nt*

**5.1. Handling of line-travel delays** 

multiples of

here.

( ) ,0 , 1 1 *Ng Nh k aux k i k i* **i x**

$$N\_b = \left\lfloor \frac{\tau\_{\text{max}}}{\Delta t} \right\rfloor + 1\tag{52}$$

**Figure 4.** Line/Cable model's complete flow diagram.

If a propagation delay is an integer multiple of *t*, the required value of *irfl* can be readily retrieved from the memory buffer. This is illustrated by figure 5 where the simulation time step is *t=0.03* ms and the travel time is =0.10 ms. It can be seen that at simulation time *t= 0.24* ms the required history value at 0.09ms is available from the table.

On a multiphase system, nevertheless, it is highly improbable that all the propagation times can be made integer multiples of a single value of *t* suitable for transient simulations. Thus, the required values must be obtained through interpolation. Figure 6 illustrates this case, where a simulation time step *t=0.04* ms is assumed instead of the *t=0.03* ms one at figure 5. Notice that now the required history value, for a time delay of 0.09ms, is not readily available.

Suppose now that the required value *irfl(t–)* is between the *k–th* and the *(k+1)–th* stored samples of *irfl*. Let be the fractional part of */t*, that can be obtained as follows:

$$
\zeta = \frac{\pi}{\Delta t} - \left\lfloor \frac{\pi}{\Delta t} \right\rfloor,\tag{53}
$$

with, 0 < <1. The estimated value of *irfl(t–)* by linear interpolation is thus:

$$\mathbf{\dot{i}}\_{r\|,L}(t-\tau\_k) \equiv \mathbf{\dot{i}}\_{r\|,L}(t-r\Delta t) + \boldsymbol{\zeta}[\mathbf{\dot{i}}\_{r\|,L}(t-k\Delta t) - \mathbf{\dot{i}}\_{r\|,L}(t-k\Delta t-\Delta t)].\tag{54}$$

Figure 7 illustrates the memory buffer management, either for *irfl,0* or for *irfl,L*. At the first simulation time step corresponding to time *t=0t*, calculated *irfl* is stored at memory *1*, and so on until step *Nb* which is the buffer size limit. Beyond this limit, memory cells *1*, *2*, *3* and on, are overwritten as figure. 7 shows, since their previously stored values are not needed any longer.

**Figure 5.** Interpolation scheme: *t* integer multiple.

**Figure 6.** Interpolation scheme: *t* non integer multiple.

**Figure 7.** History buffer management.

be the fractional part of

where a simulation time step

Suppose now that the required value *irfl(t–*

<1. The estimated value of *irfl(t–*

simulation time step corresponding to time *t=0*

available.

with, 0 <

any longer.

**Figure 5.** Interpolation scheme:

**Figure 6.** Interpolation scheme:

*t* integer multiple.

*t* non integer multiple.

samples of *irfl*. Let

the required values must be obtained through interpolation. Figure 6 illustrates this case,

5. Notice that now the required history value, for a time delay of 0.09ms, is not readily

*/*

 

,, , , ( ) ( ) [ ( ) ( )]. *rfl L k rfl L rfl L rfl L* **ii i i** *t t rt t kt t kt t*

Figure 7 illustrates the memory buffer management, either for *irfl,0* or for *irfl,L*. At the first

so on until step *Nb* which is the buffer size limit. Beyond this limit, memory cells *1*, *2*, *3* and on, are overwritten as figure. 7 shows, since their previously stored values are not needed

> *t(ms) irfl,0 irfl,L* 0.0 \* \* 0.03 \* \* 0.09 \* \* 0.12 \* \* 0.15 \* \* 0.18 \* \* 0.21 \* \* 0.24 \* \*

*t(ms) irfl,0 irfl,L* 0.0 \* \* 0.04 \* \* 0.08 \* \* 0.12 \* \* 0.16 \* \* 0.20 \* \* 0.24 \* \*

*t=0.04* ms is assumed instead of the

, *t t*

 *)* is between the *k–th* and the *(k+1)–th* stored

*t*, calculated *irfl* is stored at memory *1*, and

*t*, that can be obtained as follows:

*)* by linear interpolation is thus:

(53)

*t=0.03* ms one at figure

(54)

Linear interpolation is an order 1 numerical procedure and the trapezoidal rule used for the rest of the line model is of order 2. The question arises as to whether or not the order 2 quadratic interpolation should be adopted instead. This has been investigated at (Gutierrez– Robles et al., 2011) and it has been found that the increase in accuracy is marginal.
