Extended Critical Damping Adjustment Method for Electromagnetic Transients Simulation in Power Systems

Qiang Li 1, Yong Wang 1,2,* ID , Huaxing Rao 1, Lei Zhang 1, Jing Ye 1, Xiaolin Lei 1 and Bingwen Chen 1 1 College of Electrical Engineering & New Energy, China Three Gorges University, Yichang 443002, China; 2016120506258@ctgu.edu.cn (Q.L.); 2016120506186@ctgu.edu.cn (H.R.); leizhang3188@163.com (L.Z.); yejing2000310@163.com (J.Y.); 2016120506216@ctgu.edu.cn (X.L.); 201708520721002@ctgu.edu.cn (B.C.) 2 China Shanghai Administration Bureau of State Grid Operation Company, Shanghai 200000, China * Correspondence: yongwang2015wy@163.com; Tel.: +86-021-5172-0605


Introduction
The phenomenon of numerical oscillations is a common problem in the electromagnetic transient simulation in power systems [1], which is due to the reason that the numerical integration algorithm in power-system transient simulation calculation is mainly an implicit trapezoidal integration method [1,2].However, the implicit trapezoidal integration method is not L-stable and does not have the ability to eliminate numerical oscillations [3].Therefore, the traditional power-system electromagnetic transient simulation software, such as EMTP 3.0, use critical damping adjustment (CDA) to eliminate numerical oscillations [4,5].
The CDA method is essentially a combined integration algorithm.The core step of the CDA method is to use two-half-step backward Euler method for numerical integrations to eliminate numerical oscillations [6,7].The backward Euler method is a one-step, low-order, L-stable numerical integration algorithm, which eliminates numerical oscillations from its source.
In addition, other schemes for solving the problem of numerical oscillations indirectly use different technical approaches to approximating non-state variables for obtaining their "real" values [1,5].For example, PSCAD version 4.0 uses a half-step interpolation algorithm to eliminate numerical oscillations, but this method is only suitable for offline simulations [8].
Based on the state space method in electromagnetic transient simulations [9,10], we proposed an extended critical damping adjustment method (ECDA), in this paper, to solve the problem of numerical oscillations directly from the origin of the numerical oscillations.The core idea of the ECDA is similar to that of CDA.The difference between these two methods lies in the use of a new class of low-order linear multistep formulae (L 2 MF) to avoid the numerical oscillations.The absolute value of the main coefficient of the local truncation error of a new four-step L 2 MF using Chebyshev grid points was smaller than that of the backward Euler method.The pure L 2 MF method and the implicit trapezoidal integration method were used to form the ECDA to solve numerical examples.The simulation results show that the L 2 MF method can effectively avoid the numerical oscillations, and also indirectly show that the ECDA has a better simulation effect than the CDA in theory, and hence has the potential applications in EMTP-type programs.

