Transient Simulation of Distribution System with Converter-Based Distributed Generation for Stability Analysis Using Frequency Response Optimized Integrators

— This paper introduces a novel transient simulation scheme and its extension for stability analysis of distribution systems with converter-based distributed generation. Efficiency and accuracy of the simulation scheme with extension are investigated. It is demonstrated that millisecond-level step sizes are applicable to efficient electromagnetic transient simulation for stability analysis of unbalanced power systems while maintaining satisfactory accuracy. The paper also supports the extension of the simulation scheme by showing that the extension has only an insignificant impact on the accuracy but significantly enhances the efficiency.


INTRODUCTION
Traditionally, distribution systems were considered not to have stability issues due to their connection to strong voltage sources and their passivity property with limited amounts of distributed generation integration [1]. Studies on distribution systems focused on the steady state behavior via power flow computation and short-circuit computation [1]- [2]. In recent years, however, power systems world-wide are going through a major change with dramatically increasing penetration of distributed generation at the distribution level [1], [3]- [8]. This change significantly impacts the dynamics of distribution systems and potentially the overall power systems. As a result, stability analysis on distribution systems is drawing more and more research interest and effort [1], [3]- [6], [8].
Distributed generation is increasingly connected to power grids through power electronic converters and associated controllers, which form converter systems. Grid-connected converter systems can be classified as grid-feeding, gridforming or grid-supporting based on their operation modes [9]. Among them, grid-feeding converter systems are dominant [9]. They are studied in this paper as representatives of grid-connected converter systems integrated at the distribution level.
In the past, stability analysis was confined to bulk transmission systems [10]- [12]. Transient stability (TS) simulation, also known as electromechanical transient simulation, is a powerful tool in this type of study [10]- [13]. In TS simulation, voltages and currents in power system networks and electric machine stators are represented as quasi-stationary phasors; dynamics are considered only for electric machine rotors and controllers [10]- [12]. Computational efficiency of TS simulation is reasonably high with the typical millisecond-level step sizes [13]- [14], making it suitable for large-scale power systems. To enable TS simulation, several simplifying assumptions are necessary, including balanced three-phase conditions [3], [10]- [13]. Unfortunately, such a simplifying assumption does not hold for distribution systems, which are inherently unbalanced. As a result, TS simulation is not directly applicable to distribution systems.
Different from TS simulation, electromagnetic transient (EMT) simulation considers detailed models of power apparatus and is applicable to arbitrary single-phase or multiphase systems [15]- [16]. Ideally, EMT simulation should be applied to time domain simulation of distribution systems for stability analysis [3]. Unfortunately, computational inefficiency significantly prevents the applications of traditional EMT simulation to large-scale systems [3]. Such computational inefficiency is mainly due to the typical small step sizes ranging from microseconds to at most hundreds of microseconds adopted in traditional EMT simulation [13]- [14], [16].
In stability analysis, fast transients are not of interest [10]- [11]; lumped parameter models of power system networks and devices are used [10]- [12]; power electronic converters are represented with the averaged model [1], [10], [17]. These simplifying considerations can relax the constraint on step sizes. The millisecond-level step sizes conventionally adopted in TS simulation may seem appealing and plausible in EMT simulation if the same simplifying considerations are made. Nevertheless, the frequency-dependent error of the commonly used implicit trapezoidal method in traditional EMT simulation still limits the step sizes [18]- [20]. In particular, the implicit trapezoidal method introduces different degrees of error to different frequency components of the differential state variables. It is accurate for slow variants but its accuracy deteriorates as the frequency increases. To cater to the utility frequency (50 or 60 Hz) which is dominant in AC voltage and current waveforms, small step sizes have to be chosen to ensure sufficient accuracy. Otherwise, the simulation results are distorted; and the simulation process may even fail.
To overcome the drawback of the implicit trapezoidal method and enable large step sizes, several EMT-like simulation schemes are proposed in the literature [18]- [23]. Shifted frequency analysis (SFA) [20] calculates the frequency shifted analytic signals rather than the natural signals. After the frequency shift, the complex-valued signals become slow variants so the implicit trapezoidal method can be applied with large step sizes without introducing significant error. Frequency-adaptive simulation of transients (FAST) [21]- [22] switches between the frequency shifted analytic signals and the analytic signals adaptively to achieve multiscale simulation. References [18]- [19] modify the error frequency response of the implicit trapezoidal method to make it exactly accurate for the utility frequency and apply the modified numerical integrator to EMT simulation. Frequency response optimized integrators considering second order derivative are proposed and applied to EMT simulation in [23]. This paper applies the novel transient simulation scheme proposed in [23] and an extension of the simulation scheme proposed in [24] to a typical distribution system integrated with converter-based distributed generation and investigates the performance. A new transient simulation platform is introduced for stability analysis of distribution systems; potential applications of the simulation scheme with extension in this area are concurrently revealed. Millisecondlevel step sizes are shown to be efficient in EMT simulation without sacrificing accuracy. The paper also provides more evidence supporting the aforementioned extension, which enhances efficiency while leaving accuracy almost unchanged.
The remainder of this paper is organized as follows. Section II reviews numerical solution to ordinary differential equations (ODEs) and the novel transient simulation scheme with extension. Section III introduces the grid-feeding converter systems considered in this paper, including the single-phase and three-phase cases. It also states two computational schemes investigated in this paper, namely the original simulation scheme [23] and the extension [24]. The distribution system studied in this paper is described in Section IV. Simulation results are presented and discussed in Section V. Lastly, Section VI concludes the paper and points out some directions for future research.

