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

© 2012 Ramos-Leaños et al., licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. An Advanced Transmission Line and Cable Model in Matlab for the Simulation of Power-System Transients


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 c an be of interest al so to e le ctronic engine ers 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.

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: . .
where V is the vector of voltages, I is the vector of currents, Z and Y are the (N X N) per unitlength 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: 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: 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: where, is the characteristic impedance matrix and its inverse is the characteristic admittance matrix

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: 0 () The application of f() to a square matrix A of order NN as its argument is accomplished as follows: where A 0 is equal to U, the NN unit matrix.Consider now the case of a diagonal matrix: 1 2 00 00 . 00 Application of (8) to  yields: 00 00 00 00 00 00 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 : 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: 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).

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: Then, the value of C1 is determined from ( 13) and ( 14): Expressions ( 5) and ( 6) are applied to the line end conditions at x = L: 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: 00 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:  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: 0, 0 , 0 sh aux 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 , ,, 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).

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.

Rational approximation of Yc
The following rational representation has been proposed for YC in (Marti, 1982) and (Morched et al., 1999): 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).

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): where Hm is a diagonal matrix of the form 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: where Di is the rank-1 matrix obtained from pre-multiplying the i-th column of M by the ith 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: ;1 , 2 , , where is a minimum phase-shift function (Bode, 1945) and i is the time delay associated to the velocity of the i-th mode.Hence: 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: 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): 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 k H  at (34) can now be considered free of delay factors and can be fitted as follows: where Nh(k) is the fitting order for the k-th term k H  , 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 k H  .As (35) is introduced in (33), the following rational form is obtained for H: 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).

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 ,0 0 0 1 where 0 1,2,..., Application of the Inverse Laplace Transform to (37) and ( 38) gives the following continuous-time-state-space (CTSS) form for ish,0: ,0 0 0 1 and 0 1,2,..., The CTSS form to evaluate iaux,0, is derived now.On replacing the fitted form of H given by ( 36) into ( 22 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(Dommel, , 1992)).
Application of this rule to (44) with t as the solution time step results in: where ai=(2+tqi)/ (2-tqi) and i G  =(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:  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: ,0, 0 , 0 , 11 1 Ny 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. In this chapter, the Matlab environment has been chosen for this end.

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(Dommel, & 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.
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 to determine line end (nodal) voltages v0 and vL.
Step 2. Shunt current due to the characteristic admittance of the line is calculated by ( 49) repeated here for convenience: ,00 , 0 sh y aux   iG v i Step 3. Auxiliary current source value, due to the reflected traveling waves at the remote line end, is updated by ( 43): () ,0, 11 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: ,0, 0, 0 2 rfl sh aux

 ii i
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): 00, 0 hist

 iG vi
Step 5. Internal states inside the line model are updated by ( 48) and ( 45): 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 Ntt spans the total simulation time of interest.

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.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.03ms 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.03ms 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: , tt with, 0 < <1.The estimated value of irfl(t-) by linear interpolation is thus: ,, , , 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.

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.

Aerial line case
The transversal geometry of this test case is shown in figure 8. Phase conductors are 1192.5ASCR 54/19 and ground wires are 7 No 5 AWLD.This case consists of three coupled threephase 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.

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

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: 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 unknown parameters.
Consider the following relation (Gustavsen & Semlyen, 1999) An over-determined least squares equation-system is then obtained by evaluating (57) at a number M of specific frequencies, with M>2N: where A is the M2N 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: ( 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): , where W is a diagonal matrix containing previously calculated poles n 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: 1 1 It is clear in (61) that the two polynomials containing the poles n p cancel each other, and that the zeros n z  become the poles of f(s).Notice further that the denominator on the lefthand-side of ( 61) can be written as follows: Zeros n z  are then obtained by finding the roots of (Gustavsen & Semlyen, 1999) which is equivalent to finding the eigenvalues of Q i n ( 6 0 ) ( G u s t a v s e n & S e m l y e n , 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).

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.Semlyen, A., Abdel-Rahman, M. (1982).A state variable approach for the calculation of switching transients on a power transmission line, Circuits and Systems, IEEE Transactions on, vol.29, no.9, (September 1982), pp. (624-633) Wedepohl, L. M., (1965).Electrical characteristics of polyphase transmission systems with special reference to boundary-value calculations at power-line carrier frequencies, Electrical Engineers, Proceedings of the Institution of, vol.112, no.11,(November 1965), pp.(2103-2112)

Figure 1 .
Figure 1.Multi-conductor line segment of length L.
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.

Figure 2 .
Figure 2. Frequency domain circuit representation of a multi-conductor line.
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 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.

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

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

Figure 10 .
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.
Line and Cable Model in Matlab for the Simulation of Power-System Transients 285model results.Notice that both sets of results overlap and that no difference can be seen by eye.Figure13also 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 15 .
Figure 15.Voltage responses at sending and receiving ends.Unit step excitation.

Figure 16 .
Figure 16.Voltage responses at sending and receiving ends.Sinusoidal excitation.
An Advanced Transmission Line and Cable Model in Matlab for the Simulation of Power-System Transients 277

'k,i= xk,i(t-t).
are discrete-state variables and primed variables denote their value at one previous time step xThe discretetime version of (43) maintains its original form: , ki R  =(tRk,i)/ (2-tpk,i).xk,i An Advanced Transmission Line and Cable Model in Matlab for the Simulation of Power-System Transients 279 An Advanced Transmission Line and Cable Model in Matlab for the Simulation of Power-System Transients 281