Open access peer-reviewed chapter

# An Advanced Transmission Line and Cable Model in Matlab for the Simulation of Power-System Transients

By Octavio Ramos-Leaños, Jose Luis Naredo and Jose Alberto Gutierrez-Robles

Submitted: December 2nd 2011Reviewed: April 8th 2012Published: September 26th 2012

DOI: 10.5772/48530

## 1. Introduction

The design and operation of power systems, as well as of power apparatuses, each time depends more on accurate simulations of Electromagnetic Transients (EMTs). Essential to this is to count with advanced models for representing power transmission lines and cables. Electromagnetic Transients Program (EMTP), the most used EMT software, offer various line models. Among these, the most important ones are: 1) the Constant Parameters Line model (CP), 2) the Frequency Dependent or J. Marti Line model (FD) and 3) the Universal Line Model (ULM). The CP Line model is the simplest and most efficient one from the computational point of view. Nevertheless, it tends to overestimate the transient phenomena as it considers that line parameters are constant. Thus, it is recommended only for modeling lines on zones distant to an area where a transient event occurs. The FD Line model (Marti, 1982) evaluates multi–conductor line propagation in the modal domain and takes into account effects due to frequency dependence of the line parameters. Nevertheless, as the transformations between the modal and the phase domains are approximated by real and constant matrices, its accuracy is limited to cases of aerial lines which are symmetric or nearly symmetric. The FD model tends to underestimate the transient phenomena. ULM (Morched et al., 1999) takes into account the full-frequency dependence of line parameters. ULM works directly in phase domain, thus avoiding simplifying assumptions regarding modal–to–phase transformations. So far it is the most general model, capable to accurately represent asymmetric aerial lines as well as underground cables.

The development of ULM is fairly recent and these authors consider that it still is a subject for further research and development. The authors believe also that researchers and power system analysts will benefit considerably from the full understanding of the theoretical basis of the ULM, as well as from counting with a ULM–type code that is easy to understand and modify. One problem with this is that the theoretical basis of ULM includes various topics and subjects that are scattered through several dozens of highly specialized papers. Another difficulty with this is the high complexity of the code for a ULM–type model. This chapter aims at providing a clear and complete description of the theoretical basis for this model. Although this description is intended for power engineers with an interest in electromagnetic transient phenomena, it can be of interest also to electronic engineers involved in the analysis and design of interconnects. The chapter includes as well the description of Matlab program of a ULM–type model, along with executable code and basic examples.

## 2. Multi-conductor transmission line analysis

### 2.1. Telegrapher’s Equations

Electromagnetic behavior of transmission lines and cables is described by the Modified Telegrapher Equations, which in frequency domain are expressed as follows:

dVdx=ZI.E1
dIdx=YV.E2

where V is the vector of voltages, I is the vector of currents, Z and Y are the (N X N) per unit-length series impedance and shunt admitance matrix of a given line with N conductors, repectively. To solve equations (1) and (2), let equation (1) be first differentiated with respect to x; then, (2) is used to eliminate the vector of currents at the right hand side. The resulting expression is a second order matrix ODE involving only unknown voltages:

d2Vdx2=ZYV.E3

In the same way, equation (2) can be differentiated with respect to x and (1) can be used to eliminate the right-hand-side voltage term. The resulting expression involves unknown currents only:

d2Idx2=YZI.E4

Solution to (4) is:

I(x)=C1eYZx+C2eYZx,E5

where C1 and C2 are vectors of integration constants determined by the line boundary conditions; that is, by the connections at the two line ends. In fact, the term including C1 represents a vector of phase currents propagating forward (or in the positive x–direction) along the line, whereas the one with C2 represents a backward (or negative x–direction) propagating vector of phase currents. Expression (5) is an extension of the well–known solution for the single–conductor line. Note that this extension involves the concept of matrix functions. This topic is explained at section 2.2.

The solution to (3) takes a form analogous to (5) and it is obtained conveniently from (5) and (2) as follows:

V(x)=Y1dIdx=Zc[C1eYZxC2eYZx]E6

where,ZC=Y1YZis the characteristic impedance matrix and its inverse is the characteristic admittance matrix YC=YZZ1.

### 2.2. Modal analysis and matrix functions

Matrix functions needed for multi-conductor line analysis are extensions of analytic functions of a one–dimensional variable. Consider the following function and its Taylor expansion:

f(x)=k=0akxkE7

The application of f() to a square matrix A of order NN as its argument is accomplished as follows:

f(A)=k=0akAk,E8

where A0 is equal to U, the NN unit matrix. Consider now the case of a diagonal matrix:

Λ=[λ1000λ2000λN].E9

Application of (8) to yields:

f(Λ)=k=0akΛk=[kakλ1k000kakλ2k000kakλNk]=[f(λ1)000f(λ2)000f(λN)]E10

This expression thus shows that the function of a diagonal matrix is simply obtained applying the one–dimensional form of the function to the matrix nonzero elements. Consider next the function of a diagonalizable matrix A; that is,a matrix A that is similar to a diagonal one :

A=MΛM1.E11

where M is the nonsingular matrix whose columns are the eigenvectors of A, while is the matrix whose diagonal elements are the eigenvalues of A (Strang, 1988).

Application of f() as in (10) to A yields:

f(A)=k=0ak(MΛM1)k=M(k=0akΛk)M1=Mf(Λ)M1E12

