Full electromagnetic transient simulation for large power systems

1 Introduction

The electromagnetic transient theory [1]based on nodal analysis,proposed first in the 20th century,is still the basic method of the most commonly used electromagnetic transient simulation programs.As the electromagnetic transient simulation is generally smallscale,it is traditionally employed in equipment-level simulation such as the overvoltage and control protection simulation.With the booming development of the DC transmission and new energy power generation in China,coupling between the power electronic system and the AC power grid is increasingly strong.Actual faults occurred in the China power grid in recent years,which suggest that non-fundamental transient processes,as well as the dynamic behavior of the control and protection system become main causes of mutual effects between AC and DC power systems,which has changed the traditional characteristics of power grids and brought challenges to the simulation for large power grids.Owing to the limitation of the electromechanical transient method in simulating power electronics and their mutual effect with the AC power grid [2],it is necessary to study the electromagnetic transient simulation for system-level applications.It is a considerable challenge to realize a full electromagnetic transient simulation for the entire grid,given the difficulty of ensuring simulation accuracy,stability,and efficiency,when the simulation scale is largely expanded.

There are many research documents on the electromagnetic transient simulation for large power systems.In terms of simulation accuracy and stability,the rotating machine model is a key issue.For the dq model,Park Transformation is performed when the machine interfaces with the network and non-linear equations related to the rotor angle are introduced.When explicit forecast or compensation interface methods are adopted [3,4],there will be some interface error,and the simulation is sensitive to the time step.For a relatively large system,the interface error can deteriorate calculation convergence [5].To tackle this problem,a phase domain(PD)machine model [6,7]needs to be employed.The PD model is based on an abc coordinate system wherein the machine and network form simultaneous equations without interface transformation,thus ensuring good numerical stability.The drawbacks of the PD model are that the inductance matrix of the machine is time-varying,meaning that LU decomposition is required at each time step,which lowers calculation efficiency.In [7],a virtual rotor winding was introduced to transform the time-varying matrix into a constant one,and the simulation error was found to be tolerable.This way,the problem is solved.

In terms of efficiency,parallel computing technology is indispensable in electromagnetic transient simulation.A traditional parallel algorithm includes a node tearing method and the branch partition method [8].Marti et al.[9]propose a multi-areas Thevnin equivalent(MATE)method for the network equation.Reference [10],based on MATE,proposes a a more efficient parallel algorithm by optimizing the calculation process and distributing resources in a reasonable manner.Reference [11]proposes the GENE method to enhance the efficiency of power systems with power electronic elements.The algorithm proposed in reference [12]shares similar ideas for the GENE method,which calculates AC and DC systems separately,avoiding frequent LU decomposition.Another approach to improve efficiency is to increase the simulation time step.According to the characteristics of power systems,real waveforms can be regarded as low-frequency signals modulated by a 50 Hz AC carrier.With the introduction of a frequency variable,only low-frequency signals are left after real waveforms are demodulated.Thus,a large step-size algorithm,i.e.,the frequency-adaptive algorithm,can be used. [13]and [14]describe the frequency-adaptive method in the frequency and time domains,respectively.The latter method features a clearer physical conception and relatively simple element model.In addition,the initialization of electromagnetic transient simulation is also an important issue that needs to be considered [15].Initialization is considerably more difficult for large power systems.

Based on the current level of electromagnetic transient simulation technology,this paper describes the main technical ideas of realizing full electromagnetic transient simulation for large power systems.The combination of small and large time steps algorithm is proposed and the grid is divided into part A and part B.In the small time step system,a two-level gird division structure is proposed,in which a multi-rate parallel structure is used among subsystems and an element-level parallel is used within subsystems.For power electronic elements,the method of integrating Norton Equivalence to the subsystem equation is proposed.Possible initialization ideas of large power systems are also described.According to the above ideas and methods,the top technical framework for full electromagnetic transient simulation of large power systems is designed.

2 Main technical ideas

2.1 Basic algorithm