Construction of a New Class of Linear Multistep Formulae by Classic Differential Quadrature
The initial-value problem of an ordinary differential equation was shown in the following: .
where t is time and x is a one-dimensional variable; f (t, x) is a function of the time t and the variable x; and x 0 is the initial value of an unknown state variable.One of main methods for solving the initial-value problem, as indicated in Equation ( 1), is the use of classical linear multistep formulae (LMF) [11], and one kind of the LMF is the following backward differentiation formula (BDF), which was described as following: where s is the number of steps; α i is integration coefficient; h denotes the integration step; x n+i ≈ x(t n+i ) is the approximate value of a state variable at a time point t n+i ; and f n+s = f (t n+s , x n+s ).Obviously, when s takes one, Equation ( 2) is the first-order backward Euler method that is A-stable and L-stable.Based on the traditional BDF, new linear multistep methods were constructed as follows: where α i and β s are integration coefficients of new LMF, respectively.In order to obtain the unknown integration coefficients of new LMF, we similarly used s-stage s-order time-domain differential quadrature method (DQM) [12][13][14] to solve Equation (1).For any time-step (h = t n+1 − t n ), we have the following equation: . . .
where c i = ( t i − t n )/h, i ∈ (0, s) is the non-equidistant grid point, which is distributed over the regularization interval [0,  where G is called the weighted coefficient matrix of DQM, and the calculation formulae of each element in G is given in Reference [15] G 0 = −Ge, where e is the column vector with each element equaling to one.In Equation ( 4), the grid point distribution has to be a non-equidistant grid, including Chebyshev grid, Chebyshev-Gauss-Lobatto grid, and Legendre grid [16].In this paper, we used the DQM with Chebyshev grid points [16] for constructing the new LMF.Suppose: where α s = 1.
As an example, the new three-step LMF using Chebyshev grid points was shown below: The new four-step LMF using Chebyshev grid points is as follows: In addition, the backward Euler method is obtained when s takes one in Equation ( 3) regardless of the kind of the used grid.Similarly, when s takes two, we get the second-order Gear method.

The Compatibility and Convergence of the L 2 MF
The new linear multistep methods have practical applications only if they satisfy compatibility and convergence.Therefore, according to Equation (3), the first-characteristic polynomial and the second-characteristic polynomial of z were as follows [17]: where ρ(1) = 0 and ρ (1) = σ (1).Therefore, the new linear multistep methods satisfy the compatibility.Suppose ρ(z) = 0.A unique root modulus is |z 1 | = 1, and the rest of the roots' modulus satisfies |z i | < 1, i ∈ (2, s).Therefore, it also meets the root condition.Since the new linear multistep formulae satisfy both the compatibility and the root condition, the new linear multistep formulae are convergent according to the Dahlquist equivalence theorem [17,18].

Stability of the L 2 MF
For the stability analysis of the new linear multistep method, the following standard test equations can be considered: where Re(λ) < 0. According to Equation (3), the corresponding characteristic equation was obtained: where the stability domains of different s values in Equation (3), supposing, were shown in Figures 1 and 2. In Figures 1 and 2, the outside regions of boundary curves were the stability domain of the new LMF.Obviously, the numerical stability regions of the new linear multistep methods contained the entire left-half complex plane, so the new LMF were A-stable.
According to the Dahlquist barriers theorem [19], new two-step LMF were second-order methods, and other LMF were first-order methods.The main coefficient of the local truncation error for the three-step LMF and the four-step LMFs are −0.6683 and −0.3009, respectively.
where the stability domains of different s values in Equation ( 3), supposing, were shown in Figures 1 and 2. In Figures 1 and 2, the outside regions of boundary curves were the stability domain of the new LMF.Obviously, the numerical stability regions of the new linear multistep methods contained the entire left-half complex plane, so the new LMF were A-stable.
According to the Dahlquist barriers theorem [19], new two-step LMF were second-order methods, and other LMF were first-order methods.The main coefficient of the local truncation error for the three-step LMF and the four-step LMFs are 6683 .0 − and 3009 .0 − , respectively.In addition, Equation ( 9) was substituted into Equation (3); hence the following equation was obtained: . Therefore, the new linear multistep formulae were L-stable [19].
In summary, the new linear multistep methods are low-order methods.So it is named as low-order LMF (L2MF) method.In numerical stability, the L2MF method was both A-stable and L-stable.Therefore, compared with the A-stable second-order implicit trapezoidal integration where the stability domains of different s values in Equation (3), supposing, were shown in Figures 1 and 2. In Figures 1 and 2, the outside regions of boundary curves were the stability domain of the new LMF.Obviously, the numerical stability regions of the new linear multistep methods contained the entire left-half complex plane, so the new LMF were A-stable.
According to the Dahlquist barriers theorem [19], new two-step LMF were second-order methods, and other LMF were first-order methods.The main coefficient of the local truncation error for the three-step LMF and the four-step LMFs are 6683 .0 − and 3009 .0 − , respectively.In addition, Equation ( 9) was substituted into Equation (3); hence the following equation was obtained: . Therefore, the new linear multistep formulae were L-stable [19].
In summary, the new linear multistep methods are low-order methods.So it is named as low-order LMF (L2MF) method.In numerical stability, the L2MF method was both A-stable and L-stable.Therefore, compared with the A-stable second-order implicit trapezoidal integration In addition, Equation ( 9) was substituted into Equation (3); hence the following equation was obtained: where |x n+s /x n+s−1 | → 0 if Re(λh) → −∞ .Therefore, the new linear multistep formulae were L-stable [19].In summary, the new linear multistep methods are low-order methods.So it is named as low-order LMF (L 2 MF) method.In numerical stability, the L 2 MF method was both A-stable and L-stable.Therefore, compared with the A-stable second-order implicit trapezoidal integration method, the L 2 MF method has strong damping characteristics, which can effectively suppress the numerical oscillations, and hence is suitable for electromagnetic transient simulations in power systems.

Extended Critical Damping Adjustment Scheme
From the above analysis, it can be seen that the L 2 MF method belongs to the numerical methods with low-order and strong numerical stability.For this purpose, it is necessary to combine the implicit trapezoidal integration method, which is widely used in electromagnetic transient calculation, to calculate the implicit trapezoidal integration method under normal circumstances, and to switch to L 2 MF only when the system is mutated, such as switch operation and line fault.In practice, we can use the same method of auxiliary detection mutation as the CDA method and the algorithm switching strategy.Because the absolute value of the local truncation error of the four-step L 2 MF method is smaller than that of the backward Euler method (the absolute value of the main coefficient of local truncation error is 0.5), it was recommend that a four-step L 2 MF method and an implicit trapezoid integration method are used to form that combined integration algorithm for the electromagnetic transient simulation in power systems.
In the new combination integration scheme, the implicit trapezoidal integration method had the ability of self-starting computation, and the L 2 MF method did not have the self-starting calculation ability, which needed a startup algorithm for independent calculation.The starting numerical method used in this paper was explicit Euler method (forward Euler method).Figure 3 is a flow chart of electromagnetic transient simulation in power systems on a basis of extended CDA method (ECDA).
numerical oscillations, and hence is suitable for electromagnetic transient simulations in power systems.

Extended Critical Damping Adjustment Scheme
From the above analysis, it can be seen that the L2MF method belongs to the numerical methods with low-order and strong numerical stability.For this purpose, it is necessary to combine the implicit trapezoidal integration method, which is widely used in electromagnetic transient calculation, to calculate the implicit trapezoidal integration method under normal circumstances, and to switch to L2MF only when the system is mutated, such as switch operation and line fault.In practice, we can use the same method of auxiliary detection mutation as the CDA method and the algorithm switching strategy.Because the absolute value of the local truncation error of the four-step L2MF method is smaller than that of the backward Euler method (the absolute value of the main coefficient of local truncation error is 0.5), it was recommend that a four-step L2MF method and an implicit trapezoid integration method are used to form that combined integration algorithm for the electromagnetic transient simulation in power systems.
In the new combination integration scheme, the implicit trapezoidal integration method had the ability of self-starting computation, and the L2MF method did not have the self-starting calculation ability, which needed a startup algorithm for independent calculation.The starting numerical method used in this paper was explicit Euler method (forward Euler method).Figure 3  The combined integration algorithm proposed in this paper was a class of algorithms, and the combination algorithm was the CDA method when s takes one.Therefore, the new combination algorithm proposed in this paper was also named as the extended CDA method.The combined integration algorithm proposed in this paper was a class of algorithms, and the combination algorithm was the CDA method when s takes one.Therefore, the new combination algorithm proposed in this paper was also named as the extended CDA method.

Simulation of Half-Wave Rectifier Circuit
The example is half-wave rectifier circuit.As shown in Figure 4, the excitation source was an ideal sinusoidal voltage source E(t) = 100 √ 2 sin(100πt).The value of the DC voltage source is 20 V .The first-order differential equation for the circuit is as follows: where i is the current, R is the resistance, and L is the inductance.A three-step L 2 MF method and an implicit trapezoidal integration method were used to solve Equation (13), and the integration step in this case was set to h = 0.1 ms.The simulation results are shown in Figure 5.

Simulation of Half-Wave Rectifier Circuit
The example is half-wave rectifier circuit.As shown in Figure 4, the excitation source was an ideal sinusoidal voltage source . The value of the DC voltage source is 20 V .The first-order differential equation for the circuit is as follows: where i is the current, R is the resistance, and L is the inductance.A three-step L2MF method and an implicit trapezoidal integration method were used to solve Equation (13), and the integration step in this case was set to 0.1 ms h = . The simulation results are shown in Figure 5.As shown in Figure 5, since the three-step L2MF method was L-stable and had a strong damping characteristic, which could suppress the numerical oscillations of non-state variables when the diode D turned on.In contrast, the implicit trapezoidal integration method used in the style of EMTP produces a serious numerical oscillation when the network changed suddenly.It is worth mentioning that the implicit trapezoidal integration method directly used for solving Equation (13) did not produce numerical oscillations.

Simulation of Half-Wave Rectifier Circuit
The example is half-wave rectifier circuit.As shown in Figure 4, the excitation source was an ideal sinusoidal voltage source . The value of the DC voltage source is 20 V .The first-order differential equation for the circuit is as follows: where i is the current, R is the resistance, and L is the inductance.A three-step L2MF method and an implicit trapezoidal integration method were used to solve Equation (13), and the integration step in this case was set to 0.1 ms h = . The simulation results are shown in Figure 5.As shown in Figure 5, since the three-step L2MF method was L-stable and had a strong damping characteristic, which could suppress the numerical oscillations of non-state variables when the diode D turned on.In contrast, the implicit trapezoidal integration method used in the style of EMTP produces a serious numerical oscillation when the network changed suddenly.It is worth mentioning that the implicit trapezoidal integration method directly used for solving Equation (13) did not produce numerical oscillations.As shown in Figure 5, since the three-step L 2 MF method was L-stable and had a strong damping characteristic, which could suppress the numerical oscillations of non-state variables when the diode D turned on.In contrast, the implicit trapezoidal integration method used in the style of EMTP produces Appl.Sci.2018, 8, 883 7 of 13 a serious numerical oscillation when the network changed suddenly.It is worth mentioning that the implicit trapezoidal integration method directly used for solving Equation (13) did not produce numerical oscillations.

Simulation of High-Voltage Transmission Line without Load Being Switched-In Suddenly
This example is an electromagnetic transient simulation of a single-phase middle-long transmission line with distributed parameters.The initial phase of the sinusoidal voltage source was ϕ 0 = 90 • .The other parameters of the transmission line were given in Figure 6.This example is an electromagnetic transient simulation of a single-phase middle-long transmission line with distributed parameters.The initial phase of the sinusoidal voltage source was 0 90 ϕ = °.The other parameters of the transmission line were given in Figure 6.The mathematical model of the electromagnetic transient process for the transmission line was described by the telegraph equation as follows: where y and t are the spatial variables of the length and the time variable, respectively; ( , ) v y t and ( , ) i y t denote the changes of the voltage and the current with respect to the space and the time, respectively; 0 l , 0 c , 0 r and 0 g are the unit parameters of inductance, capacitance, resistance, and susceptance, respectively.In this paper, we adopt the time-domain method to settle the telegraph Equation (14).First, in order to obtain the voltage and current variables along the transmission line, we decompose the whole line into N sub-intervals uniformly, as is shown in Figure 7. Supposing the total length of the transmission line is l , we define: the length of each small interval is the voltage at the point, where is the current from the point to the point k y ; and 0 i and are the currents of two terminals in the line, respectively.
Furthermore, we apply the fourth-order accuracy interpolation formula to do spatial discretization of the telegraph equation [20]: where λ denotes ν or i ; 1 θ and 2 θ are the constant coefficients decided by the interpolation formula; and Substituting Equation ( 15) into the telegraph equation (Equation ( 14)), we obtained the spatial discretization form of the telegraph equation as follows: The mathematical model of the electromagnetic transient process for the transmission line was described by the telegraph equation as follows: where y and t are the spatial variables of the length and the time variable, respectively; v(y, t) and i(y, t) denote the changes of the voltage and the current with respect to the space and the time, respectively; l 0 , c 0 , r 0 and g 0 are the unit parameters of inductance, capacitance, resistance, and susceptance, respectively.In this paper, we adopt the time-domain method to settle the telegraph Equation ( 14).First, in order to obtain the voltage and current variables along the transmission line, we decompose the whole line into N sub-intervals uniformly, as is shown in Figure 7. Supposing the total length of the transmission line is l, we define: the length of each small interval is ∆y = l/N.Let v k denote the voltage at the point, where y = y k = k • ∆y, k ∈ (0, N); i k , k ∈ (1, N) is the current from the point y k−1 to the point y k ; and i 0 and i N+1 are the currents of two terminals in the line, respectively.Furthermore, we apply the fourth-order accuracy interpolation formula to do spatial discretization of the telegraph equation [20]: where λ denotes ν or i; θ 1 and θ 2 are the constant coefficients decided by the interpolation formula; and Substituting Equation ( 15) into the telegraph equation (Equation ( 14)), we obtained the spatial discretization form of the telegraph equation as follows: For the head-end and the tail-end of the line, we used the second-order accuracy interpolation formula [20] to do spatial discretization as follows: where Considering the boundary conditions of both ends of the transmission line, we obtained the following equations by combining Equations ( 16)-( 18): where K, ρ ∈ R (2N+1)×(2N+1) is the constant coefficient matrix; ω(t) is the input excitation source for the transmission line; and where: 20) where The four-step L2MF method and the implicit trapezoidal integration method were used to solve Equation (19), and the integration step of the simulation test was 1 s h = μ for both L2MF and the implicit trapezoidal integration method.Moreover, the transmission line was uniformly divided into N = 50 sections .The simulation results are shown in Figure 8, and Figure 9 is an enlargement of Figure 8.As shown in Figures 8 and 9, the calculation results of the four-step L2MF method were basically consistent with the second-order the implicit trapezoidal integration method.In the moment of the circuit closing, the implicit trapezoidal integration method had a continuous numerical oscillation due to the lack of L-stability.On the contrary, the four-step L2MF method was L-stable, and therefore had strong damping property, which effectively inhibited the numerical oscillations.The four-step L 2 MF method and the implicit trapezoidal integration method were used to solve Equation (19), and the integration step of the simulation test was h = 1 µs for both L 2 MF and the implicit trapezoidal integration method.Moreover, the transmission line was uniformly divided into N = 50 sections.The simulation results are shown in Figure 8, and Figure 9 is an enlargement of Figure 8.The four-step L2MF method and the implicit trapezoidal integration method were used to solve Equation (19), and the integration step of the simulation test was 1 s h = μ for both L2MF and the implicit trapezoidal integration method.Moreover, the transmission line was uniformly divided into N = 50 sections .The simulation results are shown in Figure 8, and Figure 9 is an enlargement of Figure 8.As shown in Figures 8 and 9, the calculation results of the four-step L2MF method were basically consistent with the second-order the implicit trapezoidal integration method.In the moment of the circuit closing, the implicit trapezoidal integration method had a continuous numerical oscillation due to the lack of L-stability.On the contrary, the four-step L2MF method was L-stable, and therefore had strong damping property, which effectively inhibited the numerical oscillations.The four-step L2MF method and the implicit trapezoidal integration method were used to solve Equation (19), and the integration step of the simulation test was 1 s h = μ for both L2MF and the implicit trapezoidal integration method.Moreover, the transmission line was uniformly divided into N = 50 sections .The simulation results are shown in Figure 8, and Figure 9 is an enlargement of Figure 8.As shown in Figures 8 and 9, the calculation results of the four-step L2MF method were basically consistent with the second-order the implicit trapezoidal integration method.In the moment of the circuit closing, the implicit trapezoidal integration method had a continuous numerical oscillation due to the lack of L-stability.On the contrary, the four-step L2MF method was L-stable, and therefore had strong damping property, which effectively inhibited the numerical oscillations.As shown in Figures 8 and 9, the calculation results of the four-step L 2 MF method were basically consistent with the second-order the implicit trapezoidal integration method.In the moment of the circuit closing, the implicit trapezoidal integration method had a continuous numerical oscillation due to the lack of L-stability.On the contrary, the four-step L 2 MF method was L-stable, and therefore had strong damping property, which effectively inhibited the numerical oscillations.