Therefore, the function of a diagonalizable matrix is conveniently obtained first by factoring A as in (10), then by applying the function to the diagonal elements of and, finally, by performing the triple matrix product as in (11) and (12).

It is clear from subsection 2.1, that multi-conductor line analysis requires evaluating matrix functions of YZ. To do so, it is generally assumed that YZ always is diagonalizable (Wedephol, 1965; Dommel, 1992). Although there is a possibility for YZ not being diagonalizable (Brandao Faria, 1986), occurrences of this can be easily avoided when conducting practical analysis (Naredo, 1986).

## 3. Line modelling

Figure 1 shows the representation of a multi-conductor transmission line (or cable) of length L, with one of its ends at x = 0 and the other at x = L. Let I0 be the vector of phase currents being injected into the line and V0 the vector of phase voltages, both at x=0. In the same form, IL and VL represent the respective vectors of injected phase currents and of phase voltages at x=L. Line equation solutions (5) and (6) are applied to the line end at x=0:

I0=I(0)=C1+C2E13
V0=V(0)=ZC[C1+C2].E14

Then, the value of C1 is determined from (13) and (14):

C1=I0+YCV02.E15

Expressions (5) and (6) are applied to the line end conditions at x = L:

IL=I(L)=C1eYZLC2eYZLE16

and

VL=V(L)=ZC[C1eYZLC2eYZL].E17

After doing this, (17) is pre–multiplied by YC and subtracted from (16) to obtain:

ILYCVL=2C1eYZL.E18

Finally, (15) is introduced in (18) rendering

ILYCVL=eYZL[I0+YCV0]E19

Expression (19) establishes the relation among voltages and currents at the terminals of a multi-conductor line section. Its physical meaning follows from realizing that the term I0+YCV0 at its right hand side represents a traveling wave of currents leaving the line end at x = 0 and propagating in the positive x–axis direction, whereas IL–YCVL at the left hand side is the traveling wave of currents leaving the line end at x = L.

By a similar process as the previous one for deriving (19), it is possible to show also that the following relation holds as line equation solutions (5) and (6) are applied to line end conditions at x=0:

I0YCV0=eYZL[IL+YCVL]E20

Note however that this relation can also be obtained by simply exchanging at (19) sub–indexes 0 and L. This exchange is justified by the input/output symmetry of the line section. Expressions (19) and (20) provide a very general mathematical model for a multi-conductor transmission line. This is a model based on traveling wave principles. Let (19) and (20) be rewritten as follows:

IL=Ish,LIaux,LE21

where, Ish,L =YCVL is the shunt currents vector produced at terminal L by injected voltages VL. Iaux,L =HIrfl,0 is the auxiliary currents vector consisting of the reflected currents at terminal 0, Irfl,0= I0+ YCV0 and the transfer functions matrix H=e-(YZ)L.

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

I0=Ish,0Iaux,0E22

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

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 auxiliary sources Iaux,0 and Iaux,L.

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 instance, the transformation of (21) to the time domain yields:

i0=ish,0iaux,0E23

with

ish,0=yC*vLE24

and

iaux,0=h*irfl,0E25

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 currents can be represented as

irfl,L=2ish,Liaux,LE26

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 convolutions that are carried out conveniently by means of State–Space methods (Semlyen & Abdel-Rahman, 1982). State–Space equivalents of (23) and (24) arise naturally as YC and H are represented by means of fitted rational functions (Semlyen & Dabuleanu, 1975).

## 4. Phase domain line model

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):

YC=G0+i=1NyGisqiE27

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 single-delay terms is obtained from the following modal factorization (Wedepohl, 1965):

H=MHmM1E28

where Hm is a diagonal matrix of the form

Hm=diag[eγ1L,eγ2L,...,eγNL]E29

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:

H=i=1NDieγiLE30

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γiL=eγ˜iLesτi;i=1,2,,NE31

where exp(γ˜iL)is a minimum phase–shift function (Bode, 1945) and (i is the time delay associated to the velocity of the i–th mode. Hence:

H=i=1NDieγ˜iLesτiE32

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:

H=k=1NgH˜kesτkE33

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

H˜k=i=1IkDieγ˜iLesτik=1,2,...,NgE34

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 termH˜kat (34) can now be considered free of delay factors and can be fitted as follows:

H˜k=i=1Nh(k)Rk,ispk,ik=1,2,...,NgE35

where Nh(k) is the fitting order for the k–th termH˜k, 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 matrixH˜k. As (35) is introduced in (33), the following rational form is obtained for H:

H=k=1Ngesτki=1Nh(k)Rk,ispk,iE36

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 τk 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

Ish,0=G0V0+i=1NyWiE37

where

Wi=GisqiV0i=1,2,...,NyE38

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

ish,0=G0v0+i=1NywiE39

and

dwidt=qiwi+Giv0 i=1,2,...,NyE40

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

Iaux,0=k=1Ngi=1Nh(k)Xk,iE41

with

Xk,i=Rk,ispk,iIrfl,Lesτk;i=1,2,...,Nyk=1,2,...,NgE42

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

iaux,0=k=1Ngi=1Nh(k)xk,iE43
dxk,idt=pk,ixk,i+Rk,iirfl,L(tτk);i=1,2,...,Nyk=1,2,...,NgE44

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:

xk,i=ak,ix'k,i+R˜k,i[irfl,L(tτk)+i'rfl,L(tτk)];i=1,2,...,Nyk=1,2,...,NgE45

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:

iaux,0=k=1Ngi=1Nh(k)xk,iE46

Transmission line simulation of EMTs requires the use of time steps (t smaller than any of the travel times (k 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):

wi=aiw'i+G˜i(v0+v'0);i=1,2,...,NyE47

where ai=(2+(tqi)/ (2-(tqi) and G˜i=((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:

yi=(G˜i1wiv0)/(ai+1);i=1,2,...,NyE48

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

yi=aiy'i+v'0;i=1,2,...,NyE49
ish,0=Gv0+iyaux,0E50

where

iyaux,0=i=1NyG^iyi;G^i=(ai+1)G˜i;
G=G0+i=1NyG˜iE51

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:

i0=Gv0+ihist,0E52

with

ihist,0=iyaux,0iaux,0=i=1NyG^iyik=1Ngi=1Nh(k)xk,iE53

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.

## 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.

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 to determine line end (nodal) voltages v0 and vL.

2. Shunt current due to the characteristic admittance of the line is calculated by (49) repeated here for convenience: ish,0=Gv0+iyaux,0

3. Auxiliary current source value, due to the reflected traveling waves at the remote line end, is updated by (43): iaux,0=k=1Ngi=1Nh(k)xk,i

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: irfl,0=2ish,0iaux,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): i0=Gv0+ihist,0

1. Internal states inside the line model are updated by (48) and (45): yi=aiy'i+v'0xk,i=ak,ix'k,i+R˜k,i[irfl,L(tτk)+i'rfl,L(tτk)]

2. 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 NtΔt 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 τ1, τ2, …, τ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

τmax=max{τ1,τ2,...,τNg}E54

and buffer length Nb is calculated as follows:

Nb=τmaxΔt+1E55

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:

ζ=τΔtτΔt,E56

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

irfl,L(tτk)irfl,L(trΔt)+ζ[irfl,L(tkΔt)irfl,L(tkΔtΔt)].E57

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=0(t, 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.

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.

## 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.

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.

### 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.

 Radius of inner solid conductor (m) 0.015 Resistivity nuclei/sheath (ohm/m) 4.25e-8/2.84e-8 Inner/Outer radius of sheath (m) 0.0258/0.0263 Relative permittivity of 1st & 2nd insulation 2.5

Cable data.

## 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(s)n=1Nrns+p¯nE58

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):

n=1Nrns+p¯nf(s)(1+n=1Nr˜ns+p¯n),E59

where, N is the order of approximation, p¯nrepresents the unknown poles and r^nand r˜nare unknown residues. Poles are initialized by values distributed logarithmically over the frequency range of interest. Expression (56) is now rewritten as follows:

n=1Nrns+p¯n(n=1Nr˜ns+p¯n)f(s)f(s).E60

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

Ax=b,E61

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:

A(j1)x(j)=b,E62

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):

Q=Wgx˜T,E63

where W is a diagonal matrix containing previously calculated polesp¯n, 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:

n=1Nr^ns+p¯nn=1Nr˜ns+p¯n+1=n=1N(s+z^n)n=1N(s+z˜n)f(s).E64

It is clear in (61) that the two polynomials containing the poles p¯ncancel each other, and that the zeros z˜nbecome the poles of f(s). Notice further that the denominator on the left–hand–side of (61) can be written as follows:

n=1Nr˜ns+p¯n+1=n=1N(s+z˜n)n=1N(s+p¯n).E65

Zeros z˜nare then obtained by finding the roots of (Gustavsen & Semlyen, 1999)

i=1N(r˜in=1,niN(s+p¯n))+n=1N(s+p¯n)=0,E66

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).

## 8. Conclusions

Proper design and operation of present-day power systems and apparatuses each time require accurate simulations of their electromagnetic transient performance. An important aspect of these simulations is the realistic representation of transmission lines by digital computer models. The ULM is the most general line model available today, mostly with EMTP–type programs. By being of relatively recent creation, this model still is a subject for substantial improvements in accuracy, stability and computational efficiency. It has been postulated in this work that, both, researchers and power system analysts will benefit considerably from the full understanding of the theoretical basis of the ULM, as well as from counting with a ULM–type code that is easy to understand and modify. It has been contended also that the best way to carry out ULM research and development is by providing a model version in an interpretive environment and Matlab has been the platform chosen for this. This chapter provides a comprehensive description of the theoretical basis of ULM, phase domain line model. In addition to this, full description of a ULM prototype in Matlab has been provided here, along with executable code and typical application examples.

## Appendix

CODE EXECUTION