Traditionally,a three-phase instantaneous value is used to represent electrical quantities.As the frequency scale is large,a small integration step is employed.For power systems with synchronous generators such as the main power supply,in most cases,the voltage and current are standard sine waves whose amplitude and phase position change at low frequency,far below the synchronous speed.Efficiency is low when a small time step is employed to simulate natural waves.For large power systems,if the amplitude and phase are extracted from the sine waves,the time step can be increased considerably,and thus,the calculation efficiency can be improved.The time domain method in [14]can be used as an idea for realizing the large time step simulation,enabling a smooth switch between the amplitude/phase and the natural waveforms.According to this idea,the entire power grid can be divided into different parts:(1)system A-power electronic elements and AC power grids with waveform distortion after disturbances use the small-time step algorithm.(2)System B-the remaining systems of the standard sine wave use the large time step algorithm.Systems A and B are interfaced through buses,as shown in Fig.1.The interface rules are:System B transmits three-phase instantaneous bus voltages to System A,while System A fits the three-phase instantaneous currents into amplitudes and phase angles,and then,they inject them to System B.The interface method is similar to that of the electromechanical-electromagnetic transient hybrid simulation [16].However,the difference is that System B is not an electromechanical transient simulation based on a fundamental sequence network; however,the electromagnetic transient simulation is based on a threephase instantaneous value.The interface is stricter in theory.

System B requires a low calculation amount,and therefore,the scale should be as large as possible to improve efficiency.First,the grid is far away from the power electronic element or fault location,and they can be selected as System B.When the system becomes stable after disturbances,the stable part of System A can be transformed into System B,and the interface buses will be changed accordingly.The operation above will be further elaborated in 2.4.

Fig.1 Schematic diagram of interface of system A and system B

2.2 Machine-network interface

A machine-network interface is one of the main problems restricting efficiency and stability of the electromagnetic transient simulation.The interface error of traditional dq models may affect simulation stability.In PD models,after certain transformations to machine windings,time-varying items in the inductance matrix are removed,and the matrix becomes constant,which solves the interface problem.Reference [7]proposed a specific method and the following is a brief description.

The equation set for the machine voltage and flux linkage is as follows:

Where,Ls,Lsr,Lrs,and Lr are the inductance matrices of stator/rotor windings; Lr is a constant value,while the remaining three are relevant to rotor position(i.e.,rotor angle).Substitute(2)into(1)to eliminate rotor current.With differencing of the implicit trapezoidal method,we get

Equation(3)shows the equivalent circuits of the machine,and they can be integrated into the network matrix easily for simultaneous calculation.The impedance matrix Zeq is:

Zeq consists of the constant part Z0 and the time-varying part relevant to rotor position.Unfortunately,owing to the nonlinear characteristics of the rotor angle,the iteration is required in the process of calculation.Add a virtual damper winding at the q axis of the rotor to realize Zq = Zd,which means that Zeq = Z0.In that way,the time-varying part of Zeq is removed,by avoiding matrix decomposition and iteration.The parameters of the resistance and the leakage reactance of the virtue winding can be set reasonably,whose decay time constant is very small,and thus,introducing an acceptable high frequency error outside of the frequency range interested .Therefore,the calculation efficiency and stability can be significantly improved.The PD model with constant inductance matrix should be the main model employed in electromagnetic transient simulation for large power systems.

2.3 Power electronics simulation

The core part of the power electronic element is the converter.With on-off state changes,the admittance matrix changes constantly as well.To minimize the influence,the basic idea is to separate the elements from external system,and therefore,they calculate them independently.

First,the system with transmission lines is decoupled.However,due to the time delay characteristics of transmission lines,power electronic elements are completely independent from external system.In this case,small time step is adopted for power electronic elements as their high switching frequency,while external system does not need such small time steps.Therefore,multi-rate simulation is necessary.In other words,a small time step is adopted for the power electronic elements and a large time step is employed for the external system.They exchange data at the time of each large time step.

Second,if the conditions for the transmission line separation are not met,the Norton equivalence of the power electronic elements is integrated into the external system,and the interface voltages are obtained.Then,the internal state of power electronic elements can be obtained.The converter bus or other PCC points can be selected as the interface,as shown in Fig.2.

Fig.2 Schematic of the combined calculation of power electronic equivalence and power grid

In this method,data exchange is necessary in each time step,and therefore,the time step of the external system is the same as that of the power electronic elements.This method is actually a special case of node tearing.The specific solving steps are as follows:

Step 1:Generate Norton equivalent circuits for the power electronic element:

Where,subscript d represents internal nodes,e represents interface nodes,id is the injection current of internal nodes,idh is the historical current of the internal nodes,and ieh is the historical current of the interface nodes.Eliminate ud,which is the voltage of internal nodes,with the Kron method.We then get:

Through equation(6)and Fig.2,the Norton equivalent parameters are:

Step 2:Equivalent circuits of power electronic elements and the external system are subject to the simultaneous solution so as to obtain ue ,which is the voltage of interface nodes.

Step 3:Calculate the internal state of power electronic elements.Based on equation(5)and as ue is known,

Step 4:Update the admittance matrix and the historical current items according to the on-off state and proceed to Step 1.

2.4 Parallel calculation

A parallel calculation is necessary to improve the efficiency of the electromagnetic transient simulation.For System B,with a large time step,the efficiency is no longer a main constraint.Therefore,parallel calculation is mainly used in System A.In terms of space,when the grid is divided using the transmission lines,subsystems are totally independent and an iteration is not needed,thus making this method the first option in full electromagnetic transient simulations for large power systems.However,there are requirements on matching the relation between the time step and the line length.For short lines not meeting time step requirements,the node tearing method can be used for grid division.The equivalent network composed of the set of tearing nodes is responsible for coordination and communication among subsystems.Therefore,a twolevel subsystem structure takes shape:level I is divided with a transmission line; level II is divided according to nodes within a single subsystem.Levels I and II nest with each other and parallel solving can be employed for level II subsystems.A different time step can be employed among different level I subsystems.Different time steps for the calculation can also be employed among different level II subsystems without power electronic elements.Level II subsystems with power electronic elements should share the same time step,solved according to 1.3.Fig.3 provides an example of grid division; the entire structure is as follows:System A is divided into four level I subsystems(S1-S4),with the grid division transmission line,b1,b2,and b3,among them.S1 is further divided into four level II subsystems(S1.1-S1.4).System B has just one subsystem S5,which interfaces with subsystems S3 and S4 of System A through nodes n1 and n2,respectively.The rules of interfacing are shown in Fig.1.When subsystem S3 of System A becomes stable,transform it and the connecting line b2 into System B to improve calculation efficiency.At this time,n2 and n3 become the interfaces of Systems A and B.n3 is the interface node between b2 and subsystem S2,while n1 is transformed into an internal interface of System B,as shown in Fig.4.

Fig.3 Example of two-level grid division

Fig.4 Schematic diagram of transformation of system A to system B

Grid division parallel is of relatively coarse granularity.Multi-core architecture of the high performance computer group can be employed for further parallelization for solving one subsystem.Take the node voltage method as an example.The network equation is

Where,Ib is the column vector of the branch currents,I is the column vector of injection currents,and A is the nodebranch incidence matrix.The branch represents the features of the elements; for the kth branch of Ib,it meets

where gk refers to branch admittance,u is the column vector of the node voltages,ik is the branch current,and ihk is the historical items of the branch current.Substitute(10)into(9),we get

The parallel strategy is as follows:

Step 1:Update matrix A and g according to network topology changes ➔ It is only necessary to find out the part that has changed; parallel calculation is not needed;

Step 2:Solve the historical item ihk of branch currents(k=1,…,b.Suppose there are b branches)➔ The branches are not coupled and parallel calculation can be conducted;

Step 3:Calculate the injection current of nodes,namely the item at the right of(11)➔ It is simple matrix multiplication and subtraction,and parallel calculation is not needed;

Step 4:Solve equation(11)➔ It is a solving process of linear simultaneous equation,and parallel solving can be conducted;

Step 5:Update the other state quantity of each element/branch ➔ Parallel solving can be conducted.

Step 6:Back to Step 1 and repeat the steps.

There are three calculation steps in which parallel solving may be used with no communication among the processes.The combination of parallel solving in the grid and parallel calculations of the grid division realize high degrees of parallelism,greatly increasing the calculation efficiency.In terms of hardware,a multicore computer group system and high-speed communication network are used to realize the parallel calculation strategy.

2.5 Initialization

The initialization of electromagnetic transient simulation is the assignment of an initial value of state variable for differential equations.Starting from the zero state is long process and the efficiency is low.Therefore,the load flow solution is used for initialization.Owing to the differences in the methods and the models between the electromechanical and the electromagnetic transient simulation,certain errors will occur when the power flow is used to derive the initial value of the state variable for differential equations of the power grid.Moreover,it is impossible to derive the internal state of power electronic elements with power flow.For initializing the electromagnetic transient for large power systems,the idea is to combine the power flow derivation with startup path optimization.

First,we conduct power flow derivation.Transform the power flow phasor into the three-phase node voltages and the branch currents,in order to clamp the port voltages of elements such as the generator,transformer,line,and power electronics.The body of elements will be initialized independently to remove the mutual effect among elements and shorten the time.