Very Fast Transient Overvoltage Calculation
In this case, we simulate a situation when the disconnector of 550 kV gas insulated switchgear (GIS) acts on the no-load bus.During this process, very fast transient overvoltage (VFTO) is caused by slow contact action and poor arc-extinguishing performance of switch [21,22].Figure 10 is a simplified equivalent circuit diagram for this process.In Figure 10, e(t) is the AC power supply; L t and C t are used to simulate the lumped inductance and capacitance parameters of the transformer and other equipment.DS is the isolation switch, and its mathematical model can be established by a nonlinear time-varying resistance.C r is the lumped ground capacitance with initial voltage to simulate the residual charge; L 1 and L 2 are import and export short bus, respectively.All parameters of this example in Figure 10 were taken from Reference [23].
Appl.Sci.2018, 8, x FOR PEER REVIEW 10 of 13 In this case, we simulate a situation when the disconnector of 550 kV gas insulated switchgear (GIS) acts on the no-load bus.During this process, very fast transient overvoltage (VFTO) is caused by slow contact action and poor arc-extinguishing performance of switch [21,22].Figure 10 is a simplified equivalent circuit diagram for this process.In Figure 10, ) (t e is the AC power supply; t L and t C are used to simulate the lumped inductance and capacitance parameters of the transformer and other equipment.DS is the isolation switch, and its mathematical model can be established by a nonlinear time-varying resistance.r C is the lumped ground capacitance with initial voltage to simulate the residual charge; L1 and L2 are import and export short bus, respectively.All parameters of this example in Figure 10 were taken from Reference [23].The mathematical model of the VFTO simulation was established by using the telegraph equation for the short bus 1 L and 2 L .The short bus, L1, was uniformly divided into N1 = 20 sections, and N2 = 7 of the short bus L2 was selected.The following ordinary differential equations were obtained by the interpolation formulae [20]: As the simulation time-step of VFTO was nanosecond scaled, it could be considered that A was the constant matrix in a time-step and ) (t δ was the incentive source of VFTO. The four-step L2MF and the implicit trapezoidal integration method were used to solve Equation (16).The integration step of the implicit trapezoidal integration method was 1ns h = and a smaller integration step of L2MF was 0.1ns h = . The simulation results in this case are shown in Figure 11.
In this simulation test, the simulation time-  The mathematical model of the VFTO simulation was established by using the telegraph equation for the short bus L 1 and L 2 .The short bus, L 1 , was uniformly divided into N 1 = 20 sections, and N 2 = 7 of the short bus L 2 was selected.The following ordinary differential equations were obtained by the interpolation formulae [20]: As the simulation time-step of VFTO was nanosecond scaled, it could be considered that A was the constant matrix in a time-step and δ(t) was the incentive source of VFTO.
The four-step L 2 MF and the implicit trapezoidal integration method were used to solve Equation (16).The integration step of the implicit trapezoidal integration method was h = 1 ns and a smaller integration step of L 2 MF was h = 0.1 ns.The simulation results in this case are shown in Figure 11.
In this simulation test, the simulation time-step of VFTO was nanosecond scaled, which accurately simulated the change rule of the open-circuit voltage u(t) at the end of the line.Neither of these methods produced numerical oscillations, and the waveforms of VFTO for both the methods were almost the same.In this case, we simulate a situation when the disconnector of 550 kV gas insulated switchgear (GIS) acts on the no-load bus.During this process, very fast transient overvoltage (VFTO) is caused by slow contact action and poor arc-extinguishing performance of switch [21,22].Figure 10 is a simplified equivalent circuit diagram for this process.In Figure 10, ) (t e is the AC power supply; t L and t C are used to simulate the lumped inductance and capacitance parameters of the transformer and other equipment.DS is the isolation switch, and its mathematical model can be established by a nonlinear time-varying resistance.r C is the lumped ground capacitance with initial voltage to simulate the residual charge; L1 and L2 are import and export short bus, respectively.All parameters of this example in Figure 10 were taken from Reference [23].The mathematical model of the VFTO simulation was established by using the telegraph equation for the short bus 1 L and 2 L .The short bus, L1, was uniformly divided into N1 = 20 sections, and N2 = 7 of the short bus L2 was selected.The following ordinary differential equations were obtained by the interpolation formulae [20]: As the simulation time-step of VFTO was nanosecond scaled, it could be considered that A was the constant matrix in a time-step and ) (t δ was the incentive source of VFTO. The four-step L2MF and the implicit trapezoidal integration method were used to solve Equation (16).The integration step of the implicit trapezoidal integration method was 1ns h = and a smaller integration step of L2MF was 0.1ns h = . The simulation results in this case are shown in Figure 11.
In this simulation test, the simulation time-step of VFTO was nanosecond scaled, which accurately simulated the change rule of the open-circuit voltage ( ) u t at the end of the line.Neither of these methods produced numerical oscillations, and the waveforms of VFTO for both the methods were almost the same.

RC Series Circuit Simulation Calculation
This example is the series of a RC circuit, as shown in Figure 12. Figure 12a is a basic linear RC series test circuit.Figure 12b is the excitation voltage source waveform.The exact solution of the capacitor voltage in the circuit was given as follows: where ∆t = 0.1 ms; a time constant, τ = RC, τ is 0.1 ms; and u(t − ∆t) is a delayed unit-step signal.

RC Series Circuit Simulation Calculation
This example is the series of a RC circuit, as shown in Figure 12. Figure 12a is a basic linear RC series test circuit.Figure 12b is the excitation voltage source waveform.The exact solution of the capacitor voltage in the circuit was given as follows:   In this paper, the four-step L2MF and the backward Euler method were used to solve the voltages at both the ends of the capacitor.The calculated time-step of the two algorithms was both 0.1ms h = . The absolute error curve of the calculation results solved with the use of the two algorithms is shown in Figure 13.As shown in Figure 13, the calculation error of the four-step L2MF was smaller than that of the backward Euler method.This is because, although both of them were first-order, the absolute value of the local phase error of the four-step L2MF method was smaller than that of the backward Euler method.This example also indirectly verifies that the calculation effect of the ECDA combined with the four-step L2MF and the implicit trapezoidal integration method was better than that of the CDA.

Conclusions
Based on the traditional LMF, a class of low-order LMF was derived by the classic differential quadrature method.The L2MF method had L-stability, and the absolute value of the basic coefficient of the truncation error of four-step L2MF was smaller than that of the backward Euler method.Inspired by the CDA method, we proposed ECDA for electromagnetic transient simulation, which used the implicit trapezoidal integration method to perform integration operation under normal circumstances and switched the integration algorithm to the four-step L2MF method only when the network changed.The simulation results of the simple L2MF method and the implicit In this paper, the four-step L 2 MF and the backward Euler method were used to solve the voltages at both the ends of the capacitor.The calculated time-step of the two algorithms was both h = 0.1 ms.The absolute error curve of the calculation results solved with the use of the two algorithms is shown in Figure 13.As shown in Figure 13, the calculation error of the four-step L 2 MF was smaller than that of the backward Euler method.This is because, although both of them were first-order, the absolute value of the local phase error of the four-step L 2 MF method was smaller than that of the backward Euler method.This example also indirectly verifies that the calculation effect of the ECDA combined with the four-step L 2 MF and the implicit trapezoidal integration method was better than that of the CDA.

RC Series Circuit Simulation Calculation
This example is the series of a RC circuit, as shown in Figure 12. Figure 12a is a basic linear RC series test circuit.Figure 12b is the excitation voltage source waveform.The exact solution of the capacitor voltage in the circuit was given as follows:   In this paper, the four-step L2MF and the backward Euler method were used to solve the voltages at both the ends of the capacitor.The calculated time-step of the two algorithms was both 0.1ms h = . The absolute error curve of the calculation results solved with the use of the two algorithms is shown in Figure 13.As shown in Figure 13, the calculation error of the four-step L2MF was smaller than that of the backward Euler method.This is because, although both of them were first-order, the absolute value of the local phase error of the four-step L2MF method was smaller than that of the backward Euler method.This example also indirectly verifies that the calculation effect of the ECDA combined with the four-step L2MF and the implicit trapezoidal integration method was better than that of the CDA.

Conclusions
Based on the traditional LMF, a class of low-order LMF was derived by the classic differential quadrature method.The L2MF method had L-stability, and the absolute value of the basic coefficient of the truncation error of four-step L2MF was smaller than that of the backward Euler method.Inspired by the CDA method, we proposed ECDA for electromagnetic transient simulation, which used the implicit trapezoidal integration method to perform integration operation under normal circumstances and switched the integration algorithm to the four-step L2MF method only when the network changed.The simulation results of the simple L2MF method and the implicit

Conclusions
Based on the traditional LMF, a class of low-order LMF was derived by the classic differential quadrature method.The L 2 MF method had L-stability, and the absolute value of the basic coefficient of the truncation error of four-step L 2 MF was smaller than that of the backward Euler method.Inspired by the CDA method, we proposed ECDA for electromagnetic transient simulation, which used the implicit trapezoidal integration method to perform integration operation under normal circumstances and switched the integration algorithm to the four-step L 2 MF method only when the network changed.The simulation results of the simple L 2 MF method and the implicit trapezoidal integration method showed that the L 2 MF method could avoid the numerical oscillation completely, and the ECDA had a better simulation effect than the CDA.The next step will apply the ECDA in EMTP-type programs.

Figure 1 .Figure 2 .
Figure 1.Stability region of new linear multistep formulae (LMF) for even numbers of s.

Figure 1 .
Figure 1.Stability region of new linear multistep formulae (LMF) for even numbers of s. 0

Figure 1 .Figure 2 .
Figure 1.Stability region of new linear multistep formulae (LMF) for even numbers of s.

Figure 2 .
Figure 2. Stability region of new LMF for odd numbers of s.

Figure 3 .
Figure 3. Flow chart of electromagnetic transients simulation in power systems on a basis of extended critical damping adjustment method (ECDA).

Figure 3 .
Figure 3. Flow chart of electromagnetic transient simulation in power systems on a basis of extended critical damping adjustment method (ECDA).

Figure 4 .
Figure 4. Half-wave rectifier circuit.D represents the diode, which is an ideal switching device with voltage described as d(t).

Figure 5 .
Figure 5. Load voltage calculated by the three-step low-order linear multistep formulae (L2MF) method and the implicit trapezoidal integration method(TR).

Figure 4 .
Figure 4. Half-wave rectifier circuit.D represents the diode, which is an ideal switching device with voltage described as d(t).u RL (t) denotes a load voltage.

Figure 4 .
Figure 4. Half-wave rectifier circuit.D represents the diode, which is an ideal switching device with voltage described as d(t).

Figure 5 .
Figure 5. Load voltage calculated by the three-step low-order linear multistep formulae (L2MF) method and the implicit trapezoidal integration method(TR).

Figure 5 .
Figure 5. Load voltage calculated by the three-step low-order linear multistep formulae (L 2 MF) method and the implicit trapezoidal integration method(TR).

Figure 6 .
Figure 6.Schematic diagram of a single-phase high-voltage overhead transmission line without load being switched-in suddenly.

Figure 6 .
Figure 6.Schematic diagram of a single-phase high-voltage overhead transmission line without load being switched-in suddenly.

Figure 7 .
Figure 7.A piecewise model of a transmission line.

Figure 8 .Figure 9 .
Figure 8. Terminal voltage of the transmission line calculated by the four-step L2MF method and the implicit trapezoidal integration method.

Figure 7 .
Figure 7.A piecewise model of a transmission line.

Figure 7 .
Figure 7.A piecewise model of a transmission line.

Figure 8 .Figure 9 .
Figure 8. Terminal voltage of the transmission line calculated by the four-step L2MF method and the implicit trapezoidal integration method.

Figure 8 .
Figure 8. Terminal voltage of the transmission line calculated by the four-step L 2 MF method and the implicit trapezoidal integration method.

Figure 7 .
Figure 7.A piecewise model of a transmission line.

Figure 8 .Figure 9 .
Figure 8. Terminal voltage of the transmission line calculated by the four-step L2MF method and the implicit trapezoidal integration method.

Figure 11 .
Figure11.VFTO calculated by the four-step L2MF method and the TR method.

Figure 11 .
Figure 11.VFTO calculated by the four-step L2MF method and the TR method.Figure 11.VFTO calculated by the four-step L 2 MF method and the TR method.

Figure 11 .
Figure 11.VFTO calculated by the four-step L2MF method and the TR method.Figure 11.VFTO calculated by the four-step L 2 MF method and the TR method.

Figure 12 .
Figure 12.RC (resistance and capacitor) series test circuit and its excitation source.(a) RC series test circuit; and (b) excitation source waveform in the RC test circuit.

Figure 13 .
Figure 13.Computation errors of the L2MF and backward Euler algorithms as a function of time.

Figure 12 .
Figure 12.RC (resistance and capacitor) series test circuit and its excitation source.(a) RC series test circuit; and (b) excitation source waveform in the RC test circuit.

Figure 12 .
Figure 12.RC (resistance and capacitor) series test circuit and its excitation source.(a) RC series test circuit; and (b) excitation source waveform in the RC test circuit.

Figure 13 .
Figure 13.Computation errors of the L2MF and backward Euler algorithms as a function of time.

Figure 13 .
Figure 13.Computation errors of the L 2 MF and backward Euler algorithms as a function time.
1]; t 0 = t n is a time-starting point; t s = t n+1 is a time-end point; t i = t n + c i h, i ∈ (1, s) are internal points; and x n+i ≈ x( t i ) is the approximate value of a state variable at a time-internal point t i .