II. NUMERICAL SOLUTION TO ODES AND THE NOVEL TRANSIENT SIMULATION SCHEME WITH EXTENSION
A. Numerical Solution to ODEs Consider a general ODE of the form ( , , ) where t is the time instant; x is the state variable; u is the input; f is a function depending on t, x and u. A single input is considered in (1) for clear presentation, but the same idea applies to multiple inputs. A single-step numerical integrator considering first order derivative may be applied to discretize (1) as where h is the step size; b 0 and b -1 are coefficients of the numerical integrator. The commonly used implicit trapezoidal method and backward Euler method are of this type of numerical integrators [25]. For the implicit trapezoidal method, b 0 = h/2 and b -1 = h/2 . For the backward Euler method, b 0 = h and b -1 = 0. Similarly, a single-step numerical integrator considering second order derivative may be applied to discretize (1) as where b 0 , b -1 , c 0 and c -1 are coefficients of the numerical integrator. The second order derivative of the state variable required in (3) is calculated as The coefficients in (3) can be chosen in such way that the resulting numerical integrator introduces no error at a specified nonzero angular frequency ω select despite the specific value of the step size. It is very common that the state variable of an ODE from the transient model of power systems is dominated by a nonzero frequency. For example, AC voltages and currents are typically state variables of ODEs describing capacitance and inductance; and they are dominated by the utility frequency. In these cases, ω select can be specified at the dominant frequency so as to significantly reduce the overall error introduced by the numerical integrator in the discretization process. The frequency response of the numerical integrator is thus optimized from the numerical error viewpoint. As a result, large step sizes are enabled to enhance efficiency in the simulation without sacrificing accuracy. Similarly, some numerical integrators are accurate for slow variants based on their frequency response of the error. They can equally be understood as frequency response optimized if they are applied to an ODE of which the state variable is slowly varying.
Reference [23] proposes and reviews four frequency response optimized integrators considering second order derivative, namely Integrators A-D. These numerical integrators are listed in Table I. Detailed investigations on the properties of the numerical integrators are performed in [26].
Their key features are summarized as follows: 1) Integrator A is highly accurate for state variables dominated by a nonzero frequency and slow variants. Provided the step size is smaller than one period defined by ω select , it behaves like an A-stable numerical integrator. Its ability to track fast decaying transients is weak. It should be replaced temporarily by Integrator B immediately after discontinuities. 2) Integrator B is accurate for state variables dominated by a nonzero frequency and slow variants. Some high frequency oscillations may cause numerical blow-up in transient simulation carried out with Integrator B. Its ability to track fast decaying transients is strong. 3) Integrator C is highly accurate for slow variants. It is Astable. Its ability to track fast decaying transients is weak. It should be replaced temporarily by Integrator D immediately after discontinuities. 4) Integrator D is relatively accurate for slow variants. It is Lstable. Its ability to track fast decaying transients is strong. Numerical integrators considering first order derivative can be inspected from the viewpoint of error frequency response as well. In fact, the implicit trapezoidal method and the backward Euler method are accurate for slow variants but their accuracy deteriorates as the frequency increases [18]- [20].

