Multi-Rate and Parallel Electromagnetic Transient Simulation Considering Nonlinear Characteristics of a Power System

Circuit (a) (b) Figure 6. The circuit used to increase node of network: (a) actual circuit; and (b) abstract circuit. 5.2. Analysis of Algorithm Efficiency 5.2.1. Influence of Subnet Partition on Efficiency In order to verify the influence of subnet numbers on efficiency, the power system is divided in different ways. Apparently, the scale of the AC/DC hybrid system shown in Figure 4 is rather small and unnecessary to divide it further into more than three subnets. Therefore, the circuit shown in Figure 6a is added to Figure 4 in order to increase the node number, and Figure 6b is the abstract expression of Figure 6a. The system after processing is shown in Figure 7.


Introduction
The electromagnetic transients program (EMTP) has been widely used for real-time simulation of power systems [1], verifying the control strategy [2], and detecting the operational status [3].However, one of the limitations of conventional EMTP-type simulators is their inability of performing efficient simulations.The scale of power systems is expanding [4], and the efficiency of electromagnetic transient simulations cannot meet the needs of scientific tests and theoretical research.
The concept of "diakoptics" was first introduced by Kron in the 1950s [5].Kron's diakoptics was aimed at the solution of large power system networks with increased computational speeds.Jose.R.Marti proposed the MATE concept in order to solve the network in partitioned form [6].The main philosophy of MATE is to split the solution of a large system into the solution of a number of smaller subsystems, plus the solution of the links joining the subsystems.Additionally, a multilevel MATE [7] and latency technique [8,9] are used in MATE, and these methods aim at simulating nonlinearities and solving the differential equations with different integration steps.
In this section, the fundamentals of the MATE-based parallel algorithm considering nonlinear characteristics of a power system is described, and the physical significance of the parallel algorithm is analyzed.

Fundamentals
Taking an abstract AC/DC hybrid power system as an example, the parallel algorithm considering nonlinear characteristics of the power system is described.The partition of the sample power system is shown in Figure 1.Subnet B is a DC grid, including rectifiers, inverters, DC lines; the components in branch set s-u are switches, so the admittance matrix of subnet B is changing.Both subnet A and C are AC grids, and supposing the component parameters of subnet A are unchanging, which is changing in branch set q-r of subnet C under different operation states (such as an on-load tap changer, OLTC), the admittance matrix of subnet A is unchanging, which is opposite to that for subnet C. It is worth noting that the node in Figure 1 does not specifically refer to a single node, but represents some node set, and we also suppose that "node set-node set" represents a branch set, such as l-y.In these branch sets, k-m, l-y, and n-x are used for network partitioning and are used to divide the network into different parts.s-u and q-r are the branch set with changing parameters, and they do not participate in the formation of admittance matrices.
The corresponding nodal equations of each subnet are as follows [6]:

+( )
( + ) where Gi, Vi, and hi denote the admittance matrix, node voltage, and accumulated current of subnet i (i denotes A, B, and C), separately.In particular, the elements of GB and GC don't include s-u and q-r parameters, respectively.Therefore, GA, GB, and GC are constant matrices.iα, iβ, and iγ denote the links' current between subnets.iρ and iη denote the links' current of s-u and q-r.pk, pl, pm, pn, py, px, pq, pr, and pu denote connectivity arrays, which reflect the relationship between the nodes and links' current.In particular, let pη denote pq + pr, pρ denote ps + pu.
As for s-u and q-r, the branch equations are: where zρ, zη, and eρ, eη denote branch links' Thevenin impedances and Thevenin voltages of s-u and qr, separately.From Equations ( 4) and ( 5), the dimension of iρ and iη are proportional to the number of changing branches.In general, the number of these changing branches in the AC subnet is small and is It is worth noting that the node in Figure 1 does not specifically refer to a single node, but represents some node set, and we also suppose that "node set-node set" represents a branch set, such as l-y.In these branch sets, k-m, l-y, and n-x are used for network partitioning and are used to divide the network into different parts.s-u and q-r are the branch set with changing parameters, and they do not participate in the formation of admittance matrices.
The corresponding nodal equations of each subnet are as follows [6]: where G i , V i , and h i denote the admittance matrix, node voltage, and accumulated current of subnet i (i denotes A, B, and C), separately.In particular, the elements of G B and G C don't include s-u and q-r parameters, respectively.Therefore, G A , G B , and G C are constant matrices.i α , i β , and i γ denote the links' current between subnets.i ρ and i η denote the links' current of s-u and q-r.p k , p l , p m , p n , p y , p x , p q , p r , and p u denote connectivity arrays, which reflect the relationship between the nodes and links' current.In particular, let p η denote p q + p r , p ρ denote p s + p u .As for s-u and q-r, the branch equations are: where z ρ , z η, and e ρ , e η denote branch links' Thevenin impedances and Thevenin voltages of s-u and q-r, separately.From Equations ( 4) and (5), the dimension of i ρ and i η are proportional to the number of changing branches.In general, the number of these changing branches in the AC subnet is small and is uncertain for the DC subnet.For example, changing branches of twelve pulse rectifier circuits is much less than modular multilevel converter (MMC).As for the proposed method, the number of changing branches should not be too large, or the calculation amount of the links' current would be extremely large.
As for k-m, l-y, and n-x, the branch equations are [6]: where z α , z β , z γ , and e α , e β , e γ denote the branch links' Thevenin impedances and Thevenin voltagse of k-m, l-y, and n-x, separately.By combining Equations ( 1)-( 8), Equation ( 9) can be obtained: i α , i β , i γ , i ρ , and i η could be obtained by solving Equation (9), and V A , V B , and V C could be solved via putting them into Equations (1)- (3).From Equation (9), solving the inversion of the admittance matrix is inevitable when obtaining the links' current.From above, G A , G B , and G C are constant matrices, and their inversion is needed to be solved once during the entire simulation.Additionally, both links' current of the changing and partition branches are solved together, which is different from the literature [7].In [7], changing branches are dealt with in their subnet calculation unit.

Physical Significance of the Parallel Algorithm
In the electromagnetic transient parallel algorithm based on MATE, the power system is divided into several subnets.The influences among different subnets are reflected through the links' current obtained by the Thevenin equivalent.The Thevenin equivalent regards some nodes as the research object, and acquires equivalent parameters in view of these nodes.Taking node k and l of subnet A as an example, the Thevenin voltage and impedance are described as follows: In Equation ( 9), p T k G −1 A h A and p T l G −1 A h A denote the Thevenin voltages of nodes k and l, respectively.p T k G −1 A p k and p T l G −1 A p l denote the Thevenin self-impedance of nodes k and l, respectively.p T k G −1 A p l and p T l G −1 A p k denote the Thevenin mutual-impedance between node k and l.

Multi-Rate Simulation in Parallel Algorithm
In this section, the concept of non-synchronized time and synchronized time is described, and the way that multi-rate simulation is applied to the above parallel algorithm is described.

Concept of Non-Synchronized and Synchronized Time
Supposing the simulation step of AC and DC subnets are ∆T and ∆t, and ∆T = λ∆t, in which λ is an integer, the relationship of ∆T and ∆t is shown in Figure 2.During the entire ∆T, only the DC In terms of network modeling, it is feasible to adopt different integral method.Common methods include the implicit trapezoidal integral method, Euler method, Gear integral method, and the Adams-Basfoss integral method.Considering the calculation amount and truncation error, this paper selects the implicit trapezoidal integral method to model the network.

Calculation Method at Non-Synchronized and Synchronized Time
Supposing the network is solved before 0 t , the following text describes the calculation method in non-synchronized and synchronized time among a T  .

