A Time-Domain Explicit Integration Algorithm for Fast Overvoltage Computation of High-Voltage Transmission Line

This paper presents a multiple time-step solver based on a time-domain explicit integration algorithm for improving the computational speed of high-voltage transmission line electromagnetic transient (EMT) simulation. For weakly rigid EMTmodels of high-voltage transmission lines, the previous precise Runge-Kutta integration method with a small time step is adopted; for strongly rigid nonlinear EMTmodels of high-voltage transmission lines, a an improved precise integration method using a large time step is used for the solver of EMTsimulation. Practical simulations for overvoltages of high-voltage transmission line show that multiple time-step EMT solver can be applied to diﬀerent cases of EMT simulation for transmission line.


Introduction
As is known to all, the EMT simulation program (EMTP) [1,2] is widely used to calculate the overvoltage of power systems. Overvoltage calculation often is implemented by solving different ordinary differential or partial differential equations that are used to describe different electrical phenomena or processes. erefore, the numerical algorithms of the differential equation become a key part of the EMT simulation. In EMT-type simulation tools including EMTP, PSCAD/EMTDC [1], NETOMAC, and RTDS, EMTP uses the implicit trapezoidal integration method with second-order accuracy and A-stability to discretize the differential equation of basic electrical components (inductors and capacitors), and the discrete time-domain solution of state variables is obtained by solving the algebraic equation [2,3]. e implicit trapezoidal integration method is limited in integration accuracy and numerical stability. On the one hand, when the power grid topology changes abruptly in EMT simulation, voltage or current simulation waveform will produce a numerical oscillation phenomenon [2,3]. A typical and simple simulation example is the calculation of overvoltage of overhead transmission line without load being switched suddenly. When the trapezoidal method is used to calculate the overvoltage at the end of the transmission line, serious numerical oscillation occurs in the simulation result. is simple case shows that the trapezoidal method cannot accurately simulate the EMT process for the change of topology in power grids. On the other hand, with the continuous access to distributed clean energy, the number of power grid nodes has increased dramatically, resulting in a very large calculation scale for EMT simulation. For such systems, the computation effect of implicit numerical integration algorithms does not perform as fast as explicit integration algorithms. is is mainly because the mathematical model of distributed circuit elements contains many power electronic devices, limiting links, and multiscale transient processes [4,5].
To solve this problem, multirate simulation and parallel computing are used to perform EMTs of large power grids [6,7]. A Parareal [8] parallel computation method is proposed to be applied to EMT simulation, and it is verified that the parallel algorithm can obtain an effective acceleration ratio, but the algorithm cannot suppress the numerical oscillation. Graphics processing units are used to accelerate EMTs of power systems and adopt thread-oriented conversion and automatic code generation technology [9]. Reference [10] proposed the multirate EMTsimulation using latency technology. is method uses the extended implicit trapezoidal integration method to realize the synchronization process between subnets with different simulation steps, but its computational amount is large. Reference [11] proposed the concept of "double-layer grid separation" for the topology characteristics of AC-DC hybrid system. is method greatly improved the simulation efficiency of AC-DC hybrid power grids.
To replace the role of implicit trapezoidal integration in EMT simulation, researchers have begun to seek the application of higher-order and A-stable or L-stable numerical algorithms in EMT simulation. e two-stage diagonally implicit Runge-Kutta method (2S-DIRK) [12] is applied to EMT simulation. e 2S-DIRK is L-stable as same as the implicit Euler method, so it can effectively avoid numerical oscillation [2]. Two-stage three-order Radau IIA method [13] is used for a new parallel method for solving nonlinear EMT state-space equation at multiple time points, and then the Newton method is used to solve it iteratively. e algorithm not only has high computational efficiency but also can avoid numerical oscillation.
Both the precise integration method [14] and the differential quadrature (DQ) method [15] have high accuracy and strong stability, and both of them have been widely used in the field of computational mechanics [16,17]. In order to improve the applicable scope of the precise integration method, the precise Runge-Kutta method is proposed, which is a new algorithm formed by the precise integration method and the explicit Runge-Kutta method with the fourth order [18]. e calculation efficiency and precision of the precise Runge-Kutta method need to be improved, so it cannot be directly used in the EMT simulation of the power systems. In order to improve the EMT simulation efficiency via an improved precise integration method, this paper makes full use of the advantages of the DQ and precise integration method to construct a new method to solve the problem. e improved precise integration method has high accuracy and good stability, and it can use a large time step to improve the simulation efficiency of EMTs. Compared with the implicit trapezoidal integration using a small time step, the calculation error of this method meets the needs of practical engineering calculation and the calculation efficiency is considerable.
Compared with the implicit numerical algorithm, the explicit numerical algorithm has higher execution efficiency. However, the existing explicit methods, such as the explicit Euler method, have too low numerical precision to be directly used in power system transient simulation. Fortunately, the new method proposed in this paper is essentially an explicit single-step high-precision integration method. e idea of applying the improved precise integration method to electromagnetic transient simulation proposed in this paper is the biggest advantage over the existing implicit electromagnetic transient simulation algorithm [19,20] in calculation efficiency, which will be confirmed later in this paper. e remainder of this paper is organized as follows. In Section 2, an improved precise integration method based on DQ with different grids is presented in detail. In Section 3, numerical simulations for overvoltages of high-voltage transmission lines with different time scale and complexity are given to verify the accuracy and effectiveness of the improved precise integration method. Finally, the conclusions derived from the study are presented in Section 4.