B. The Novel Transient Simulation Scheme with Extension
The transient model of power systems is typically a set of differential-algebraic equations (DAEs). The ODEs in the DAE set are discretized and converted to algebraic ones by numerical integrators to be solved numerically. To this end, the implicit trapezoidal method is extensively used in power system transient simulation [15]- [16]. Discontinuities have to be dealt with carefully to avoid numerical oscillation [27]- [28]. A common practice is to temporarily use the backward Euler method in place of the implicit trapezoidal method when a discontinuity takes place [20], [28].
The novel transient simulation scheme is proposed in [23] which is initially purely based on frequency response optimized integrators considering second order derivative. The possibility to extend the simulation scheme, by applying the commonly used numerical integrators considering first order derivative to the slow part of a power system model, is explored and justified in [24].
The novel transient simulation scheme with extension utilizes multiple numerical integrators to discretize the ODEs, rather than using a single numerical integrator for the entire system. The numerical integrators are selected in a case-bycase manner for individual ODEs according to the frequency spectrum of the state variables. Integrator A is used for the ODEs with state variables that have a nonzero dominant frequency component, which is typically the utility frequency or double the utility frequency. Integrator C is used for the ODEs with state variables that are relatively slow variants. The implicit trapezoidal method is used for the ODEs with state variables that vary slowly, if at all. The dominant frequency component of individual differential state variables can be determined based on engineering experience.
A fixed step size is used by the novel transient simulation scheme with extension during a simulation run, as is commonly adopted in power system transient simulation. The equations of the whole power system under study are simultaneously solved at each time step via iteration. This simultaneous solution process ensures consistency during the time-stepping to achieve high fidelity. When discontinuities are encountered, Integrator A, Integrator C and the implicit trapezoidal method are temporarily replaced by Integrator B, Integrator D and the backward Euler method respectively.

III. GRID-FEEDING CONVERTER SYSTEMS
The grid-feeding converter systems considered in this paper are briefly introduced in this section. The single-phase case is based on [29] while the three-phase case is based on [24]. More detailed descriptions and how the different numerical integrators are applied to different parts of the converter systems are given in those two references.

A. Overview
A high-level diagram of a grid-feeding converter system is shown in Fig. 1. For the three-phase case, Fig. 1 should be understood as the single-line diagram. The main objective of the converter system is to regulate its real and reactive power outputs. Externally given real and reactive power references are sent to the power controller. The power controller accordingly generates a control order for the remaining part of the converter system, which is usually expressed as a current phasor in the device reference frame. The overall behavior of the remaining part is similar to an equivalent controlled AC current source. In order for the equivalent controlled AC current source to feed the correct current waveforms into the external power grid, the cumulative phase angle of the terminal voltage is required, which provides information about the system reference frame.
The equivalent controlled AC current source is decomposed, in Fig. 2, into a filter, a voltage source converter (VSC) and a current controller. The current order phasor from the power controller is sent to the current controller, which accordingly generates a voltage order phasor. According to the voltage order phasor and the cumulative phase angle of the terminal voltage, the VSC generates the required voltage waveforms. The filter processes the output waveforms of the VSC so that the currents fed into the external power grid are close to sinusoids at the utility frequency.
In Figs. 1 and 2, a measurement and synchronization mechanism converts voltage and current waveforms into phasors for other controllers. It also calculates the real and reactive power outputs of the converter system. Moreover, the mechanism measures the cumulative phase angle of the terminal voltage, which is known as grid synchronization [9].

