Implementation of the frequency dependent line model in a real-time power system simulator

In this paper is described the implementation of the frequency-dependent line model (FD-Line) in a real-time digital power system simulator. The main goal with such development is to describe a general procedure to incorporate new realistic models of power system components in modern real-time simulators based on the Electromagnetic Transients Program (EMTP). In this procedure are described, firstly, the steps to obtain the time domain solution of the differential equations that models the electromagnetic behavior in multi-phase transmission lines with frequency dependent parameters. After, the algorithmic solution of the FD-Line model is implemented in Simulink environment, through an S-function programmed in C language, for running off-line simulations of electromagnetic transients. This implementation allows the free assembling of the FD-Line model with any element of the Power System Blockset library and also, it can be used to build any network topology. The main advantage of having a power network built in Simulink is that can be executed in real-time by means of the commercial eMEGAsim simulator. Finally, several simulation cases are presented to validate the accuracy and the real-time performance of the FD-Line model.


Introduction
Real-Time digital power system simulators are a novel technology of recent development based on supercomputers for testing physical equipment such as protective relays and control systems (Dargahi et al., 2012;Dufour et al., 2014;Guillaud et al., 2015;Kuffel et al., 1995).Actually, this technology is spreading rapidly all over the world.This is mainly due to several reasons.Firstly, actually most of the commercial real-time platforms allow users to create new models of power system components as well as control systems for carrying out many studies and analysis in the real-time simulation.The main requirement to validate the performance of physical equipment such as protective relays with HIL tests is to count with realistic models of power system components and control systems.However, it is very common that these models are not yet available for such platforms.As usual, the implementation of these models may also become a challenging task for engineers that do not count with the theoretical background in electromagnetic transient simulations.
In this paper is described the implementation of the frequency-dependent line model (FD-Line) in a real-time digital power system simulator.The idea behind this is to describe a general procedure for incorporating new realistic models of power systems components into modern EMTPbased real-time simulators (Dufour et al., 2010(Dufour et al., , 2011(Dufour et al., , 2014;;Guillaud et al., 2015, Gustavsen et al., 2013).The first step to achieve this goal is to present a comprehensive description of FD-Line model.This means to obtain the time domain solution of differential equations that models the electromagnetic behavior of the transmission line as do the EMTP.Then, the algorithmic solution of the FD-Line model is fully implemented in Simulink for running off-line simulations of electromagnetic transients.This implementation allows the free assembling of the FD-Line model with any element of the Power System Blockset (PSB) library and ARTEMIS software in Simulink as well as to run real-time simulations of electromagnetic transients via RT-Lab platform in the PC-Based real-time simulator eMEGAsim (Dufour et al., 2010(Dufour et al., , 2011(Dufour et al., & 2014)).This paper also provides some test simulation cases of power systems for assessing the accuracy and the real-time performance of the FD-Line model.The parameters used for assessing the real-time performance are the real-time step, which is provided by the simulator, and the accuracy of the transient waveform among off-line and real-time simulations.

FD-Line Model
The FD-Line is well-known as one of the most used and reliable line models in EMTP-based programs (Dommel, 1985;Marti, 1982).This model provides accurate representations for a large class of aerial lines.This section presents a comprehensive frequency and time analysis of the FD-Line equations to obtain the suitable time domain solution, as in the EMTP, for simulating electromagnetic transients in power networks.The algorithmic solution of these equations is implemented in the PSB/Simulink platform to allow the assembling of power networks with other builtin power components.

FD-Line Equations in the Phase Domain
The electromagnetic behaviour of uniform multi-conductor transmission lines and cables is described, in the frequency domain, by the Modified Telegrapher Equations (1a) and (1b) (Dommel, 1985): where V and I are the vectors of voltages and currents at the independent conductors or phases of the line; Z and Y are the matrices of series impedance and parallel admittance, respectively, both are in per unit length; x is the longitudinal distance and w is the angular frequency in radians-per-second.Matrix Z in (1c) is formed by inductive L and resistive R parts In aerial lines the conductive term G usually is negligible and the approximation indicated in (1d) is used almost always.
The traveling wave solution of (1a) and (1b) for a two port representation of a transmission line is given by ( 1e) and ( 1f) and where Yc is shown in (1g) is characteristic admittance matrix and H in (1h) is the matrix of propagation function.It should be note that Y c and H are frequency dependent functions.The methods for computing these functions have been described by Marti (1981), Marti et al. 2017 andDommel (1985).
Traveling wave equations of currents (1e) and (1f) involve matrix relations to 2N c scalar equations, with N c being the number of conductors in the line.Moreover, these line equations are said to be in the phase domain due to all the conductor voltages and currents are coupled.The twoport representation of a multi-conductor transmission line is illustrated in Figure 1.It can be seen that I a and I b are vectors of injected currents at nodes "a" and "b", and V a and V b are vectors of nodal voltages.
The FD-Line model equations were formulated on the basis that a Quasi-TEM mode is propagating along a uniform and imperfect transmission line surrounding by a homogeneous medium (Clayton, 1994;Faria, 1993).It is assumed that the effects of resistive components of transmission lines are small.Also, high-order modes are not considered in

IRACHETA-CORTEZ, FLORES-GUZMAN, AND HASIMOTO-BELTRAN
this line model due to the operation of a transmission line model occurs at low frequencies.

FD-Line Equations in the Modal Domain
The FD-Line model is based on the solution of traveling wave equations (1e) and (1f) in the modal domain.This line model assumes constant and real transformations matrices for decoupling phase-domain line equations.In general, however, these modal transformation matrices are computed as complex functions of frequency.Modal relations between currents and voltages are given by ( 2a) and (2b) (Dommel, 1985;Gustavsen, 2012;Marti, 1981;Wedepohl, 1963) T I is the matrix whose columns are the eigenvectors of YZ.These columns also correspond to the line modes of currents.T V is the matrix of eigenvectors of ZY, or of voltage line modes.The columns of both matrices, T I and T V , can be normalized in such form that the following relation holds: It should be noted that equations (2i) and (2j) represent 2Nc decoupled equations.
As for the shunt currents, I m,sh-a and I m,sh-b , these are defi ned as shown in (2k) and (2l): Figure 2 depicts the frequency-domain representation of the two line terminals with two Norton equivalents.Both equivalents conform a complete single-phase line model.Note from this model that there is no direct connection between the two line ends.Their relation is only through the auxiliary current sources that involve a time delay.The two line ends are thus said to be decoupled in time.

Time Domain Solution
The FD-Line equations previously described need to be solved in the time domain for simulating electromagnetic transients as do the EMTP.Thus, the general time domain form for both, (2g) and (2h) result in (3a) where the lowercase variables are the time domain form of their uppercase counterparts and the symbol * indicates a convolution.The solution of (3a) is carried out very effi ciently through state-space analysis.This is done when Y cm and H m are represented as rational functions in (3b) and ( 3c) Sub-index m in (2a) and (2b) denotes modal quantities.On applying (2a) and (2b) into equations ( 1e) and (1f), relations (2c) and (2d) are obtained: 2c)   and Y cm and H m are the modal matrices of characteristic admittance and of propagation.Since these two matrices are diagonal, relations (2c) and (2d) are said to be uncoupled.
Diagonal matrix H m is related with H and H T as shown in ( 2e) where Λ is the diagonal matrix of eigenvalues for both, ZY and YZ, and ℓ is the length of the transmission line being modeled.
As for the modal characteristic admittance matrix, it is related to Y c as shown in (2f) (Gustavsen, 2012, Wedephol, 1963): where K denotes the fi tting orders, c m,k and r m,k are vectors of residues, p m,k and q m,k are vectors of poles, d m is the vector of residues and τ m is the vector of modal time delays.
The state-space form for i m,sh in (3a) is derived when (3b) is multiplied by V m according to expressions (2k) and (2l).The result is (3d) and where Dt is the fi xed simulation time step, the upper bound of summations K in (3h) and (3j) denote fi tting orders and the primed variables x' mk , x' mk , v' m denote past values at the previous time step t-Δt.Far currents i m,far and i' m,far denote also past values at t-τ m and t-τ m -Δ t , respectively.It is important to note that K in (3h) and ( 3j) is directly related to the number of computations in the line model.Hence, lower values of K result in higher computational effi ciencies.Low fi t orders can be obtained with the Vector Fitting program which is available in the public domain (Gustavsen, 1999(Gustavsen, , 2017;;Kocar et al., 2016;Sintef.no, 2017).
The time domain solution for i m,sh in (3h) is not yet in proper in DSS form since the variable "v m " appears at the right hand side with its time value.The pr oper DSS form is obtained by redefi ning the state variable x m,k = {x m,1 , x m,2 , …, etc.} in (3h) as shown in (3l): By replacing this equation in (3h), expressions (3m) and (3n) are obtained These expressions are now proper for DSS equivalents for i m,sh in (3a).Finally, the FD-Line equations can be written in the following generic EMTP form as shown in (3o) where G m = diag{d m } is the matrix line conductance.Figure 3 provides the DSS circuit representation of the FD-Line model.It should be noticed that, at each side of the line, there is an equivalent modal conductance due to the shunt current term.Next to each conductance there are two sources of history currents.These are due to shunt currents (hist m,sh-a ) and traveling waves of currents originating at the remote line ends (i m,aux ).

RTFD-Line model
The RTFD-Line model is called hereafter the implementation of the frequency-dependent line model in the PSB/ Simulink library for its parallel and real-time execution in a hardwareplatform based on multiprocessors as is shown in Figure 4a ( Dufour et al., 2014;).The state-space form for i m,aux in (3a) is derived when (3c) is multiplied by I m,far according to expressions (2i) and ( 2j).
The result is (3e) Equations ( 3d) and (3e) can be now written in time domain by using the inverse Laplace transform (assuming zero initial conditions) in (3f) and (3g) Expressions (3f) and (3g) are in the continuous-time and their implementation in a digital processor still requires their conversion into a discrete-time form.Discrete-time state-space forms (DSS) of ( 3f) and (3g) are obtained by replacing differentials with a numerical integration rule, such as Euler, Mid-point, Gear, etc.The Mid-point rule is equivalent to trapezoidal integration and, since it is the most used in EMTP-based programs, it is the one adopted here.On applying the trapezoidal rule to differential equations (3f) and (3g), DSS expressions (3h) to (3k) are obtained

Main Features
Figure 4b shows the implementation of the RTFD-Line in PSB/Simulink through two S-Function blocks (Dufour et al., 2014).Each block representing a line-end by a Norton equivalent circuit whose DSS form is shown in Figure 3.The algorithmic solution of each S-Function is programmed in C language for executing faster electromagnetic simulations than the ones performed with interpreted languages such as Matlab M-code.
The auxiliary currents i m,aux , in modal quantities, are the interface variables between the two S-Function blocks of a transmission line.
The main advantage of the RTFD-Line model is that it can be used as a means to split large power networks for their execution in parallel.The resulting power sub-networks are always simulated with the same time step ∆t.For the proper operation of this line model, the modal travel times τ m must be larger than the simulation time step ∆t.More specifi cally, the constraint when the RTFD-Line model is executed sequentially is given by (4a) τ m_series ≥ Δt Thus, let suppose that a given power network is executed in Simulink with a time step of 50 μs.The minimal lengths obtained after replacing this time step in (4d) and (4e) are respectively 15 and 30 km.The algorithmic solution of one FD-Line block is explained as follows: Step I.-The term i m,far is evaluated with (5a) (5a)   where the history current terms hist m-sh and i m,aux must be initialized before starting the simulation and their updating is according to the generic EMTP form (3o).The modal voltages v m are evaluated from phase voltages provided by the external network solver as described in (2b).
Step II.-The term i m,aux-b is evaluated with DSS expressions (3j) and (3k).
The term i' m,far , at the previous time step t-Δt, must be initialized before starting the simulation.
Step III-Saving and updating i m,aux in the history buffer A history buffer, such as the one depicted in Figure 6, is used to save new samples of i m,aux (pointer p1) during a discrete time interval equal to τ m .The oldest samples of i m,aux (pointers p2 and p3) are extracted from the history buffer at each time step.These samples are sent to the opposite FD-Line block for updating i m,far during the next time step as shown in (5a).The length of the buffer must be calculated before starting the simulation as shown in (5b) where M is the total number of modes.It is assumed that modal time delays must be larger than the simulation time step Δt.The process to extract the oldest samples of i m,aux from the history buffer requires interpolations as depicted in Figure 6.This is due to modal time delays τ m are not a multiple integer of the time step Δt.The expression (5c) is a vector of real and fractional numbers.This expression is used to perform interpolations between two consecutive samples of i m,aux in (5d) i m,aux = i m,aux p 2 )   where p 2 and p 3 are pointers for the oldest samples of i m,aux .Three pointers per mode are required for saving and updating the values of i m,aux .All pointers in the history buffer are moving one position at each time step Δt.Each pointer comes back to the initial position after a discrete time interval of length fl oor(τ m /Δt).The process for saving and extracting samples of i m,aux is illustrated in Figure 6.
Step IV-Updating internal DSS state variables w m,k in (3j) and x m,k in (3m).
Step V-The term hist m is calculated according to the generic EMTP form (3o) and is transformed into phase quantities in (5e) Finally, steps I to V are repeated until the simulation is completed.

Interface between the RTFD-Line model and the SSN solver
The State Space Nodal solver (SSN) is a methodology of recent development for the simultaneous interfacing between nodal and state-space equations (Dufour et al., 2010(Dufour et al., , 2011;;Iracheta et al. 2010Iracheta et al. , 2015)).This solver is part of the Opal-RT's Advanced Real-Time Electro-Mechanical Simulation software package (ARTEMIS) and its primarily application is for the use with the Power System Blockset (PSB) or SimPowerSystems.The SSN solver allows the partitioning of electrical networks into state-space groups; that is, in subsystems that can be solved in parallel and for real-time applications.The general equation for interfacing RTFD-Line with SSN is given by (6a): where G is the global conductance matrix, i(t) is the vector of known nodal current injections and the v(t) is the vector of unknown nodal voltages.
The general approach of SSN prescribes that, before solving (6a), each state-space group must have contributed to the assembly of G and the calculation of i(t).
Interface variables from RTFD-Line to SSN must be transformed into phase quantities.The current injection of the transmission line to the SSN is given by (5e).
The contribution of the transmission line conductance to the SSN global conductance matrix is as shown in (6b): where G is the line-model conductance in phase domain and G m is the diagonal matrix of the modal conductance.

Real-Time Simulations
In this section, two confi gurations of lines, a three-phase and a six-phase line, are used to validate the accuracy and the real-time performance of the RTFD-Line model in the eMEGAsim simulator.The target operating system for this simulator is Linux RedHawk 2.6.18.8 (Dufour et al., 2014).

Three-phase line
The transversal geometry used for this transmission line is depicted in Figure 7a and the connection layout in PSB/ Simulink is depicted in Figure 7b.The line length for this example is 191,3 km.The experiment consists in energizing three phases at t = 0 s by means of a 60 Hz AC source with 1 p.u., phase-to-phase rms voltage, and a source resistance of 1e -3 Ω. Phases A and B at the end of the line are opencircuited throughout the simulation time.However, phase C starts in short-circuit at t = 0 s with an RL branch fault of 1 Ω and 1mH.The fault in Phase C is cleared at t = 0,2 s.
The aim of this experiment is to compare execution times among low and high order fi ts.For doing this, the complete three phase network is executed two times in one processor on eMEGAsim.At each time, all transmission lines are executed with the same fi t orders for both Y cm and H m .First, it is considered a low fi t order when the three-phase line is fi tted with 11 poles for Y cm and 6 poles for H m .Second, it is considered a high fi t order when the line is fi tted with 69 poles for Y cm and 67 poles for H m .
The execution time attained with the eMEGAsim simulator for both fi t orders is 6 μs.Basically, there is no noticeable difference between low and high order models due to execution times were rounded to the nearest integer in ms.Table I summarizes the execution times for three-phase and double-circuit lines.The complete power network of Figure 8b is executed in one processor of the eMEGAsim simulator with low and high order fi ts for Y cm and H m .The execution times are summarized in Table I.At this table, the low-order fi t is when the double-circuit line is fi tted with 15 poles for Y cm and 21 poles H m .The high-order fi t is when the line is fi tted with 134 poles for Y cm and 134 poles for H m .The execution times for both low-order and high-order fi ts are respectively 8 and 9 μs.The PSB/Simulink diagram of a three-phase transmission network with three nodes provided in Figure 9 is now used as tested.The transversal geometry of the three power lines is given by Figure 7a.

Double-circuit Line
The network is energized at node 1 with a 60 Hz sinusoidal voltage of magnitude 1 p.u. and its internal source impedance has a magnitude of R = 1 Ω and L = 25 mH.Nodes 2 and 3 are connected to a RL load of magnitude R = 500 Ω and L = 25 mH.Further, a fault module is placed at node 2 for causing switching events during the real-time simulation.
The following test cases are executed on the eMEGAsim simulator for comparing the real-time performance of a three-phase power network.The details for each case are described next: Case A, the three phase RL load is disconnected at node 2 at t = 0 s.Then, a three-phase fault occurs at node 2 at t = 1,05 s for then being cleared at 1,15 s.
Case B, the three phase RL load disconnected at node 2 at t = 0 s.Then, a single-phase fault occurs at node 2 at t = 1,1 s for then being cleared at t = 1,2 s.
Both cases are executed fi rst in a single processor on the eMEGAsim simulator with the minimum time step for achieving its real-time performance.After this, cases A and B are executed with three processors on eMEGAsim as depicted in Figure 9.
Table II shows the time steps attained on the eMEGAsim for these cases.From this table, it follows that the minimum time step to run the complete network is 15,8549 μs that corresponds to a bandwidth of 31,53 kHz.Then, with the parallel execution of this network with three processors, the minimum real-time step is reduced approximately by a factor of three, while the bandwidth is increased by the same factor.Thus, it can be said that parallelization helps to analyze faster transients in large networks.
The waveforms of voltages and currents captured from the real-time execution at node 2, for cases A and B, are depicted in Figures 10 and 11.These waveforms were compared against the off-line P SB/Simulink simulation, with a time step of 16 μs, but no difference was found.This means that the deviations between real-time and off-line simulations are zero.The real-time waveforms match perfectly with the behaviour observed in off-line simulations. a.

Computational times
This section provides several execution times attained for lines with 1, 2, 3 and 6 phases on one processor of eMEGAsim and also, with different fi t orders.These results are summarized in Figure 12 by plotting the execution times in μs versus the sum of poles per mode (Y cm and H m ).In this fi gure, the slopes of execution times for each transmission line are almost constant.There is an execution time difference of around 20 percent between the one with the lines fi tted with the lowest-order and the one with the highest-order.
A detailed analysis is shown in Table III concerning the execution times for lines with 1, 2, 3 and 6 phases when using the RTFD-Line model.The second and third columns of this table specify the computation times for updating the history currents in the S-function.The rates in μs/pole represent the slopes for the execution times depicted in Figure 12.The constant factor in the third column involves the calculations for converting values between modal and phase quantities and vice versa, the linear interpolation, the extracting and updating of history buffers, etc.This factor is obtained when the line is considered as a Constant Parameter Line (CP-Line).This also means that functions Y cm and H m are considered as constant.The fourth column shows the execution time required by the external network solver (SSN).Clearly, it can be observed that the execution times are almost proportional to the number of phases.
The following empirical formula is derived from the results provided in Table III for estimating the execution times in multi-phase lines for several order fi ts: where Nph is the number of phases, NpolY and NpolH are the number of poles for Y cm and H m , respectively.smaller time steps.This results in higher effective bandwidths of frequency for analyzing more accurately EMTs.
The magnitudes of voltage and current waveforms captured directly from the real-time simulation on eMEGAsim are in good agreement with those obtained from off-line simulations in PSB/Simulink simulation.The effective bandwidths obtained from the minimal execution times attained by the RT-FD-Line model on eMEGAsim simulator are adequate to simulate real-time simulations of switching and faults occurring in power networks.

Conclusions
A comprehensive description of the real-time version of the frequency-dependent line model (RTFD-Line) model has been presented in this paper.The RTFD-Line model can be used for running off-line simulations of electromagnetic transients (EMTs) in the PSB/Simulink toolbox or for running real-time simulations through the distributed realtime platform (RT-LAB) on the eMEGAsim simulator.In both cases, the frequency-dependence of line parameters is included for simulating EMTs in power systems.
The main contribution of this paper has been the fully integration of the RTFD-Line model within the Simulink environment.This implementation allows the free assembling of the line model with any power network topology and with other built-in power elements of the PSB/ Simulink toolbox.Moreover, through the parallelization of this line model, it is possible to execute more efficient offline simulations of EMTs in the PSB/Simulink toolbox.
The main feature of the RTFD-Line is the parallelization of power networks to execute real-time simulations with 2f)Norton equivalents for the two line ends are derived from (2c) and (2d) as the sum of vectors of two currents in (currents at each line-end, I m,aux-a and I m,aux-b , are defi ned in terms of traveling waves originating at the respective opposite line-end; that is shown in (

Figure 1 .
Figure 1.Two-port representation of a multiconductor transmission line.Source: Authors.

Figure 2 .
Figure 2. Circuit representation of the FD-Line model in the frequency-domain.Source: Authors.
hand, when the RTFD-Line model is executed in parallel it is necessary a second time step ∆t in (4a) to compensate the time delay generated by the parallelization of power sub-networks on the eMEGAsim simulator.So, expression (4a) becomes (4b)τ m_parallel ≥ 2Δt (4b)By considering the fact that modal velocities in aerial lines are relatively close to the speed of light, it can be assumed that the length of the line ℓ is related to (4a) and (4b) as shown in (4c) 108 m/s is the speed of light.Now, from expressions (4a), (4b) and (4c) one can state that the minimal lengths for the proper operation of the RTFD-Line model, both sequentially and parallel, are defi ned by (4d) and (4e)

Figure 3 .
Figure 3. Norton equivalent of the FD-Line model in the time domain.Source: Authors

Figure 4 .
Figure 4. RTFD-Line model in PSB/Simulink.Source: Authors Algorithmic solution of the RTFD-Line modelFigure 5 depicts a fl owchart with the algorithmic solution of the RTFD-Line model for its parallel execution in a PC-cluster.The parallelization of this line model is made possible by the natural time decoupling between the two Norton equivalents of each line end.In this manner, each Norton block can be assigned to a different processor.

Figure 5 .Figure 6 .
Figure 5. Flowchart with the algorithmic solution of the FD-Line model.Source: Authors

Figures
Figures 8a and 8b provide the line confi guration and the power network layout of a double-circuit transmission line in PSB/Simulink.

Figure 7 .*
Figure 7. Three-phase circuit line, a) Line confi guration, b) Power network layout in PSB/Simulink.Source: Authors Table I.Execution Times Attained on eMEGAsim with the RTFD-Line Model Number of phases

Figure 8 .
Figure 8. Double-circuit line in PSB/Simulink, a) Line confi guration, b) Power network layout.Source: Authors Three phase transmission network: single processor vs. multi-processor

Figure 10 .
Figure 10.Real-time simulation for case A. a) Voltage waveforms at node 2, b) Current waveforms.Source: Authors The accuracy of the transient responses captured from the real-time execution on eMEGAsim simulator is compared against the off-line PSB/Simulink simulation.

*
Execution times are the minimum time steps for running the simulation in eMEGAsim.

Figure 11 .
Figure 11.Real-time simulation for case B. a) Voltage waveforms at node 2, b) Current waveforms at the RL load.Source: Authors

Table 2 .
Execution Times Attained on the eMEGAsim Simulator

Table 3 .
Analysis of Execution Times in the RTFD-Line Model

of phases Computational time RTFD-Line (S-function) Network Solver (SSN) Total Time
* Total number of poles is the sum of order fits for Ycm and Hm.Source: Authors