Improved Precise Integration Method
Generally, the state-space model of EMTs in power systems can be described as an initial value problem of the first-order ordinary differential equation: where t represents time variable; f(t, x(t)) is a function of time t and state variable x; x 0 is the value of state variable at initial time t 0 . Formula (1) can always be rewritten into ordinary differential equations consisting of homogeneous and nonhomogeneous parts by a certain transformation, and where H ∈ R n×n is the coefficient matrix related to system parameters. To the right of the medium sign in equation (2), the first item is called homogeneous part, and the second item is called nonhomogeneous term. e exact solution of the homogeneous equation of equation (2) on the computer is as follows: In the single-step integration, equation (3) has the following discrete recursive calculation format: e precise integration method (PIM) [14] proposed by Zhong uses additive theorem and incremental storage to calculate exponential matrix T � exp(h × H): where m � 2 N (where N � 15). For very small η � h/m, we can use Taylor series to expand the exponential matrix T in Substituting formula (6) into formula (5) yields According to the matrix multiplication, the following recursive operations are performed iteratively: where T a,0 � T a . At the end of the cycle, In the single integration step [t k , t k+1 ], the solution of the nonhomogeneous dynamic equation (2) can be expressed as follows: 2 Mathematical Problems in Engineering In formula (10), the exponential matrix T � exp(h × H) in the first term on the right side of equation (10) has been obtained by the precise integration method, while the second integration is related to the characteristics of the power system, which is called the Duhamel integration term. Because the calculation of the first term can be achieved by PIM on computer, the numerical error mainly comes from the numerical calculation error of Duhamel integration term. In this paper, Duhamel integration terms are calculated by time-domain differential quadrature scheme based on different grids.

Uniform Grid.
Taking the time-domain differential quadrature method [15] as an example, the integration format of Duhamel integration term in equation (10) is derived. For the initial value problem _ z � w(t, z), the s-th formula of the DQ scheme with s-stage s-order can be expressed as follows: where h is the time step; c j is the grid point; b j are the integration coefficients related to the grid points of DQM. e DQM uses uniform grid, and the formula of Duhamel integration term x k+1 in equation (10) is as follows: where t i � t k + (i × h/s). e explicit numerical method can be used for estimation and then the second equation in formula (11) can be directly used for single-step numerical integration. Specifically, when s � 4, the four-order explicit Runge-Kutta method is used to calculate the estimated value x k+i/s (i � 1, 2, . . . , s) as follows: where S 1 � Hx k + g(t k , x k ); e other values of x k+1/2 , x k+3/4 , and x k+1 are calculated in turn according to formula (13), and the approximate value of Duhamel integration term x k+1 in equation (12) is obtained as follows: where It can be seen from formula (16) that only the first exponential matrix T 3 can be calculated by the PIM in the calculation of this method, so that the amount of calculation will not be too large because of the increase in the number of integration nodes. After finding the approximate value x k+1 of the Duhamel integration term, the approximate value of x k+1 can be obtained by substituting it into equation (10).