B. Filter, VSC, Power and Current Controllers
The filter and VSC compose the power part of the converter system. For a grid-feeding converter system, the filter is typically an inductor (L filter). In the single-phase case, the inductance of the L filter is denoted by L f . In the three-phase case, the inductance of each phase is assumed to be equal and decoupled; and L f is understood as the per-phase inductance.
In this paper, a VSC is simplified as a controlled AC voltage source. Its DC side is supposed to be strong enough to realize the functionalities and is not further modeled. This is a commonly accepted simplifying consideration in stability analysis [1], [17]. In the three-phase case, the output voltages of the VSC are balanced in positive sequence.
The power controller and the current controller mentioned in the previous subsection are shown in Figs. 3 and 4 respectively. They are applicable to both the single-phase and three-phase cases.

C. Measurement and Synchronization Mechanism
The in-phase and quadrature signals of a waveform are required in calculating the corresponding phasor. In the single-phase case, an improved second-order generalized integrator based quadrature signal generator (SOGI-QSG) [30] is used to generate the in-phase and quadrature signals of a given waveform x, which can be the terminal voltage or the current through the L filter. The improved SOGI-QSG is shown in Fig. 5. In the three-phase case, the in-phase and quadrature signals of x are calculated by Clarke's transformation. Here x is understood as vector-valued consisting of three-phase components. In particular 2 2 cos(0) cos( ) cos( ) Applying a phase shift factor based on the cumulative phase angle δ of the terminal voltage to the in-phase and quadrature signals of x, the corresponding phasor X in the device reference frame is obtained where X d and X q are the real part and imaginary part of X respectively. A phase-locked loop (PLL) is shown in Fig. 6. The imaginary part of the terminal voltage phasor is used as an input; ω syn is the synchronous angular frequency. The output of the PLL is the cumulative phase angle of the terminal voltage.
The real power output, reactive power output and terminal voltage phasor magnitude are first pre-calculated as The measured values are then obtained by passing these precalculated values to low-pass filters. The time constant of the low-pass filters for real power output, reactive power output and terminal voltage phasor magnitude measurement is denoted by T P , T Q and T V respectively.

D. Two Computational Schemes
The signals inside the power controller are slow variants. In this paper, two computational schemes are implemented. In Scheme 1 (the original simulation scheme proposed in [23]), Integrator C is used to discretize the ODEs of the power controller. In Scheme 2 (the extension proposed in [24]), the implicit trapezoidal method is used to discretize the ODEs of the power controller. The implementation of other parts of the studied power systems in both computational schemes is the same. Fig. 7 shows a diagram of the distribution system studied in this paper, which is based on the famous 13-bus test feeder [31]. A three-phase electricity source is connected to Bus 650 representing a conceptual transmission system. Three threephase grid-feeding converter systems are connected to the distribution system at Busses 675, 680 and 692 respectively through three-phase transformers. Two single-phase gridfeeding converter systems are connected to the distribution system at Busses 611 and 652 respectively through singlephase transformers.

