Xin Ma, Xiao-Ping Zhang
©SHUTTERSTOCK.COM/VECTORMINE
Electromagnetic transient (EMT) simulation plays an important role is power system operation, control, and planning. Considering that the digital computer is unable to give a continuous history of transient phenomena, it is essential to solve the time response at discrete intervals reasonably. In the 1960s, Dommel first started an EMT solution by a nodal matrix method. To reduce solution time for nonlinear and time-varying components, Dommel proposed a compensation method and developed equivalent circuits for those components in 1971. To improve the calculation speed, Dommel introduced a new technique for solving EMTs with the implicit trapezoidal rule of integration in 1972. In 1986, Dommel published his EMT Program (EMTP) theory book, where all of the EMT component modeling and overall solution methods are detailed.
EMT simulation tools can be implemented either offline or in real time. The well-known Electromagnetic Transients Program (EMTP®) development started in 1966 by Dommel. In addition to EMTP, power system computer aided design (PSCAD) was first conceptualized by Manitoba Hydro International Ltd. in 1988 and has become a competitive alternative to network transient analyzers. Using a parallel processing structure, the real-time digital simulator (RTDS) was first introduced by RTDS Technologies Inc. in 1992 to simulate EMT in real time with a time step of 50–100 ${\mu}$s.
In these commercial programs, a library of components is available and a circuit model can be constructed by dragging, dropping, and connecting the appropriate model block on the drawing canvas. Although drawing a schematic diagram is simple, a typical student registered in undergraduate or postgraduate courses is subjected to learn exclusively from fundamental operations. Because of huge pressure and limited time, it is difficult to be familiar with a large quantity of materials in the curriculum. Nearly all EMT simulation tools start from modeling basic components. A poorly designed EMT solution can result in matrix singularity and nonconvergence, which will make tasks become unsolvable for the next time step. With the very strict requirement of accuracy, it is essential for engineers to be aware of EMT issues and how to avoid them.
The purpose of this article is to provide a practical and concise EMT modeling guideline to avoid common mistakes. The contributions can be summarized as follows:
The conclusions are supported for both offline simulators and real-time simulators. It is designed for rapid test of research and development ideas, especially for use in an advanced undergraduate and graduate level.
The article is organized as follows: For numerical stability, the “Integration Methods” section compares numerical integration methods. For discrete-time modeling, the “EMT Component Models” section provides a step-by-step guide to develop ordinary differential equations for different EMT components. To speed up EMT computation, the “Network Solution” section provides efficient methods. For verifying these models and methods, the “Case Studies” section gives a four-machine 11-bus case study. The last section summarizes conclusions.
Basically for EMT studies, all system components are normally represented by difference equations. EMT simulators are expected to provide fast and accurate solutions for power grids with different components. Digital techniques are unable to give a continuous calculation of the transient phenomena, but rather a sequence of discrete intervals. A continuous period is required to be discretized from the current time step to the next time step until the maximum time is reached. To meet the simulation requirement, a simulation time step needs to be extremely small for discretization. For solving ordinary differential equations, common numerical methods include back Euler, forward Euler, and the trapezoidal rule. Figure 1 details a comparison among back Euler, forward Euler, and the trapezoidal rule:
Figure 1. (a)–(c) Integration methods comparison.
By comparison, the trapezoidal rule is simple, numerically stable, and accurate enough to integrate the ordinary differential equations in EMT. Dommel discussed that discretization using forward Euler or backward Euler causes truncation errors and can lead to numerical instability, in the EMTP Theory Book in 1986. Therefore, the trapezoidal method is the most frequently used in EMT simulators, for instance EMTP. In the “Network Solution” section, how this is implemented in EMT simulation will be elaborated.
Based on arithmetic operations, EMT component models can be divided into linear and nonlinear components:
By applying the trapezoidal rule, a linear component can be converted into a Norton equivalent circuit, which has an equivalent conductance in parallel with a known history current source. This is known as the companion circuit model.
Linear lumped components are resistances, inductances, and capacitances (RLC) branches.
To derive an accurate model for an arbitrary RLC branch between node k and node m, we can start from a single inductor of the circuit. Any current change through an inductor creates a changing flux and induces voltage. Figure 2 details the derived procedure and equivalent circuit for the inductor.
Figure 2. (a) and (b) EMT solution for RLC branch.
According to Faraday’s law of induction, inductance between node k and m starts from an ac steady-state solution. After the trapezoidal rule is applied to the integral of the inductor, the exact differential equation is replaced by the approximate discrete equation, and can be rewritten as an equivalent conductance in parallel with a history term, already known from the solution at the previous solution point at time ${t}{-}\Delta{t}$. Therefore, once currents and voltages are obtained at time ${t}{-}\Delta{t}$, the history terms can be updated for each individual inductor at time t.
To simplify modeling, it is essential to get general representation of an RLC branch. Similar to a single inductor, the EMT solution of an arbitrary RLC branch can be obtained if the trapezoidal rule is applied to the integral. Figure 2 details the mathematical models and equivalent circuits of an RLC branch. In spite of series or parallel connections, any RLC branch can be represented by a conductance GRLC in parallel with a history term IRLC. For the same RLC branch between node k and node m, parameters GRLC, HRLC, and JRLC remain constant for the whole calculation.
Table 1 summarizes the detailed parameters and physical representation for common RLC branches. For example, The RC branch can be treated as a two-terminal network by combining unessential buses; therefore, the proposed mathematical models and parameters can effectively reduce computational burden.
Table 1. Comparison between different RLC branches.
The transmission lines modeling can be classified into a lumped parameter model and distributed parameter model. Based on the traveling wave model, the distributed parameter model implements a three-phase distributed line model where the RLC are uniformly distributed along the line. Contrary to the distributed parameters model, a three-phase ${\pi}$-section line lumped parameter model implements a balanced three-phase transmission line model with lumped parameters in a ${\pi}$-section circuit. For transient simulation, the traveling wave model is faster and more accurate than the lumped parameter model, as shown in Dommel’s EMTP Theory Book. In 1995, Zhang proposed the theoretic analysis and guidelines on the selection of a different number of equivalent ${\pi}$ sections of the transmission line model according to the frequency range of interest. A lumped model can be well represented by combining RLC branches. Table 2 details the comparison between the lumped parameter model and distributed parameter model.
Table 2. Comparison between lumped parameter line model and distributed parameter traveling wave model.
To get an accurate model of three-phase distributed transmission lines with loss, we start from single-phase lossless transmission lines. Figure 3 details that an EMT solution for transmission lines needs to go through four steps. At the first step, the surge impedance and travel time can be obtained based on traveling wave theory. At the second step, the distributed parameter model of lossless transmission line can be represented as an equivalent resistance and history term, similar to the RLC branch. The difference is that the history term is not instantly integrable because of travel time ${\tau}$. At the third step, a transmission line with loss is modeled with lumped R/4 in the terminals of transmission lines and R/2 in the middle of the transmission line, in which great accuracy were proved in experiments in Dommel’s article, published in 1969. At the fourth step, mutual three phases are converted to uncoupled modal domains to avoid an off-diagonal conductance matrix. By removing mutual couplings, each modal domain equation can be solved independently using surge impedance and travel time.
Figure 3. (a) and (b) EMT solution for a transmission line.
Contrary to linear components, nonlinear components cannot be well represented by a simple Norton circuit. This is because additional nonlinear calculations, such as sine wave function and limiter, are involved in modeling nonlinear components. Typical nonlinear components include the synchronous machine and control system.
Synchronous machines play a significant role in controlling power system stability. A universal model is built up to represent synchronous machines, induction machines, and dc machines. Based on electric machine theory, synchronous machines consist of magnetic rotors and stators. It is always challenging to model electrical and mechanical parts simultaneously in a large system.
The key solution is to integrate differential equations for a whole synchronous machine. Because of the coupling relationship, the electrical part and mechanical part are unable to be solved simultaneously. To overcome this, we start from predicting the mechanical part because electrical changes are much faster. Figure 4 presents the detailed calculation of synchronous machines. The first step is to predict mechanical angular speed and angle by linearization so as to eliminate numerical errors. The second step is to convert physical quantities from abc reference of frame to dq0 reference of frame using Park’s transformation. Based on the trapezoidal rule, the third step is to update currents in armature and field windings according to a voltage–current–flux relationship. The fifth step is to update electrical torque. By back substitution, the sixth step is to solve the question of the exact mechanical part model, and update history values and then go back to the first step until the maximum solution time duration is reached.
Figure 4. EMT solution for synchronous machine.
The modeling of a control system, such as an excitation system, is the combination of basic arithmetic operations, limiters, and a proportional-integral-derivative (PID) controller. To represent a wide variety of excitation systems, IEEE has standardized 12 models for the excitation system currently in use, such as ac1A, dc1A, and ac4A, etc. To obtain an accurate model of an excitation system, we start with an ac4A excitation system to provide field voltage to control cd and terminal voltage. Figure 5 presents the mathematical model and block diagram of the ac4A system. At the first step, the output can be obtained by adding the referenced stator terminal voltage, stabilization voltage, and compensated voltage. At the second step, the limiter block ensures that the output will not exceed the upper bound and lower bound. At the third step, transient gain reduction is implemented by a PID control block with transient gain reduction lead and lag time constants, which can be discretized as a history term and current term. At the fourth step, the main regulator is implemented by a nonwindup limiter. With the nonwindup limits, the output comes off the limits as soon as input re-enters the range within limits. To overcome the coupling relationship between input and output, predicted value is linearized by history values. Applying the trapezoidal rule, nonwindup limiter can be partitioned as linearization, PID control block, and limiter in the discrete time domain.
Figure 5. (a) and (b) EMT solution for ac4A excitation system.
In addition to excitation systems, prime movers and governor systems are to provide mechanical torques for synchronous machines. Steam turbines, as one of the most common prime movers, can convert high temperature and high pressure into rotating energy.
Figure 6 shows the block diagram of a steam prime mover and a governor system. The first step is adder. The second step is multiplier. The third step is subtracter. At the third step, the load setup point is obtained by a subtracter. Based on the trapezoidal rule, a low-pass filter is implemented at the fourth and fifth steps. Therefore, steam governor and prime mover can be discretized step by step.
Figure 6. (a) and (b) EMT solution for steam governor and prime mover system.
Node voltage solution is able to straightforwardly provide EMT solutions for a network only consisting of linear components. This is because all linear components, such as the RLC branch, can be represented with an equivalent resistance in parallel with a history term. Equivalent resistance is involved in the admittance matrix, while the history term is considered as an injected current source. For a known node m, such as ground and voltage source, row m in the Y matrix can be set zero except set 1 for diagonal components. By Gaussian elimination, the remaining unknown node voltage can be obtained by ${U} = {\text{inv}}{(}{Y}{)}\,{\ast}\,{I}$. Therefore, node voltage and branch current at time t can be updated consistently using a node voltage solution.
When nonlinear components are involved, the compensation method is required to deal with nonlinearity. The key of the compensation method lies in isolating the nonlinear part from the linear parts. In the compensation method, nonlinear components are treated as current injections, which are superimposed on the rest of the linear network. This treatment was initiated by Dommel and has been detailed in his 1986 EMTP Theory Book.
Figure 7 details applying a compensation method to the simultaneous solution of nonlinear and linear equations. At the first step, nonlinear components are excluded from the original network to get a Norton circuit of the rest of the open-circuit linear network. At the second step, open-circuit voltage can be obtained by setting the current injection of the nonlinear component as zero. At the third step, equivalent conductance can be calculated by replacing the current vector with a new vector whose components are all zero, except for +1 in row k and −1 in row m. At the third step, closed-circuit terminal voltage can be obtained after nonlinear components calculation. At the fourth step, the closed-circuit voltage superimposes the response into the current vector to complete the node voltage solution. Therefore, the compensation method can effectively get the simultaneous solution of both linear equations and nonlinear equations, simplifying computation burden.
Figure 7. (a) and (b) Application of compensation method.
Power system transients are characterized by containing very high-frequency components. In 1999, Henschel proposed the frequency-shift method in his Ph.D. thesis. Basically, If the system returns to steady state, the signal spectrum becomes around the system frequency, which is called a bandpass signal. Frequency shift is to analyze a time function in the shifted time domain.
The basic procedures for frequency shift as shown by Henschel in his Ph.D. work in 1999 include: At the first step, a power system signal and dynamic phasor is represented with in-phase and quadrature. At the second step, an analytic signal is defined based on a Hilbert transform. At the third step, an analytical signal is derived using Fourier transform and system frequency is shifted by angular frequency ${-}{\omega}_{0}$.
This method allows a much larger simulation step to be applied to improve simulation efficiency without loss of accuracy.
To explore the detailed behaviors of power systems, it is necessary to construct circuits and observe simulation results. Depending on execution time and the simulation step, EMT solutions consist of offline simulation methods and real-time simulation methods. Table 3 illustrates the difference between real-time simulation and offline simulation.
Table 3. Comparison between real-time and offline simulation.
All of the case studies are implemented using a single Virtex6 ML605 FPGA board. This board has been equipped with 768 digital signal processors (DSPs), 416 RAMs, 301,440 registers, and 150,720 look-up tables (LUTs). A 100-MHz system clock is selected for global synchronization. All of the floating-point calculations are using intellectual property cores (IP Cores). The simulation step for EMT programs is 50 ${\mu}$s.
Real-time performance is validated in three aspects, including resource utilization, accuracy, and timing.
In summary, to ensure real-time implementation, resource utilization and timing constraints should be operated. In other words, the bitstream file, executable program is unable to be generated when overused resources or failed timing constraints occur.
To verify the effectiveness of EMT mathematical models, a four-machine 11-bus system given in Figure 8(a) is implemented on a single ML605 FPGA board. A single-phase grounding fault and three-phase fault are applied at bus 7 from 1 s to 1.1 s, respectively. Figure 8 presents simulation results of electric torque, mechanical torque, and current and time schedule. Figure 8(c)–(e) presents simulation results with a single-phase grounding fault, and Figure 8(f)–(h) presents simulation results with a three-phase grounding fault.
Figure 8. (a)–(g) Four-machine 11-bus topology and simulation results.
It can be seen from Figure 8 that the relative error of electric torque, currents, and mechanical torque can be maintained less than 5%. With the trapezoidal rule, this relative error can be maintained in a limited range for both single-phase fault and three-phase fault. In particular, the d axis current and electric torque simulation results illustrate the proposed models can well represent dynamics of nonlinear calculations, for instance Park’s transform. It can be demonstrated that the proposed mathematical models, such as synchronous machines and RLC branches, can well represent the dynamic behaviors under both symmetric and asymmetric grounding faults.
Figure 8(b) shows the total calculation only costs 25.4 ${\mu}$s at a simulation time step 50 ${\mu}$s. Therefore, further calculation blocks are still allowed at available time slots. Seen from Table 4, it is observed that only 35% of registers and 87% of LUTs are utilized for the four-machine 11-bus system. Therefore, a larger system is still allowed on the same FPGA board, which shows high portability and flexibility.
Table 4. Resource utilization of 11-bus four-machine system.
To trigger higher frequency transients, the standard deviation is used to measure the variation of electric torque. In the referenced model, standard deviations of electric torque are 0.3318 and 0.0892, respectively, where the impact of the single-grounding fault is much more significant than that of the three-phase grounding fault. This is solid evidence as to why a single-phase grounding fault is more preferable for EMT simulation verification than a three-phase grounding fault.
This article has presented EMT models and solutions of power systems with components, such as RLC branches and synchronous machines. Theoretical analysis and graphical comparison can help beginners develop their fast familiarity with EMT models and solutions. Based on theoretical analysis and simulation results, the following conclusions can be given:
H. W. Dommel, EMTP Theory Book. Portland, OR, USA: Bonneville Power Admin., Aug. 1986.
H. W. Dommel, “Digital computer solution of electromagnetic transients in single-and multiphase networks,” IEEE Trans. Power App. Syst., vol. PAS-88, no. 4, pp. 388–399, Apr. 1969, doi: 10.1109/TPAS.1969.292459.
H. W. Dommel, “Nonlinear and time-varying elements in digital simulation of electromagnetic transients,” IEEE Trans. Power App. Syst., vol. PAS-90, no. 6, pp. 2561–2567, Nov. 1971, doi: 10.1109/TPAS.1971.292905.
H. W. Dommel and N. Sato, “Fast transient stability solutions,” IEEE Trans. Power App. Syst., vol. PAS-91, no. 4, pp. 1643–1650, Jul. 1972, doi: 10.1109/TPAS.1972.293341.
M. Gole, O. B. Nayak, T. S. Sidhu, and M. S. Sachdev, “A graphical electromagnetic simulation laboratory for power systems engineering programs,” IEEE Trans. Power App. Syst., vol. 11, no. 2, pp. 599–606, May 1996, doi: 10.1109/59.496082.
RTDS Technology Inc., Winnipeg, MB, Canada. Real-Time Simulation. (2014). [Online] . Available: http://www.rtds.com
X.-P. Zhang and H. Chen, “Analysis and selection of transmission line models used in power system transient simulations,” Int. J. Electr. Power Energy Syst., vol. 17, no. 4, pp. 239–246, Aug. 1995, doi: 10.1016/0142-0615(95)00039-S.
S. Henschel, “Analysis of electromagnetic and electromechanical power system transients with dynamic phasors,” Ph.D. dissertation, Dept. Elect. Eng., Univ. British Columbia, Vancouver, BC, Canada, 1999.
Xin Ma (x.ma.1@bham.ac.uk) is with the Department of Electronic, Electrical, and Systems Engineering, University of Birmingham, B15 2TT Birmingham, U.K.
Xiao-Ping Zhang (x.p.zhang@bham.ac.uk) is with the University of Birmingham, B15 2TT Birmingham, U.K.; the Birmingham Energy Institute, Birmingham B25 8FJ, U.K.; and the Birmingham Center for Energy Storage, University of Birmingham, B15 2TT Birmingham, U.K.
Digital Object Identifier 10.1109/MELE.2023.3320485
2325-5897/23©2023IEEE