Non-Synchronized Time
During non-synchronized time ti(i = 1,2,…,λ − 1), it is unnecessary and cannot calculate the Thevenin parameter of the AC subnet because of its simulation step In spite of this, the AC subnet's influence on the DC subnet exists during these periods, and it is obligatory to take this influence into consideration.In order to describe this influence, the paper uses Lagrange interpolate method to calculate Thevenin voltage of AC subnet.
As for AC subnet A or C, if its network solution of is solved at t0, history current source at the next synchronized time tλ of each element could be calculated according to the principle of electromagnetic transient simulation [16].Then, the current source vector hA(tλ) and hC(tλ) could be formed.Additionally, the Thevenin voltage, in view of the link nodes, can be obtained, and the expression of them are shown in Equations ( 10) and ( 11): where "~" denote node set k, l, x, y and branch set η. EA~(tλ) and EC~(tλ) denote the Thevenin voltage in view of "~" at subnets A and C, separately.
From above, the Thevenin voltage of subnets A and C at, and before, 0 t is known.Based on the Lagrange interpolation method, interpolation of the Thevenin voltage at any given moment can be obtained.Taking subnet A as an example, its interpolation Thevenin voltage is as follows: where K + 1 denotes the number of interpolated nodes.
() k t  is the interpolated basis function which takes t as an independent variable, and t  , tT   ,…, t k T   as interpolated nodes.The expression of () The Thevenin voltage of subnet A would be obtained when taking into Equation (12).The calculation for subnet C is completely identical to subnet A. In terms of network modeling, it is feasible to adopt different integral method.Common methods include the implicit trapezoidal integral method, Euler method, Gear integral method, and the Adams-Basfoss integral method.Considering the calculation amount and truncation error, this paper selects the implicit trapezoidal integral method to model the network.

Calculation Method at Non-Synchronized and Synchronized Time
Supposing the network is solved before t 0 , the following text describes the calculation method in non-synchronized and synchronized time among a ∆T.

Non-Synchronized Time
During non-synchronized time t i (i = 1, 2, . . ., λ − 1), it is unnecessary and cannot calculate the Thevenin parameter of the AC subnet because of its simulation step ∆T (∆T = λ∆t).In spite of this, the AC subnet's influence on the DC subnet exists during these periods, and it is obligatory to take this influence into consideration.In order to describe this influence, the paper uses Lagrange interpolate method to calculate Thevenin voltage of AC subnet.
As for AC subnet A or C, if its network solution of is solved at t 0 , history current source at the next synchronized time t λ of each element could be calculated according to the principle of electromagnetic transient simulation [16].Then, the current source vector h A (t λ ) and h C (t λ ) could be formed.Additionally, the Thevenin voltage, in view of the link nodes, can be obtained, and the expression of them are shown in Equations ( 10) and ( 11): where "~" denote node set k, l, x, y and branch set η. E A∼ (t λ ) and E C∼ (t λ ) denote the Thevenin voltage in view of "~" at subnets A and C, separately.
From above, the Thevenin voltage of subnets A and C at, and before, t 0 is known.Based on the Lagrange interpolation method, interpolation of the Thevenin voltage at any given moment can be obtained.Taking subnet A as an example, its interpolation Thevenin voltage is as follows: where K + 1 denotes the number of interpolated nodes.ω k (t) is the interpolated basis function which takes t as an independent variable, and t λ , t λ − ∆T, . . ., t λ − k∆T as interpolated nodes.The Thevenin voltage of subnet A would be obtained when taking t = t 1 , t 2 , . . ., t λ−1 into Equation (12).The calculation for subnet C is completely identical to subnet A. As for the DC subnet B, the history current source at t could be calculated using network solution at t − ∆t, and then they are used to form h B at t.Then, the Thevenin voltage of subnet B could be calculated, whose expression is similar to Equations ( 10) and (11).This paper will not repeat them here.Now, the Thevenin voltage of AC and DC at non-synchronized time have been calculated.Substituting them into the right of Equation ( 9), the links' current could be calculated.Then take i α , i γ , and i ρ into Equation ( 2), the subnet solution V B can be calculated.

Synchronized Time
Both AC and DC subnets need to be solved in synchronized time t λ .As for the AC subnet, its Thevenin voltage could be calculated using the last synchronized time; this part is already discussed in Section 3.2.However, as for the DC subnet, its network solution has been solved λ − 1 times during the last ∆T.Additionally, the topology of the DC subnet may have changed because of rapid switch events.Therefore, the influence of DC to AC could not be evaluated according to a certain non-synchronous simulation result of the DC subnet, and all the solution should be taken into account in order to improve simulation accuracy.Apparently, the traditional implicit trapezoidal integral method is not applicable because it only takes one calculation result to form the history source.Thus, the paper below takes an inductor in the DC subnet as an example, and illustrates its integral method.
The differential equation of inductor is shown as follows: We divide the interval [t − ∆T, t] into λ sub-intervals with length ∆t and applying the traditional implicit trapezoidal integral method to Equation ( 13) with integral step ∆t.Equation ( 14) can be obtained: Substituting all the equations in Equation ( 14) together, Equation (15) can be obtained: where R L is the equivalent resistance of inductor, h L (t) is the history current source of inductor at synchronized time, the expression of both are as follows: Now, h A , h B , and h C at synchronized time could be calculated.Inputting them into Equation ( 9), the links' current could be calculated, then they are substituted for corresponding items in Equations ( 1)-(3) to calculate V i (i denoting A, B, and C).