IV. DISTRIBUTION SYSTEM
Detailed information about the distribution system can be found in [31]. The regulator is not implemented considering its slow response. The Phase A, B and C turn ratio of the transformer between Busses 650 and 632 is set to 1.0625, 1.05 and 1.075 respectively. When calculating per-unit values, the system power base is chosen as 1 MVA. In transient simulation, the static loads are modeled as constant impedance.
In this paper, the electricity source at Bus 650 is not assumed to be infinitely strong, so it is not an ideal AC voltage source. Instead, it is modeled as a synchronous generator. Dynamic parameters of the synchronous generator are listed in Table II. Synchronous generators have been implemented and validated with the novel transient simulation scheme in [23].
The parameters of the three-phase and single-phase transformers for converter system integration are listed in  Table III. The real and reactive power references for each three-phase converter system are 170 kW and 0 kVar respectively. Dynamic parameters of the three-phase converter systems are listed in Table IV. The real and reactive power references for each single-phase converter system are 50 kW and 0 kVar respectively. Dynamic parameters of the single-phase converter systems are listed in Table V. Note that uniform parameters for the transformers and converter systems are used in this paper only for elegance in presentation. The novel transient simulation scheme with extension does not make such unusual assumptions. In fact, each device can have its own parameters.

V. SIMULATION RESULTS
The novel transient simulation scheme with extension has been implemented with MATLAB. Its performance is investigated in this section via numerical case studies on the distribution system described in the previous section. Specifically, its accuracy and time consumption are compared to those obtained from an iterative EMT simulator [23] applied to the same test system. Accuracy of a simulation scheme given a step size can be determined by contrasting its results with the reference, which is generated by the iterative EMT simulator with a step size of 5 μs. The iterative EMT simulator is also implemented with MATLAB and has the same structure and execution as the novel transient simulation scheme with extension, except that all the ODEs of the studied system are discretized with the implicit trapezoidal method. The comparison on time consumption between the novel transient simulation scheme with extension and the iterative EMT simulator is fair because they adopt nearly the same implementations. Note that the iterative EMT simulator solves the whole system simultaneously without artificially inserted time delay [15]- [16] between different parts of the system, so consistency is achieved.
Step sizes much larger than those used in traditional EMT simulation are thus possible.
Simulation runs start from 0.0 s. Two scenarios are considered with the test system. In Scenario A, the rated power of the synchronous generator at Bus 650 is set to 500 MVA to mimic a strong transmission system. From 0.3 s to 0.5 s, a Phase-A-and-Phase-C-to-ground fault is applied at Bus 632 with a fault resistance of 0.001 p.u.. In Scenario B, the rated power of the synchronous generator is set to 5 MVA to mimic a weak transmission system. From 0.3 s to 0.5 s, a three-phase-to-ground fault is applied at the same bus with the same fault resistance as Scenario A. Later analysis in Subsection C will reveal that the system remains stable in Scenario A while it becomes unstable in Scenario B.
For simplicity in presentation, the synchronous generator rotor angle with respect to the synchronously rotating reference frame is referred to as rotor angle; the PLL output of a three-phase converter system with respect to the synchronously rotating reference frame is referred to as 3ph phase angle; the PLL output of a single-phase converter system with respect to the synchronously rotating reference frame is referred to as 1ph phase angle.

A. Qualitative Investigations
The generator speeds obtained from different simulation schemes with different step sizes in Scenario A are compared to the reference in Fig. 8. The comparison is performed for Scenario B in Fig. 9. For better visualization, only the dynamics between 1.2 s and 1.4 s are presented. It can be observed that Schemes 1 and 2 are much closer to the reference than the iterative EMT simulator with the same or one-half the step size. Their results basically overlap the reference. It is thus possible to use large step sizes with the novel transient simulation scheme and its extension to reduce the time consumption for simulation runs, while obtaining sufficiently accurate results.