The following code provides the line model described in the paper and it is embedded into an application example. It simulates the simultaneous energizing of a 150 km long aerial line. At the source side the three voltage sources have a 600 (Thevenin impedance. The program asks for the type of source (unit step or three phase sinusoids). At the load end the line is open. Figure 14 shows the geometry of the simulated line. Figure 15 shows the sending and receiving voltages for the unit step source, while Figure 16 shows the sending and receiving voltages for the sinusoidal source.

Note at Figure 15 that waveforms for phases A and C are equal and their plots are superposed. This is because the symmetry of the line and the excitation.

 clear all clc % m file to set line data LineData % Per unit length parameters [Zg,Zt,Zc,Yg,ZL,YL] = LineParameters(Mu,Eo,Rsu,Geom,Ncon,Ns,w); % Modal Parameters for k=1:Ns [Yc(:,:,k),Vm(k,:),Hm(k,:)] = ABYZLM(ZL(:,:,k),YL(:,:,k),lenght,w(k)); end % Characteristic Admittance Fitting [YcPoles,YcResidues,YcConstant,YcProportional] = YcFit(Yc,f,Ns,Ncon); % Hk fit [HkPoles,HkResidues,HkConstant,HkProportional,md] = HkFitTrace(Hm,Vm,ZL,YL,f,Ns,lenght,Ncon); % m file to execute simulation loop. SimulationLoop Code to Load LineData % Line Geometry % column 1—conductor number % column 2-- x position of each cond in m % column 3-- y position of each cod in m % column 4-- radii of each conductor % column 5-- number of conductor in bundle % column 6-- distance between conductors in bundle % column 7—conductor resistivity % column 8—conductor relative permitivity % column 9-- line lenght in m Geom=[1 0 20 0.0153 3 0.4 2.826e-8 1e3 150e3 2 10 20 0.0153 3 0.4 2.826e-8 1e3 150e3 3 20 20 0.0153 3 0.4 2.826e-8 1e3 150e3]; lenght = Geom(1,9); % Line lenght Ncon = Geom(max(Geom(:,1)),1); % # of cond Rsu = 100; % Earth resistivity Ohm-m Mu = 4*pi*1E-7; % Henry's/meters Eo = (1/(36*pi))*1E-9; % Farads/meters Rhi = 9.09E-7; % Ohm-m resistivity of the iron. Ral = 2.61E-8; % Ohm-m res of the aluminum. Rhg = 2.71E-7; % Ohm-m res of the sky wires. Ns = 500; % Number of samples f = logspace(-2, 6, Ns); % Vector of log spaced Frequencies w = 2*pi*f; % Vector of freqs in radian/sec. Function LineParameters function [Zg,Zt,Zc,Yg,ZT,YT]=LineParameters (Mu,Eo,Rsu,Geom,Ncon,Ns,w) % Function to compute the distances between conductor [Dij,dij,hij]=Height(Geom); Zg = zeros(Ncon,Ncon,Ns); Zt = zeros(Ncon,Ncon,Ns); Zc = zeros(Ncon,Ncon,Ns); Yg = zeros(Ncon,Ncon,Ns); Zcd = zeros(Ncon,Ns); Zaf = zeros(Ncon,Ns); P = (1./sqrt(j*w*Mu/Rsu)); % Complex depth Pmatrix = log(Dij./dij); % Potential Coeff. Matrix Pinv = inv(Pmatrix); % Inverse of Pmatrix % Loop to compute matrices at all frequencies for kl = 1:Ns % Geometric impedance Zg(:,:,kl) = (j*w(kl)*Mu/(2*pi))*Pmatrix; % Earth impedance for km = 1:Ncon for kn = 1:Ncon if km == kn Zt(km,km,kl) = (j*w(kl)*Mu/(2*pi))* log(1+P(kl)./(0.5*hij(km,km))); else num = hij(km,kn)^2 + 4*P(kl)*hij(km,kn) + 4*P(kl)^2 + dij(km,kn)^2; den = hij(km,kn)^2 + dij(km,kn)^2; Zt(km,kn,kl) = (j*w(kl)*Mu/(4*pi))* log(num/den); end end end % Geometric admittance Yg(:,:,kl) = (j*w(kl)*2*pi*Eo)*Pinv; end % Conductor impedance for kd = 1:Ncon; Rcon = Geom(kd,4); % conductor radii in m. Nhaz = Geom(kd,5); % # of conductor in bundle Rpha = Geom(kd,7); % Resistivity in Ohm-m. Zcd(kd,:) = (1/Nhaz)*Rpha./(pi.*Rcon.^2); Zaf(kd,:) = (1/Nhaz)*(1+j).*(1./(2.*pi.*Rcon)) .* sqrt(0.5.*w.*Mu.*Rpha); Zc(kd,kd,:) = sqrt(Zcd(kd,:).^2 + Zaf(kd,:).^2); end % Outputs ZT = Zg + Zt + Zc ; % Total impedance YT = Yg ; % Total admittance Function ABYZLM function [Yc,Vm,Hmo] = ABYZLM(Z,Y,Lo,w) [M, Lm] = eig(Y*Z); % Eigenvalues of YZ Minv = inv(M); % Inverse of eigenvectors matrix Yc = inv(Z)*(M*sqrt(Lm)*Minv); % Characteristic Admittance Gamma = sqrt(diag(Lm)); % Propagation constants. Vm = w./imag(Gamma); % Modal Velocities Hmo = diag(expm(-sqrtm(Lm)*Lo)); % Modal propag. Matrix H Function YcFit function [YcPoles,YcResidues,YcConstant, YcProportional]=YcFit(Yc,f,Ns,Ncon) % Trace of characteristic admittance matrix for k = 1:Ns Ytrace(k,1) = trace(Yc(:,:,k)); end Npol = 6; % Number of poles [Ps] = InitialPoles(f,Npol); % Set initial poles s = j*2*pi*f.'; % Vector of values of variable "s" Ka=2; % 1.-Strictly proper, 2.-Proper, 3.-Improper for khg=1:20 % Fit the trace of Yc (Poles) [YcPoles]=Poles(Ytrace.',s,Ps,Ns,Ka); Ps=YcPoles; end % Residues and constant term for Yc from poles of trace of Yc for k = 1:Ncon for l = 1:Ncon Hs(:,1) = Yc(k,l,:); % k-l term of admittance [C,D,E]=Residue(Hs.',s,YcPoles,Ns,Ka); YcResidues(k,l,:) = C; % k-l residues term YcConstant(k,l) = D; % k-l constant term YcProportional(k,l)=E; %k-l proportional term end end Function HkFitTrace function [HkPoles,HkResidues,HkConstant, HkProportional,md]=HkFit(Hm,Vm,ZL,YL,f,Ns, lenght,Ncon); % Minimum phase of each mode md = ModeDelay(Hm.',f,lenght,Vm.',Ns,Ncon); % Computing Idempotents for k=1:Ns % Function to calculate Idempotents of Y*Z [Hk] = HmIdem(ZL(:,:,k),YL(:,:,k),lenght,f(k), md,Hm(k,:)); HkIdem(:,:,:,k) = Hk; % Idempotents end for m = 1:3 for k=1:Ns TraceHk(m,k) = trace(HkIdem(:,:,m,k)); end end s = j*2*pi*f.'; % Vector of the variable "s" Ka =1;%1.-Strictly proper, 2.-Proper, 3.-Improper Npol = 5; % Number of poles [Ps] = InitialPoles(f,Npol); % Set the initial poles for m = 1:3 Hk = TraceHk(m,:); for khg=1:10 [HkPol]=Poles(Hk,s,Ps,Ns,Ka); Ps=HkPol; end HkPoles(m,:)=Ps; end % Residues for Idempotent matrices of % Hm from the poles of each trace. for m = 1:3 for k = 1:Ncon for l = 1:Ncon Hs(:,1) = HkIdem(k,l,m,:); % k-l term [C,D,E]=Residue(Hs.',s,HkPoles(m,:),Ns,Ka); HkResidues(k,l,m,:) = C; % k-l-m term HkConstant(k,l,m) = D; % k-l-m constant HkProportional(k,l,m) = E; % k-l-m prop end end end SimulationLoop program Ts = 0.016; % Observation time Nt = fix(Ts/Dt); % Number of time steps t1 = fix(md./Dt); % Delay of H in samples t0 = fix(max(md)./Dt); % Maximum time delay as expressed in % number of samples t = (0:Dt:(Nt+t1)*Dt-Dt); % Vector of time Ks = menu('CHOOSE THE TYPE OF INPUT SOURCE' , '1 -unit step' , '2 -sinusoidal'); if Ks == 1 % unit step source Isr = ones(Ncon,Nt+t0); elseif Ks ==2 % sinusoidal source Isr(1,:) = sin(337*t); Isr(2,:) = sin(337*t+2*pi/3); Isr(3,:) = sin(337*t+4*pi/3); end NpYc = length(YcPoles); % Number of poles of Yc NpH = length(HkPoles); % Number of poles for the first % Idempotent matrix Ng = 3; %Number of groups % Initialize the states for both nodes ZA = zeros(Ncon,NpYc); % State variables ZB = zeros(Ncon,NpYc); % State variables YA = zeros(Ncon,NpH,Ng); % State variables YB = zeros(Ncon,NpH,Ng); % State variables IfarA = zeros(Ncon,t0+3); % Current at node A IfarB = zeros(Ncon,t0+3); % Current at node B VO = zeros(Ncon,1); % Voltage at node A Vi = zeros(Ncon,Nt+t0); % Voltage at node A VL = zeros(Ncon,1); % Voltage at node B Vf = zeros(Ncon,Nt+t0); % Voltage at node B IO = zeros(Ncon,1); % Current at node A Ii = zeros(Ncon,Nt+t0); % Current at node A IL = zeros(Ncon,1); % Current at node B If = zeros(Ncon,Nt+t0); % Current at node B Iri = zeros(Ncon,Nt+t0); % Current at Y source Irf = zeros(Ncon,Nt+t0); % Current at Y charge IfarAint = zeros(Ncon,Ng); % Current at node A IfarBint = zeros(Ncon,Ng); % Current at node B % Constants for the state ZA and ZB Ai(:,1) = (1+(Dt/2)*YcPoles)./(1-(Dt/2)*YcPoles); Au(:,1) = ((Dt/2)./(1-(Dt/2)*YcPoles)); Bi(:,1) = (Ai+1).*Au; Gy = zeros(Ncon,Ncon); for nm = 1:NpYc Di(:,:,nm) = YcResidues(:,:,nm)*Bi(nm); Gy = Gy + YcResidues(:,:,nm)*Au(nm); end % Constants for the states YA and YB for k = 1:Ng K1(:,k) = (1+(Dt/2)*HkPoles(:,k))./(1-(Dt/2)*HkPoles(:,k)); Ka(:,k) = (((Dt/2))./(1-(Dt/2)*HkPoles(:,k))); Ku(:,k) = (K1(:,k)+1).*Ka(:,k); end for k = 1:Ng for nm = 1:NpH K2(:,:,nm,k) = HkResidues(:,:,nm,k).*Ka(nm,k); K3(:,:,nm,k) = HkResidues(:,:,nm,k).*Ku(nm,k); end end Gy = Gy + YcConstant; % Admitance of the Ish Yi = diag(eye(3)*[1/600; 1/600; 1/600]); % Admittance of the source, connected at node A Gys = inv(Gy + Yi); % Impedance to calculate VO Yr =diag(eye(3)*[1/1e6; 1/1e6; 1/1e6]); % Admittance of load connected at node B Gyr = inv(Gy + Yr); % Impedance to calculate VL % Contants terms to perform the interpolation tm =md - t1*Dt; % Time for the interpolation % Linear interpolation constants c1 = tm/Dt; c2 = 1-c1; c3 = ones(Ng,1); % Pointers for the interpolation and the buffer h1 = t1+1; h2 = t1+2; h3 = t1+3; for k = t0+2:Nt+t0-3 IfarA(:,1) = IL + Gy*VL + sum(ZB(:,:),2); IfarB(:,1) = IO + Gy*VO + sum(ZA(:,:),2); % Linear interpolation for m = 1:Ng IfarAint(:,m) = c2(m)*IfarA(:,t1(m)) + c3(m)*IfarA(:,h1(m)) + c1(m)*IfarA(:,h2(m)); IfarBint(:,m) = c2(m)*IfarB(:,t1(m)) + c3(m)*IfarB(:,h1(m)) + c1(m)*IfarB(:,h2(m)); end IfarA(:,2:h3) = IfarA(:,1:h2); IfarB(:,2:h3) = IfarB(:,1:h2); for m = 1:NpYc ZA(:,m) = Ai(m)*ZA(:,m) + Di(:,:,m)*VO; ZB(:,m) = Ai(m)*ZB(:,m) + Di(:,:,m)*VL; end for l = 1:Ng for m = 1:NpH YA(:,m,l) = K1(m,l)*YA(:,m,l) + K2(:,:,m,l)*IfarAint(:,l); YB(:,m,l) = K1(m,l)*YB(:,m,l) + K2(:,:,m,l)*IfarBint(:,l); end end HistO = - sum(ZA(:,:),2) + sum(sum(YA(:,:,:),3),2); HistL = - sum(ZB(:,:),2) + sum(sum(YB(:,:,:),3),2); VO = Gys*(Isr(:,k)+HistO); VL = Gyr*HistL; IO = Gy*VO - HistO; IL = Gy*VL - HistL; Vi(:,k) = VO; Vf(:,k) = VL; Ii(:,k) = IO; If(:,k) = IL; end Iri = Yi*Vi; Irf = Yr*Vf; vt = (0:Dt:length(Vi(1,:))*Dt-(t0+4)*Dt)'; N = length(vt); a1 = t1+1; a2 = Nt+t1-3; figure(1),plot(vt,Vi(:,a1:a2),':',vt,Vf(:,a1:a2)) ylabel('Amplitude in volts') xlabel('Time in seconds') legend('Sending end phase A' , 'Sending end phase B' , 'Sending end phase C' , 'Receiving end phase A' , 'Receiving end phase B' , 'Receiving end phase C') Function Height function[Dij,dij,hij]=Height(Geom) Ls = Geom(max(Geom(:,1)),1); Req = zeros(Ls,1); % Equivalent bundle radii k4 = sqrt(2*(Geom(:,6)/2).^2); for nc = 1: Ls; if Geom(nc,5)==1 Req(nc) = Geom(nc,4); else Req(nc) = (Geom(nc,4).*Geom(nc,5).*k4(nc).^ (Geom(nc,5)-1)).^(1./Geom(nc,5)); end end % Direct and image distances among conductors for xl = 1:Ls; for yl = 1:Ls; if xl==yl dij(xl,yl)=Req(xl); y1=Geom(yl,3); hij(xl,yl)=2*y1; Dij(xl,yl)=hij(xl,yl); else x=abs(Geom(yl,2)-Geom(xl,2)); y=abs(Geom(yl,3)-Geom(xl,3)); dij(xl,yl)=sqrt(x^2 + y^2); y1=Geom(xl,3); y2=Geom(yl,3); hij(xl,yl)=y1+y2; x=abs(Geom(yl,2)-Geom(xl,2)); y=hij(xl,yl); Dij(xl,yl)=sqrt(x^2 + y^2); end end end Function InitialPoles function [Ps]=InitialPoles(f,Npol) even = fix(Npol/2); % # of complex initial poles p_odd = Npol/2 - even; % Auxiliary variable to check if number % of initial poles is odd disc = p_odd ~= 0; % 0 for even Nr of initial poles & 1 – for % odd Nr. % Set a real pole in case of disc == 1 if disc == 0 % Even Nr of initial poles pols = []; else % Odd Nr of initial poles pols = [(max(f)-min(f))/2]; end % Set the complex initial poles bet = linspace(min(f),max(f),even); for n=1:length(bet) alf=-bet(n)*1e-2; pols=[pols (alf-j*bet(n)) (alf+j*bet(n)) ]; end Ps = pols.'; % Column vector of initial poles Function Poles function [A]=Poles(Fs,s,Pi,Ns,Ka); Np = length(Pi); % Length of vector containing starting poles CPX = imag(Pi)~=0; % 0 for real pole and 1 for complex pole rp = 0; % Initialize the index for real poles cp = 0; % Initialize the index for complex poles RePole = []; % Initialize the vector of real poles CxPole=[];%Initialize the vector of complex poles % Loop to separate real poles and complex poles for k = 1:Np if CPX(k) == 0 % Real Pole rp = rp + 1; RePole(rp) = Pi(k); elseif CPX(k) == 1 % Complex pole cp = cp + 1; CxPole(cp) = Pi(k); end end Lambda = Pi.'; RePole = sort(RePole); % Sort real poles CxPole = sort(CxPole); % Sort complex poles Lambda = [RePole CxPole]; % Concatenate poles I = diag(ones(1,Np)); % Unit matrix A = []; % Poles B = ones(Ns,1); % the weight factor C = []; % Residues D = zeros(1); % Constant term E = zeros(1); % Proportional term KQA = ones(Ns,1); cpx = imag(Lambda)~=0; % 0 if pole is real and 1 if pole is % complex. dix = zeros(1,Np); % Initializes vector of pole types if cpx(1)~=0 % If the first pole is complex dix(1)=1; % real part dix(2)=2; % imag part k=3; % continue dix for third position else k=2; % If the first pole is real continue dix for the second position end % complete the classification of the poles for m=k:Np if cpx(m)~=0 % If the pole is complex if dix(m-1)==1 dix(m)=2; % If the previous position has the real part put 2 % to identifies the imag part else dix(m)=1; % 1 for the real part of a complex pole end end end % Creates matriz A divided in four parts A = [A1 A2 A3 A4] % A1 = Dk % A2 = B.*ones(Ns,1) % A3 = B.*s % A4 = -Dk*Fs Dk=zeros(Ns,Np); % Initialize matrix with zeros for m=1:Np % Iterative cycle for all poles if dix(m)== 0 % For a real pole Dk(:,m) = B./(s-Lambda(m)); elseif dix(m)== 1 % For the real part Dk(:,m)=B./(s-Lambda(m)) + B./(s-Lambda(m)'); elseif dix(m)== 2 % For the imag part Dk(:,m) = i.*B./(s-Lambda(m-1)) - i.*B./(s-Lambda(m-1)'); end end % Creates work space for matrix A A1 = Dk; A2 = B.*ones(Ns,1); A3 = B.*s; for col = 1:Np A4(:,col) = -(Dk(:,col).*Fs.'); end % Asigns values to A if Ka == 1 A = [A1 A4]; % Strictly proper rational fitting elseif Ka == 2 A = [A1 A2 A4]; % Proper rational fitting elseif Ka == 3 A = [A1 A2 A3 A4]; % Improper rational fitting else disp('Ka need to be 1, 2 or 3') end % Creates matrix b = B*Fs b = B.*Fs.'; % Separating real and imaginary part Are = real(A); % Real part of matrix A Aim = imag(A); % Imaginary part of matrix A bre = real(b); % Real part of matrix b bim = imag(b); % Imaginary part of matrix b An = [Are; Aim]; % Real and imaginary part of A bn = [bre; bim]; % Real and imaginary part of b % Routine to applies the Euclidian norm to An [Xmax Ymax] = size(An); for col=1:Ymax Euclidian(col)=norm(An(:,col),2); An(:,col)=An(:,col)./Euclidian(col); end % Solving system Xn = An\bn; Xn = Xn./Euclidian.'; % Put the residues into matrix C if Ka == 1 C = Xn(Np+1:Ymax); % Strictly proper fitting elseif Ka == 2 C = Xn(Np+2:Ymax); % Proper rational fitting elseif Ka == 3 C = Xn(Np+3:Ymax);% Improper rational fitting else disp('Ka need to be 1, 2 or 3') end % C complex when the residues are complex for m=1:Np if dix(m)==1 alfa = C(m); % real part of a complex pole betta = C(m+1); % imag part of a complex pole C(m) = alfa + i*betta; % the complex pole C(m+1) = alfa - i*betta; % the conjugate end end % Now calculate the zeros for sigma BDA = zeros(Np); KQA = ones(Np,1); % Loop to calculate the zeros of sigma which are the new poles for km = 1:Np if dix(km)== 0 % For a real pole BDA(km,km) = Lambda(km); elseif dix(km)== 1 % For a cp with - imag part BDA(km,km) = real(Lambda(km)); BDA(km,km+1) = imag(Lambda(km)); KQA(km) = 2; Aux = C(km); C(km) = real(Aux); elseif dix(km)== 2 % For a cp with + imag part BDA(km,km) = real(Lambda(km)); BDA(km,km-1) = imag(Lambda(km)); KQA(km) = 0; C(km) = imag(Aux); end end ZEROS = BDA - KQA*C.'; POLS = eig(ZEROS).'; %Forcing (flipping) unstable poles to make them stable uns = real(POLS)"/>0; POLS(uns) = POLS(uns)-2*real(POLS(uns)); % Sort poles in ascending order. First real poles and then complex poles CPX = imag(POLS)~=0; % Set to 0 for a real pole and to1 for a %complex pole rp = 0; % Initialize index for real poles cp = 0; % Initialize index for complex poles RePole = []; % Initialize the vector of real poles CxPole = []; % Initialize the vector of cp % Loop to separate real and complex poles for k = 1:Np if CPX(k) == 0 % Real Pole rp = rp + 1; RePole(rp) = POLS(k); elseif CPX(k) == 1 % Complex pole cp = cp + 1; CxPole(cp) = POLS(k); end end RePole = sort(RePole); % Sort real poles CxPole = sort(CxPole); % Sort complex poles % For conjugate pairs store first the one with positive imag part CxPole = (CxPole.')'; NewPol = [RePole CxPole]; A = NewPol.'; % Output Function Residue function [C,D,E]=Residue(Fs,s,Pi,Ns,Ka); Np = length(Pi); CPX = imag(Pi)~=0; % 0 for a rp and 1 for cp rp = 0; % Initialize the index for real poles cp = 0; % Initialize the index for complex poles RePole = []; % Initialize the vector of real poles CxPole=[]; %Initialize the vector of complex poles % Loop to separate real poles and complex poles for k = 1:Np if CPX(k) == 0 % Real Pole rp = rp + 1; RePole(rp) = Pi(k); elseif CPX(k) == 1 % Complex pole cp = cp + 1; CxPole(cp) = Pi(k); End end RePole = sort(RePole); % Sort real poles CxPole = sort(CxPole); % Sort complex poles CxPole = (CxPole.')'; Lambda = [RePole CxPole]; I = diag(ones(1,Np)); % Unit diagonal matrix A = []; % Poles B = ones(Ns,1); % weight factor C = []; % Residues D = zeros(1); % Constant term E = zeros(1); % Proportional term cpx = imag(Lambda)~=0; % 0 for rp and 1 for cp dix = zeros(1,Np); % Vto identifies poles if cpx(1)~=0 % If the first pole is complex dix(1)=1; % put 1 in dix(1) for the real part dix(2)=2; % put 2 in dix(2) for the imag part k=3; % continue dix for the third position else k=2; % If the first pole is real continue dix for the second % position end % complete classification of the poles for m=k:Np if cpx(m)~=0 % If the pole is complex if dix(m-1)==1 dix(m)=2; % If the previous position has the real part, set to % 2 to identify the imag part else dix(m)=1; % put 1 for the real part of a cp end end end % Output matrices: Dk=zeros(Ns,Np); for m=1:Np if dix(m)==0 % Real pole Dk(:,m) = B./(s-Lambda(m)); elseif dix(m)==1 % Complex pole, 1st part Dk(:,m) = B./(s-Lambda(m)) + B./(s-Lambda(m)'); elseif dix(m)==2 % Complex pole, 2nd part Dk(:,m) = i.*B./(s-Lambda(m-1)) - i.*B./(s-Lambda(m-1)'); end end % Creates work space for matrices A and b AA1=Dk; AA2=B.*ones(Ns,1); AA3=B.*s; if Ka == 1 AA = [AA1]; % Strictly proper rational fit elseif Ka == 2 AA = [AA1 AA2]; % Proper rational fit elseif Ka == 3 AA = [AA1 AA2 AA3]; % Improper fit else disp('Ka must be 1, 2 or 3') end bb = B.*Fs.'; AAre = real(AA); % Real part of matrix A AAim = imag(AA); % Imaginary part of matrix A bbre = real(bb); % Real part of matrix b bbim = imag(bb); % Imaginary part of matrix b AAn = [AAre; AAim]; % Real and imag part of A bbn = [bbre; bbim]; % Real and imag part of b [Xmax Ymax] = size(AAn); for col=1:Ymax Eeuclidian(col)=norm(AAn(:,col),2); AAn(:,col)=AAn(:,col)./Eeuclidian(col); end % Solving system X Xxn=AAn\bbn; X=Xxn./Eeuclidian.'; % Putting residues into matrix C C=X(1:Np); % C is complex when the residues are complex for m=1:Np if dix(m)==1 alfa = C(m); % real part of a complex pole betta = C(m+1); % imag part of a complex pole C(m) = alfa + i*betta; % the complex pole C(m+1) = alfa - i*betta; % the conjugate end end % Outputs if Ka == 1 A = Lambda.'; % Poles C = C; % Residues D = 0; % Constant term E = 0; % Proportional term elseif Ka == 2 A = Lambda.'; % Poles C = C; % Residues D = X(Np+1); % Constant term E = 0; % Proportional term elseif Ka == 3 A = Lambda.'; % Poles C = C; % Residues D = X(Np+1); % Constant term E = X(Np+2); % Proportional term End

Main program

## How to cite and reference

### Cite this chapter Copy to clipboard

Octavio Ramos-Leaños, Jose Luis Naredo and Jose Alberto Gutierrez-Robles (September 26th 2012). An Advanced Transmission Line and Cable Model in Matlab for the Simulation of Power-System Transients, MATLAB - A Fundamental Tool for Scientific Computing and Engineering Applications - Volume 1, Vasilios N. Katsikis, IntechOpen, DOI: 10.5772/48530. Available from:

### chapter statistics

4Crossref citations

### Related Content

Next chapter

#### A New Modeling of the Non-Linear Inductances in MATLAB

By M. Ould Ahmedou, M. Ferfra, M. Chraygane and M. Maaroufi

First chapter

#### MATLAB COM Integration for Engineering Applications

By Mariano Raboso, María I. Jiménez, Lara del Val, Alberto Izquierdo, Juan J. Villacorta and Myriam Codes

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

View all Books