Simulation Process
In this section, the optimization of the algorithm is described.Then the simulation flowchart is shown to describe the simulation procedure.

Optimization of Simulation Process
The components in Equation ( 17) include the network solution in non-synchronized and synchronized time among the last ∆T.The summation part of the history term in Equation ( 17) could be updated at the end of each small time step ∆t, which avoids a large amount of calculation if all the addition tasks are executed when all the non-synchronized solutions have been obtained.
In addition, the process of calculating V A , V B and V C is optimized.The Equations ( 1)-( 3) are simply deformed and written in matrix form, and they are shown as follows: The first term G −1 i V i on the right side of the Equation (18), could be calculated before the links' current are obtained (i denote A, B, and C).Therefore, parallel calculations are used for solving G −1 i V i and the links' current.When the links' current are solved, we substitute them into the second term on the right side of the Equation (18), and V A , V B , and V C can be solved.

Simulation Process
Figure 3 is the simulation process chart of this algorithm.The upper dashed box is the preprocess, which is mainly used for obtaining the admittance matrix and its inversion, and only needs to be calculated once in the whole simulation process.The lower dashed box is a simulation process among ∆T, which is mainly used for pre-solution of the subnet network equation, acquiring the connection line current and the calculation of the node voltage.The calculation process of other ∆T is exactly the same as the lower dashed box.

Optimization of Simulation Process
The components in Equation ( 17) include the network solution in non-synchronized and synchronized time among the last T  .The summation part of the history term in Equation ( 17) could be updated at the end of each small time step t  , which avoids a large amount of calculation if all the addition tasks are executed when all the non-synchronized solutions have been obtained.In addition, the process of calculating VA, VB and VC is optimized.The Equations ( 1)-( 3) are simply deformed and written in matrix form, and they are shown as follows: The first term GV  on the right side of the Equation (18), could be calculated before the links' current are obtained (i denote A, B, and C).Therefore, parallel calculations are used for solving GV  and the links' current.When the links' current are solved, we substitute them into the second term on the right side of the Equation (18), and VA, VB, and VC can be solved.