B. Quantitative Investigations
In order to quantitatively investigate the accuracy of the simulation schemes, the following error metric is used in this paper. Suppose that x is a variable to be considered. The relative error regarding x of a simulation scheme given a step This paper studies the relative error of rotor angle, the averaged relative error of 3ph phase angles and the averaged relative error of 1ph phase angles. These three errors are referred to as rotor angle error, 3ph phase angle error and 1ph phase angle error respectively hereafter for simplicity. In particular, 3ph phase angle error is defined as 3 3 3 , 1 where δ 3ph,i is the 3ph phase angle of a three-phase converter system. The summation in (11) is made over all three of the three-phase converter systems. Similarly, 1ph phase angle error is defined as where δ 1ph,i is the 1ph phase angle of a single-phase converter system. The summation in (12) is made over all two of the single-phase converter systems. Rotor angle error, 3ph phase angle error, 1ph phase angle error and time consumption of different simulation schemes given different step sizes in Scenario A and B are compared in Table VI and VII respectively. Note that the smaller the error, the higher the accuracy. The duration of the simulation runs is 2.0 s. The iterative EMT simulator diverges with a step size of 4000 μs in both scenarios.
In Tables VI and VII, given the same step size, Schemes 1 and 2 typically consume more time to finish the simulation run than the iterative EMT simulator but less than twice as much. Nevertheless, this higher computational cost due to the numerical integrators considering second order derivative pays off with much higher accuracy. In other words, Schemes 1 and 2 are able to obtain similar or even more accurate results using significantly larger step sizes. According to Tables VI and VII, Schemes 1 and 2 can use a 4 times larger step size and still obtain better accuracy compared to the iterative EMT simulator in general. In some cases, they can even use an 8 times larger step size. Larger step sizes lead to fewer computations and consequently less time consumption for a simulation run. Computational efficiency is thus  achieved by Schemes 1 and 2. For example, Schemes 1 and 2 complete the simulation run in about 40 s with a 2000 μs step size in Scenario A; the iterative EMT simulator, however, takes roughly twice as much time with a 500 μs step size and still cannot match the accuracy of Schemes 1 and 2. Furthermore, Schemes 1 and 2 are more robust in that they still converge with a large step size of 4000 μs while the iterative EMT simulator diverges.
In Tables VI and VII, the difference in the errors of Schemes 1 and 2 given the same step size is subtle, which indicates that the extension impacts the accuracy of the novel transient simulation scheme only slightly. However, Scheme 2 always spends less time than Scheme 1. These observations support the extension of the novel transient simulation scheme, which indeed enhances efficiency without undermining accuracy.

C. Efficient EMT Simulation
The investigations in previous subsections have shown that the novel transient simulation scheme with extension is able to achieve efficiency and accuracy simultaneously with large step sizes. In this subsection, the extended novel transient simulation scheme (Scheme 2) is used with a step size of 1 ms to perform efficient EMT simulation of the distribution system under the two scenarios to determine if the system is stable or not in the post-fault stage. The synchronous generator rotor angle is used as the angle reference.
The 3ph and 1ph phase angles relative to the synchronous generator rotor angle are plotted in Figs. 10 and 11 respectively for Scenario A. It is shown that the system regains steady state after the fault is cleared. The system is thus stable. Note that there are 2 nd harmonics in the 3ph phase angles even in steady state. This is due to the negative sequence components in three-phase unbalanced systems.
The 3ph and 1ph phase angles relative to the synchronous generator rotor angle are plotted in Figs. 12 and 13 respectively for Scenario B. The system loses stability after the fault is cleared. It does not reach a steady-state operating point, as shown in the 3ph phase angles relative to the synchronous generator rotor angle. Also, the 1ph phase angles exhibit oscillations which do not exist in pre-fault steady state.

VI. CONCLUSION AND FUTURE WORK
In this paper, it is shown that millisecond-level step sizes, typically used in TS simulation for transmission system stability analysis, are indeed applicable to efficient EMT simulation with the novel transient simulation scheme and its extension while ensuring satisfactory accuracy. The simulation scheme with extension is thus suitable for stability    analysis of unbalanced power systems. In addition, more evidence is presented to support the extension. Although the test system in this paper is a small one, simulation results still show that the performance of the simulation scheme with extension is promising and deserves further investigation.
In the future, the extended novel transient simulation scheme may be applied to power apparatus and systems of different types and scales to further explore its performance.