Nonequidistant Grids.
e commonly used nonequidistant grids are Legendre grid, Chebyshev grid, and Chebyshev-Gauss-Lobatto grid [21]. Consider the distribution of Chebyshev grid points on regularized interval [0, 1], where Similarly, when s � 4, the four-order explicit Runge-Kutta method is used to calculate the estimated value x k+c i (i � 1, 2, . . . , s) as follows: In equation (18), en the approximation x k+1 of the Duhamel integration term is as follows: where e approximation of x k+1 can be obtained by substituting x k+1 into equation (10). Obviously, using nonuniform grids needs to calculate (s − 1) exponential matrices, which requires more extra calculation than using uniform grid.

Accuracy Test for Improved Precise Integration Method.
As shown in Figure 1, this circuit represents the equivalent schematic diagram of a 515 kV bus switching AC filter in a substation. L s is the equivalent inductance of the AC grid; R s represents the equivalent resistance of the AC grid; e s denotes the equivalent electromotive force of the AC system; α is the initial phase of power supply e s . e capacitor C represents a simplified equivalent circuit of an AC filter bank.
When the circuit breaker K is closed at t � 0 s, the AC filter is fully discharged and there is no current in the branch of inductor. In Figure 1, there are d 0t In general, the system is usually an underdamped or attenuated oscillation system. At this time, the characteristic equation of differential equation (21) has a pair of conjugate complex roots with a negative real part. Under this condition, the voltage across the AC filter can be obtained as where e s (t) � U m sin(ωt + φ s ); ω � 100π rad/s ; In Figure 1, the decay constant of resistance, inductance, and capacitor series circuit is α � R s /(2L s ); resonance angular frequency of this circuit is ω 0 � 1/ ��� � L s C; and natural angular frequency of this circuit is ω α � ������ ω 2 0 − α 2 [22,23].
In this simulation, the lumped electrical parameters are R s � 1.59 Ω, L s � 35.37 mH, and C � 3.124 μF.
e proposed methods were used to solve equation (22) and implemented with the MATLAB scripting language. en, we calculated the voltage of the capacitor C by improved precise integration method (with uniform grid and s � 4), precise Runge-Kutta integration method (four-order explicit Runge-Kutta is used with precise integration method), the two-stage diagonally implicit Runge-Kutta (2S-DIRK with second order in accuracy), and the trapezoidal method (TR) or critical damping adjustment [19]. Figures 2-4 are the computational results of this simulation test calculated by precise differential integration method (PDIM, h � 0.1 ms), precise Runge-Kutta integration method [18,24] (PRKM, h � 0.1 ms), the two-stage diagonally implicit Runge-Kutta (2S-DIRK, h � 0.01 ms), and critical damping adjustment (CDA, h � 0.01 ms) or TR (h � 0.01 ms). And the curve of difference values by these three methods relative to the true solution is also given in Figures 3 and 4. e simulation starts from the zero initial state. e simulation ends at t � 0.1 s.
As shown in Figures 2-4, PDIM has a better computational result compared with TR and 2S-DIRK in the bigger time step. It seems that the PRKM is as accurate as PDIM for this case. And even with large time step, the computational accuracy of PDIM is two orders of magnitude higher than that of the trapezoidal method and 2S-DIRK. In addition, in order to further verify the effectiveness of the proposed method in suppressing numerical oscillations, a nonlinear reactor switching circuit is used for simulation analysis. In Figure 5, the lumped electrical parameters are equivalent resistance, r s � 10 Ω, and equivalent inductor, l s � 20 mH. e piecewise linear current-flux curve of l nl is given in Figure 6.

Single-Phase Transmission Line with Nonlinear Inductance Load.
is case is a high-voltage transmission line with a nonlinear inductance load as shown in Figure 8. In  Table 1. e electrical model used to describe the EMT process of the high-voltage transmission line shown in Figure 8 is the telegraph equation. e telegraph equation is a hyperbolic partial differential equation, and it is necessary to transform it into ordinary differential equations in the form of equation (2) before different numerical methods are used for EMT simulation. We divide the whole transmission line into M sections uniformly as described in Figure 9.
e Π-type cascade equivalent circuit model was carried out on the transmission lines in Figure 8. So, the resistance, inductance, and capacitance of each section are as follows:

Mathematical Problems in Engineering
As shown in Figure 9, it is easy to establish the following first-order linear ordinary differential equations (ODEs) for    Mathematical Problems in Engineering the transmission line using Kirchhoff's voltage and current law [22,23]: Considering the boundary conditions at the head and end of the transmission line in Figure 9, equation (26) is arranged into the following matrix form: where A ∈ R (2M+3)×(2M+3) is the constant sparse matrix; μ(t) is (2M + 3) sparse column vector, which is the excitation source of EMT simulation for the transmission line; .. As shown in Figures 10 and 11, computational results of terminal voltage waveform by the CDA and PDIM are in good agreement. However, CDA is an implicit method, and the Newton-Raphson formula must be used to calculate the voltage waveform of propagative transmission line. During the calculation, the Newton-Raphson formula solves the nonlinear algebraic equations with two iterations. is process is time-consuming when using small step simulations.
e simulation efficiency of the two methods is compared in Table 2. e simulation platform is MATLAB R2012a. e tablet PC processor is AMD Ryzen 5 3500U with Radeon Vega Mobile Gfx 2.10 GHz. e tablet PC uses a 64-bit operating system, and the capacity of RAM is 8 GB.
e simulation starts from the zero initial state except for i L , which is given a small initial value so that the calculation can be performed. e total simulation ends at t � 60 ms. e simulation acceleration ratio of this example is defined as the ratio of the small time-step simulation time of the trapezoidal method (TR) to the large time-step simulation time of PDIM.
As can be observed in Table 2, the PDIM is significantly more efficient at handling nonlinear EMT simulations than the trapezoidal method. Obviously, as the simulation time step of PDIM increases, the corresponding acceleration ratio also changes significantly.
In Figure 12, the calculation result of PDIM using large time step is almost consistent with the simulation waveform of trapezoid method with small time step, which shows that PDIM has good numerical stability and high precision for nonlinear EMT models.

Lightning Overvoltage Calculation of Substation Bus.
is case is a simulation example of lightning tower overvoltage calculation. Figure 13 is a simplified equivalent circuit diagram of lightning overvoltage calculation of substation bus. During lightning stroke, the lightning channel is simulated by resistor parallel ideal current source, and the resistance r e is the resistance of lightning channel; when lightning stroke hits the top of tower, the ideal switch K closes; R ch � 10 Ω represents the impact grounding resistance of tower; the tower is modeled by lossless transmission line, whose wave impedance and wave velocity are 100Ω and 2.7 × 10 8 m/s, respectively, and the transmission line length is L 1 � 50 m. Lightning current is simulated by double exponential wave. Its wavefront time and half peak time are 2.6/50 μs and the peak value of lightning current i L (t) is 100 kA.
In this case, the electrical model of lightning overvoltage simulation for tower is established by using telegraph equation. After spatial interpolation and discretization using the fourth-and second-order interpolation formulas [25,26], the following ordinary differential equations are obtained by taking the number of space segments N � 30: where constant coefficient matrices H ∈ R 61×61 and δ(t) are the input excitation sources of overvoltage at the top of lightning tower.
Equation (31) is solved by PDIM (with Chebyshev grid and s � 4) and the trapezoidal method, respectively. e simulation step of the two methods is h � 0.01 μs. e simulation results are shown in Figure 14. As shown in Figure 14, when t � 4 μs, the voltage value of u(t) at the end of the line is about 984.5 kV and 951.6 kV for the head of the line. In Figures 15 and 16, with the time prolonging, it can be seen that the voltage and current of sending end and receiving end at the end of the transient process are almost the same. And the steady-state current values are almost 99.86 and 98.43 kA, which are the current values at the beginning and end of the transmission line, respectively, which shows the correctness of the simulation results of this case that its real steady-state current value of the lossless transmission line is near 99 kA.
In the calculation of lightning overvoltage, it can be seen from Figure 16 that, because of the fast-changing rate of double exponential lightning current, the two algorithms can only accurately simulate the changing waveform of lightning overvoltage by using smaller simulation steps. As can be observed in difference value of two methods in Figure 15, the simulation results of PDIM and the trapezoidal method are almost similar. is example shows that, for systems with very fast change frequency, PDIM is also competent.

Conclusions
Aimed at the simulation efficiency of EMT simulation for overvoltages of the high-voltage transmission line, an improved precise integration method based on DQM is proposed in this paper. e improved precise integration method inherits the characteristics of high precision and strong stability of the PIM and DQM. And PDIM improves the approximate calculation method of Duhamel integration term in the calculation of nonhomogeneous differential equations by the traditional PIM. Compared with the numerical results of CDA method or the trapezoidal method with small time step, the advantage of PDIM with larger time step is verified in the simulation efficiency of EMT simulation for high-voltage transmission lines.

Data Availability
e data used to support the findings of this study are included within the article.