Simulation Process
Figure 3 is the simulation process chart of this algorithm.The upper dashed box is the preprocess, which is mainly used for obtaining the admittance matrix and its inversion, and only needs to be calculated once in the whole simulation process.The lower dashed box is a simulation process among T  , which is mainly used for pre-solution of the subnet network equation, acquiring the connection line current and the calculation of the node voltage.The calculation process of other   In the upper dashed box, network equations of the DC and AC nets need to be written, and their admittance matrices are sent to the master side before simulation.Supposing the network has been solved before t 0 , the calculation process of the lower dashed box is illustrated as follows: (1) For the AC subnet, calculate E A∼ (t λ ) and E C∼ (t λ ) based on the calculation results of t 0 , which is shown in Equations ( 10) and (11), and interpolate to obtain their interpolation Thevenin voltage, which is shown in Equation ( 12). ( 2 (9) and send them back to the AC and DC nets.(12) During procedure (11), G −1 B h B (t λ − ∆t) is calculated at the same time.(13) V A , V B , and V C are calculated according to Equation (18).

Simulation Results
In this section, the accuracy and efficiency analysis of the proposed algorithm is analyzed.The proposed method is achieved by programing in Matrix Laboratory (MATLAB).

Accuracy Analysis of Algorithm
This paper uses an AC and DC hybrid system as the simulation example.The structure of this system is shown in Figure 4, and its parameters could be found in [17].It is necessary to illustrate that the OLTC in Figure 4 is added by the author.The fault mode is that there is a three-phase metal to ground fault at the AC bus of the inverting side when the time is 1 s and lasts 100 ms.
Energies 2018, 11, x 8 of 15 In the upper dashed box, network equations of the DC and AC nets need to be written, and their admittance matrices are sent to the master side before simulation.Supposing the network has been solved before t0, the calculation process of the lower dashed box is illustrated as follows: (1) For the AC subnet, calculate Et  based on the calculation results of t0, which is shown in Equations ( 10) and (11), and interpolate to obtain their interpolation Thevenin voltage, which is shown in Equation ( 12).
(3) For the DC subnet, calculate Et  to master side, calculate link currents according to Equation ( 9) and send them back to the AC and DC nets.(12) During procedure (11),  is calculated at the same time.

Simulation Results
In this section, the accuracy and efficiency analysis of the proposed algorithm is analyzed.The proposed method is achieved by programing in Matrix Laboratory(MATLAB).

Accuracy Analysis of Algorithm
This paper uses an AC and DC hybrid system as the simulation example.The structure of this system is shown in Figure 4, and its parameters could be found in [17].It is necessary to illustrate that the OLTC in Figure 4 is added by the author.The fault mode is that there is a three-phase metal to ground fault at the AC bus of the inverting side when the time is 1 s and lasts 100 ms.In order to verify the accuracy of the proposed algorithm, the simulation result is compared with the results generated by PSCAD.The simulation step of PSCAD is 10 μs.The network needs to divide In order to verify the accuracy of the proposed algorithm, the simulation result is compared with the results generated by PSCAD.The simulation step of PSCAD is 10 µs.The network needs to divide for this algorithm.This paper divides Figure 4 into three subnets, and each subnet is surrounded by a solid line.The simulation step for the AC sub-system is 50 µs, and is 10 µs for the DC sub-system.Additionally, the number of interpolated nodes is three and the integral method is the implicit trapezoidal integral method.The following paper would choose them as fundamental simulation settings if there are no additional instructions.Figure 5 shows the comparison of two simulation methods.In Figure 6, U ainv is a phase voltage of the AC bus on the inverter side, U dinv is the voltage of the DC bus on the inverter side, U arec is the a phase voltage of the AC bus on the rectifier side, U drec is the voltage of the DC bus on the rectifier side.From Figure 5, the simulation result of the proposed algorithm is in good agreement with that of PSCAD, and it proves the accuracy of the algorithm.

DC-Filter
Energies 2018, 11, x 9 of 15 for this algorithm.This paper divides Figure 4 into three subnets, and each subnet is surrounded by a solid line.The simulation step for the AC sub-system is 50 μs, and is 10 μs for the DC sub-system.Additionally, the number of interpolated nodes is three and the integral method is the implicit trapezoidal integral method.The following paper would choose them as fundamental simulation settings if there are no additional instructions.Figure 5 shows the comparison of two simulation methods.In Figure 6, Uainv is a phase voltage of the AC bus on the inverter side, Udinv is the voltage of the DC bus on the inverter side, Uarec is the a phase voltage of the AC bus on the rectifier side, Udrec is the voltage of the DC bus on the rectifier side.From Figure 5, the simulation result of the proposed algorithm is in good agreement with that of PSCAD, and it proves the accuracy of the algorithm.
(  As for the power system in Figure 7, the partition principle lies in the scale of divided subnet is as uniform as possible, because it can decrease the synchronization waiting process.The system in Figure 7 is divided into fourteen areas, and different partition schemes are shown in Table 1.This paper compares the simulation time using the proposed algorithm with traditional serial, single-rate MATE parallel and multi-rate parallel methods based on latency simulation methods.The simulation time is 1 s, and the simulation step of these methods are shown in Table 2

Influence of Subnet Partition on Efficiency
In order to verify the influence of subnet numbers on efficiency, the power system is divided in different ways.Apparently, the scale of the AC/DC hybrid system shown in Figure 4 is rather small and unnecessary to divide it further into more than three subnets.Therefore, the circuit shown in Figure 6a is added to Figure 4 in order to increase the node number, and Figure 6b is the abstract expression of Figure 6a.The system after processing is shown in Figure 7.As for the power system in Figure 7, the partition principle lies in the scale of divided subnet is as uniform as possible, because it can decrease the synchronization waiting process.The system in Figure 7 is divided into fourteen areas, and different partition schemes are shown in Table 1.This paper compares the simulation time using the proposed algorithm with traditional serial, single-rate MATE parallel and multi-rate parallel methods based on latency simulation methods.The simulation time is 1 s, and the simulation step of these methods are shown in Table 2.As for the power system in Figure 7, the partition principle lies in the scale of divided subnet is as uniform as possible, because it can decrease the synchronization waiting process.The system in Figure 7 is divided into fourteen areas, and different partition schemes are shown in Table 1.This paper compares the simulation time using the proposed algorithm with traditional serial, single-rate MATE parallel and multi-rate parallel methods based on latency simulation methods.The simulation time is 1 s, and the simulation step of these methods are shown in Table 2.

Influence of Subnet Scale on Efficiency
In order to verify the influence of subnet scale on efficiency, on the basis of Figure 4, this paper increases the nodes' number by a paralleling circuit shown in Figure 6b on the AC bus.The system is simulated using methods in Section 5.2.1.Additionally, the AC networks on the inverter and rectifier sides are both divided into two parts with similar scales.In other word, the number of subnets is five for each method.The number of paralleled circuits on the AC bus and the simulation time of these methods are shown in Table 3. From Table 3, the simulation time of methods 1 to 3 increases rapidly with the expansion of the network scale, but for the proposed method, the network scale has a little impact on simulation time.The reason lies in that the admittance matrix in methods 1 to 3 is changeable, so the inverse operation must be applied to it at each simulation step.However, as for the proposed method, the elements with changeable characteristics do not participate in the formation of the admittance matrix, and the influence of these elements to the network is reflected in

Influence of Subnet Scale on Efficiency
In order to verify the influence of subnet scale on efficiency, on the basis of Figure 4, this paper increases the nodes' number by a paralleling circuit shown in Figure 6b on the AC bus.The system is simulated using methods in Section 5.2.1.Additionally, the AC networks on the inverter and rectifier sides are both divided into two parts with similar scales.In other word, the number of subnets is five for each method.The number of paralleled circuits on the AC bus and the simulation time of these methods are shown in Table 3. From Table 3, the simulation time of methods 1 to 3 increases rapidly with the expansion of the network scale, but for the proposed method, the network scale has a little impact on simulation time.The reason lies in that the admittance matrix in methods 1 to 3 is changeable, so the inverse operation must be applied to it at each simulation step.However, as for the proposed method, the elements with changeable characteristics do not participate in the formation of the admittance matrix, and the influence of these elements to the network is reflected in the injection current, which means that the inversion operation is just needed once for the admittance matrix during the entire simulation.Additionally, the method of this paper considers the time constant differences of the AC subnet and the DC subnet, and adopts the multi-rate simulation mode, so the simulation efficiency is better than the other methods.The method multilevel MATE parallel simulation algorithm, which is explained in [7], together with the proposed method could make the admittance matrix constant through a special operation.The difference lies in that multilevel MATE guarantees the admittance matrix is constant by eliminating the link current in the subnet, and belongs to single rate simulation.However, the proposed method guarantees the admittance matrix is constant by regarding changing branches as the link current, and belongs to the multi-rate simulation.In order to compare the efficiency of these two methods, the simulation step of the AC/DC subnets are set to be same, and the optimization of the calculation process is canceled for the method of this paper.The simulation time of these two methods for the systems in Section 5.2.2 is shown in Table 4, in which t 1 and t 2 represent the simulation time of the proposed method and multilevel MATE, respectively.In these two methods, each network is divided into three subnets, which is also the same as the partition scheme in Section 5.2.2.The method of this paper deals with branches with changing parameters by assigning them to the link current calculation part.Multilevel MATE deals with branches with changing parameters by assigning them to the network equation solving part.From Table 4, the proposed method achieves better simulation efficiency at different network scales.The reason lies in the premise of eliminating the changing branches is the secondary subnet division as that for the multilevel MATE.However, in the process of the division, some branches with constant parameters are inevitably involved in the process of network secondary division, and the introduction of these branches will seriously affect the efficiency of eliminating branches.As for the method of this paper, it is just necessary to deal with the branches with changing parameters, and there is no need to divide the subnet again in this process, this illustrates that it does not need to involve constant branches into the calculation of the link current.Therefore, the proposed method is more efficient in the above examples.If the multi-rate simulation is applied in the above examples, the faster simulation speed will be achieved.

Number of Lagrange Interpolation Points
As for the AC network, this paper uses the Lagrange interpolation method to calculate Thevenin voltage of the AC subnet in non-synchronous time.Simulating the system of Figure 4 using different Lagrange interpolation points, and the average relative tolerance of the inverter DC bus voltage is shown in Table 5.The simulation settings are the same as Section 5.1 and the benchmark uses the PSCAD simulation results.From Table 5, when the interpolation point is three, the average relative tolerance is the smallest.The reason may lie in that when the interpolation point is two, the interpolation polynomial is the result of two-point linear interpolation, and the interpolation point is not enough, which results in the interpolation polynomial being unable to reflect the dynamic characteristics of the power grid.When the interpolation point is greater than three (not including three), the structure of the interpolation polynomial is mainly influenced by the moment before the upcoming simulation, so it cannot reflect the dynamic characteristics of the upcoming simulation moment, too.Additionally, it will increase the amount of calculation.
In fact, the average relative tolerance of the proposed method is influenced by the grid structure, number of subnets, partition method, and simulation step.Apparently, an "interpolation point of three" behaves the best in this network under this simulation condition.The choice of the interpolation point number considering various factors is the future research focus.

Network Modeling Method
Simulating the network shown in Figure 4 using different integration methods, the average relative tolerance and simulation time under different methods are shown in Table 6.The simulation settings are the same as Section 5.1 and the benchmark uses the Power Systems Computer Aided Design (PSCAD) simulation results.From Table 6, the simulation efficiency is influenced slightly by different simulation methods.The implicit trapezoidal integration method obtained the highest simulation accuracy, and the second-order Gear integral method is slightly worse than the implicit trapezoidal integration method.The reason lies in that the truncation error of the two methods is proportional to the cubic of the simulation step, and the coefficient is smaller.As for other integration methods, the truncation error is proportional to the square of the simulation step, so the simulation precision is lower.

Figure 6 .
Figure 6.The circuit used to increase node of network: (a) actual circuit; and (b) abstract circuit.

Figure 8 .
Figure 8. Simulation time of different methods under different partition schemes.

Figure 8 .
Figure 8. Simulation time of different methods under different partition schemes.

3 .
Influence of Algorithm Parameters be solved in t i (i = 1, 2, ..., λ − 1), which is called non-synchronized time.Both AC and DC subnets need to be solved in t λ , which is called synchronized time.It is necessary to point out that t 0 in Figure2is the previous synchronized time.needs to be solved in ti(i = 1,2,…,λ − 1), which is called non-synchronized time.Both AC and DC subnets need to be solved in t  , which is called synchronized time.It is necessary to point out that 0 t in Figure2is the previous synchronized time.
) For the DC subnet, let k = 1.(3) For the DC subnet, calculate E B∼ (t k ) based on the calculation results of t k − ∆t.(4) Send E A∼ (t k ), E B∼ (t k ), and E C∼ (t k ) to the master side, calculate the link currents according to Equation (9) and send them back to the DC net.(5) During procedure (4), G −1 B h B (t k − ∆t) is calculated at the same time.(6) V B is calculated according to the second line of Equation (18), and h B (t k ) is updated.(7) Let k = k + 1, and judge whether k < λ. (8) If k < λ, return to procedure (3).If not, go to procedure (9).(9) For the DC subnet, calculate E B∼ (t λ ) based on the calculation results of the whole ∆T.(10) During procedures (4) to (9), G −1 A h A (t 0 ) and G −1 C h C (t 0 ) are calculated at the same time, and E A∼ (t λ ), E C∼ (t λ ) could be obtained.(11) Send E A∼ (t λ ), E B∼ (t λ ) and E C∼ (t λ ) to master side, calculate link currents according to Equation

Table 3 .
Setting of simulation step.

Table 4 .
Comparison of simulation time.

Table 5 .
Average relative tolerance under different interpolation points.

Table 6 .
Influence of the integral method on the simulation.