Second,start the system.Release the clamping of element ports step-by-step.Different release paths will cause different results,and even numerical oscillations; therefore,the startup path should be optimized.The key is to determine the startup sequence of generators and power electronic elements and ensure stability of the system.A trial-and-error mechanism can be introduced.After startup fails,determine the causes and exit.Start up the system again after path optimization till the system is started up successfully.The process will be finished automatically without manual intervention.

3 Framework design of full electromagnetic transient simulation

According to above ideas,full electromagnetic transient simulation technical framework for system-level application is designed with Unified Modeling Language(UML).The framework is as shown in Fig.5.In the figure,there is a main calculation process and nine categories of objects including system management,System B,equivalent System A,System A of level II subsystem,System A of level I subsystem,System A of power electronics,transmission lines,machines,and other elements.The object of system management controls the connection relation of the subsystems,and they are responsible for the top-level dispatch of the calculation process.The object of equivalent system A is an equivalent system formed by the level II subsystems based on node tearing,and it controls the connection relation of all level II subsystems(See 2.4 for details).The object of system A of power electronics refers to the power electronic system unqualified for the transmission line separation and solved with equivalence(See 2.3 for details).

There are three steps according to the calculation process.

Step 1:Read calculation data and implement preparatory work such as object establishment,parameter assignment and subsystem connection.

Step 2:Conduct initialization.As described in 2.5,the initialization process is divided into initial value derivation with power flow(corresponding to 2.2,2.3,and 2.4 of Fig.5)and system startup(corresponding to 2.5,2.6,2.7,and 2.8).The former is the process of deviating the initial value of the network busbar voltages and branch currents with the power flow results.For System A of power electronics,the value assignment is conducted only on the voltage of the PCC points.The latter is the process of starting up the system on the basis of clamping machine voltages and PCC voltages of power electronics.

Step 3:Calculate normal grid division.See 3.1-3.4 of Fig.5 for the calculation process of a time step.As the method is already introduced above,and further description is not provided here.Note that multi-rate simulation in the framework is naturally compatible,and different time steps can be adopted between System B and a different System A.Except the calculation between System A and System B,no iteration is needed in any other calculation.As the scope of System B is much larger than that of System A,it is expected that the total amount of calculation will be significantly lower than that of the existing electromagnetic transient program,while the calculation accuracy is still ensured.

Fig.5 Technical framework of full electromagnetic transient simulation for system-level application

4 Simulation analysis

The North China Power Grid is taken as an example for simulation analysis.The entire system contains 2,241 three-phase nodes and 365 generators.A Fault-free startup is conducted first.See Fig.6 for electromagnetic transient simulation curve.

Fig.6 Curves of electromagnetic simulation in fault-free startup

After the grid becomes stable,a three-phase short circuit fault occurs on some bus(See Fig.7 and Fig.8)for the voltage curve of the fault point and the power angle curve of the generators selected.Fig.7(a)and Fig.8(a)are the results of the electromagnetic transient simulation; Fig.7(b)and Fig.8(b)are the results of the electromechanical transient simulation with PSD-BPA [17].Fig.8 shows the power angle curve of an unstable generator and a stable one.It is clear that after the occurrence of the fault,the results of electromagnetic transient simulation and those of the mature electromechanical transient simulation are highly consistent,which verifies that the methods and technology framework of the electromagnetic transient simulation for large power systems mentioned in this paper is valid.

Fig.7 Voltage curves of the fault bus

Fig.8 Power angles curves of two generators

5 Conclusion

In this paper,key technologies of full electromagnetic transient simulation for large power systems are introduced; the general technical framework is proposed; and preliminary simulation comparison and analysis are conducted using a real example.This is the first step in realizing the goal,laying the foundation for the continuous improvement of technical details.Based on the framework proposed in this paper,the full electromagnetic transient simulation is universal and compatible.On the one hand,it can not only simulate the system behaviors when a large number of power electronics are interconnected with the AC grid,but they are also suitable for future DC grid development and simulation demand; on the other hand,it can be used for both system-level and equipment-level simulation.This paper focuses on software technologies; the hardware environment needed to realize full electromagnetic transient simulation is not discussed in detail.Under the general framework,adjustment and optimization of the algorithm as well as the model according to different hardware platforms are important and should be studied further.

Acknowledgements

This work was supported by key project of smart grid technology and equipment of national key research and development plan of China(2016YFB0900601).