Equivalent model of multi-type distributed generators under faults with fast-iterative calculation method based on improved PSO algorithm

There are various types of distributed generators (DGs) with different grid integration strategies. The transient characteristics of the fault currents provided by the DGs are different to those of conventional synchronous generators. In this paper, a distribution network with multi-type DGs is investigated, including consideration of DG low-voltage ride through (LVRT). The fault current characteristics of two typical DGs, i.e. an inverter-interfaced distributed generator (IIDG) and a doubly-fed induction generator (DFIG), are analyzed, considering the specific operation modes. Based on analysis of the fault characteristics, an equivalent model of the multi-type DGs under symmetrical/asymmetrical fault conditions is established. A fast-iterative fault calculation method for enhancing the calculation efficiency while avoiding local convergence is then proposed using an improved particle swarm optimization (PSO) algorithm. A simulation system of the distribution network with multi-type DGs is established in PSCAD/EMTDC. The simulation results validate the high accuracy and calculation efficiency of the proposed calculation method of the fault components. This can assist in the settings of the protection threshold.


Introduction
With the rapid development of the modern distribution network and the access of distributed generation, the structure of distribution network is becoming increasingly complex [1,2]. New challenges have been raised in the relay protection of a distribution network because of the integration of various types of distributed generators (DGs) and increased penetration of DGs [3,4]. To ensure safe and stable operation of the modern distribution network and prevent chain losses of DGs under fault conditions, several requirements have been formulated requiring DGs to have low-voltage ride through (LVRT) capability [5,6]. During the LVRT process, DG fault current characteristics are significantly different to those of conventional synchronous generators [7]. This makes it difficult to calculate the fault components and satisfy relay protection requirements [8][9][10]. Hence, it is necessary to investigate the modelling techniques of a DG integrated distribution network and find effective calculation methods of the fault components under fault conditions.
In general, DGs can be divided into two types according to their forms of integration: asynchronous (e.g., doubly-fed induction generator, DFIG) and inverter-based (e.g., photovoltaic, PV). The fault currents provided by the DGs are closely related to the LVRT strategies, forms of integration and the protection strategies. Research has been conducted on the analysis of fault current characteristics and establishment of equivalent models of the DGs. For the inverter-based DG, also known as inverterinterfaced distributed generator (IIDG), the nonlinear relationship between the fault current and terminal voltage is analyzed in [11], where the IIDG is represented by an equivalent impedance in series with a constant voltage. However, the corresponding calculation method is not investigated. In [12], a calculation method of a distribution network under faults is investigated and both sub-transient and steady state equivalent models are established. However, the fault current characteristics are not analyzed. In [13], the fault current characteristics of an IIDG in microgrids are studied, though the control strategy in the islanded operation mode in microgrids is different to that in a distribution network, resulting in different fault current characteristics. For asynchronous type DGs, e.g., DFIG, crowbar protection is generally considered and the back-to-back converter is rapidly blocked during most LVRT processes [14,15]. The fault characteristics of the DFIG under continuous excitation of the rotor-side inverter have also been considered, such that the DFIG can maintain grid connection and provide voltage/reactive support under specific fault conditions [16,17]. A transient equivalent model of the DFIG under continuous excitation is explored in [16] and a calculation method of the fault current is proposed with the assumption of no sudden change of the flux linkage. However, the impact of asymmetrical faults is not considered. In [17], an equivalent model of the DFIG with continuous excitation under faults is analyzed, but the reactive support of the DFIG is not investigated. In [18], the dynamic responses of the DFIG active power are analyzed and a practical equivalent method is proposed, but the fault components cannot be directly obtained from the proposed model.
For the analysis of a distribution network with multiple DGs, different types of equivalent models are proposed in [19] according to the fault location, and an iterative calculation method of the fault current is proposed based on a superposition theorem. However, only IIDG type DGs are considered. In [20], new energy sources with rotating characteristics are proposed as equivalent to the conventional synchronous generator model. However, it is only applicable in specific conditions and cannot be used under single-phase to ground fault.
Thus we see that an equivalent model of multi-type DGs considering variable control strategies has not been established. The impact of the LVRT control strategies on fault current characteristics has not been fully analyzed either, and the existing calculation methods are relatively weak and cannot suffice for practical applications. Therefore, it is necessary to establish an equivalent model of multi-type DGs and develop an effective calculation method for the distribution network with multi-type DGs.
The rest of this paper is organized as follows. The fault current characteristics of the IIDG and DFIG are analyzed in Sect. 2. They are controlled under LVRT requirements and with continuous excitation of the rotor-side inverter, respectively. Based on the analysis of the fault characteristics, an equivalent model of multi-type DGs is then established. In Sect. 3, a fast-iterative calculation method for the distribution network with multi-type DGs, with enhanced calculation efficiency while avoiding local convergence, is proposed based on a particle swarm optimization (PSO) algorithm. A simulation system of the distribution network with multi-type DGs is established in PSCAD/EMTDC in Sect. 4 and the simulation results verify the theoretical analysis and the proposed fast-iterative calculation method.

Fault current characteristics of the IIDG
The fault current characteristics of an IIDG are closely related to the LVRT strategy. In order to accurately analyze the fault characteristics, a LVRT strategy that meets the requirements should be established at the initial stage. A typical LVRT strategy [21] is depicted in Fig. 1. As shown, when a fault occurs resulting in voltage drop, an IIDG should remain connected to the grid at the initial 0.625 s. After 0.625 s, if the voltage remains lower than the diagonal line, the DG can be disconnected from the grid. In addition, to meet the LVRT requirement, the IIDG needs to adjust the reactive current according to the voltage drop to provide sufficient reactive support [22]: where i qref is the output current reference of the q-axis component, k represents the required reactive current coefficient, V pcc is the voltage at the point of common coupling (PCC), and I max is the maximum current of the inverter and is generally set to 120% of the rated current.
To ensure the inverter is operated within the maximum current limit, the active current reference is given as: where i dref is the reference of the d-axis active current, i d0ref is the current reference prior to the fault, and qref is the maximum active current allowed considering the inverter current limit. The three-phase fault currents of the IIDG under a symmetrical fault are [23]: where i a1 , i b1 and i c1 represent the positive-sequence currents of Phase A, B and C, respectively. i d1 and i q1 are the active and reactive current references after the fault, θ is the initial current angle of Phase A, and φ is the current angle after the fault, i.e. φ = arctan(i q1 /i d1 ).
The inertia of grid-integrated power-electronics-based units varies based on their size and control. The typical inertia time constant of an IIDG is relatively small, and is approximately 10 ms. As a result, the transient components of an IIDG after faults can be ignored and the steadystate fault component is the main factor that influences the protection threshold setting. It can be seen from (3) that the amplitude and phase of the IIDG output currents have changed and the specific changing ranges are related to the operating condition before the fault and the voltage drop after the fault. For asymmetrical faults, the amplitude and phase of the positive sequence voltage can be locked rapidly and accurately by the phase-locked loop [24]. Reactive and active current references are obtained according to (1) and (2), and the output currents can rapidly track their references through PI regulation of the inner current loop. Under asymmetrical faults, the IIDG output current is the same as that under symmetrical faults shown in (3), indicating that only positive sequence currents are generated under both symmetrical and asymmetrical faults. This is different to conventional synchronous generators. Hence, an IIDG can be equivalent to a positive-sequence current source, whose magnitude depends on the operating condition prior to the fault and the voltage drop at the PCC after the fault. (1)

Fault current characteristics of the DFIG
According to the grid integration requirements [25], a DFIG is required to provide reactive power during the fault period to support the distribution network. Under asymmetrical faults, unbalanced voltages have negative impact on the DFIG, e.g., increased losses on the stator and rotor windings, excessive heat, and emergence of pulsating electromagnetic torque. To ensure safe operation, a control strategy for balancing the stator current, which eliminates the side effects of the negative-sequence components, is applied in this paper. Thus, the rotor-side negative-sequence current can be reduced, and the second harmonic frequency components of the stator active power and electromagnetic torque can be suppressed [26]. The LVRT process of the DFIG under the control strategy is analyzed as follows.
The stator side applies the generator convention for derivation. Considering the positive and negative dq rotating coordinates, the DFIG stator voltage and flux linkage can be derived as [27]: where the superscripts '+' and '−' indicate a synchronous rotating dq coordinate under forward/positive and reverse/negative rotations, respectively. The subscripts s and r represent the respective stator and rotor components, while ω 1 and ω 2 represent the respective positive and negative dq angular velocities, i.e. ω 2 = − ω 1 . The stator voltage is set as the d-axis vector reference, and thus: where u s (1) and u s (2) represent the positive and negative sequence stator voltage components, respectively. During the fault period, a typical approach is to block the outer power control loop, so only the impact of the inner current control loop needs to be considered. In addition, since the inertia time constant of the rotor is large, the change of the rotor speed is much slower than those of (4)  (4) and (5) and ignoring the stator resistance yield: Since the reactive power is regulated based on the voltage drop, the stator reactive current reference can be obtained as: By substituting (8) into (7), the rotor reactive current can be derived as: When an asymmetrical fault occurs, the DFIG fault current is composed of both positive and negative sequence components. In order to mitigate the negative sequence component, which may lead to unbalanced heating of the stator winding, a control strategy is applied to eliminate the negative sequence components of the stator current by forcing their references to 0, i.e.: Under asymmetrical fault conditions, a double frequency (2ω 1 ) component will be generated by the negative sequence component in the forward rotating dq coordinate. This double frequency component can be filtered by a notch filter, represented by [27]: where ω n represents the notch frequency and is 2ω 1 in this paper, and ω c represents the cut-off frequency. A notch filter is added in the phase-locked loop to extract and filter the positive and negative sequence components. Similar to the positive sequence stator flux linkage, there is n the fault period. Thus, rotor current references can be obtained in the negative dq coordinate by considering (4)- (6) and (10), as: Because of the small capacity of the rotor inverter, the overall rotor current must be limited. Hence, the rotor active current reference is limited to the minimum value of the active current reference prior to the fault i rd0ref and the maximum active current allowed under the current limit, as: where I rmax represents the maximum current allowed by the rotor side inverter and is generally set to 120% of the rated current. According to (4)-(13), the DFIG stator currents under symmetrical and asymmetrical faults are: The rotor side inner current loop is controlled by PI feedforward control to track the current references, which can be designed with Type I/Type II controller to ensure: where i sd and i sq represent the d-axis and q-axis components of the DFIG stator current, respectively. The grid side converter is generally controlled with unit power factor operation to reduce the reactive current impact on the distribution network. In addition, the grid side converter eliminates the negative-sequence current under asymmetrical faults, so the output d-axis and q-axis current components i gd and i gq are: where p g0 is the active power of the grid side inverter and u gd is the amplitude of the grid voltage vector. As both d-axes of the grid side and rotor side inverters are aligned to the same grid voltage, the overall DFIG currents in the dq-axes i d and i q can be calculated as:

Equivalent model of multi-type DGs
Based on the above analysis, IIDG and DFIG can both be equivalent to a constant positive sequence current source with its amplitude and phase having a functional relationship with the positive sequence voltage at the PCC. Thus, the equivalent model can be expressed as İ DG = I∠ϕ iu with the magnitude I and the angle difference between the positive sequence current and voltage φ iu given as: where i dref and i qref are presented in Table 1 In Table 1, α and m represent the reactive current coefficient and voltage limit specified by the grid integration requirements. I max is the maximum inverter current limit, which is 120% of the rated value for IIDG and L m •I rmax /L s for DFIG. According to (18), an equivalent model of DG can be established, as illustrated in Fig. 2.

PSO-based fast-iterative calculation method of fault components for distribution networks
Based on the analysis in Sect. 2, a DG can be equivalent to a positive sequence current source under fault conditions. However, the amplitude and phase have a complex nonlinear relationship with the positive sequence voltage (17) This results in difficulties in determining the fault components with conventional calculation methods. For instance, in the conventional calculation method with separation of the real and imaginary components, the nonlinear equations can be solved but with low calculation efficiency. In addition, the nonlinearity increases with the increased complexity of the system, which may result in the conventional approach being unable to produce the results within limited time. In this paper, a fastiterative calculation method based on the improved PSO algorithm is proposed to make the nonlinear calculations faster with higher accuracy.

Improved PSO algorithm
PSO is an intelligent algorithm that imitates birds in search of food with the merits of fast convergence, small number of parameter settings and easy implementation [28]. However, PSO is prone to be trapped in local optima, i.e. local convergence. The improvement of the original PSO algorithm was first proposed in [29], and in this paper an adaptive weight coefficient method is applied to avoid local convergence. The standard PSO algorithm corrects individual behaviors to obtain the optimal solution through information sharing and individual experience among groups. The particle swarm consists of n particles flying in a D-dimensional space with a certain speed. The position of the ith particle is denoted as X i = [x i1 , x i2 ,…, x iD ] and the corresponding velocity is denoted as is an optimal position searched by the particle, while P g = [p g1 , p g2 ,…, p gD ] is the aggregation of the optimal positions currently searched by the particle swarm. The dth velocity v k+1 id and the position x k+1 id at the (k + 1)th iteration are: where c 1 and c 2 represent the acceleration constants, and ω is the inertia weight. r 1 and r 2 are random numbers between 0 and 1. The inertia weight ω affects the searching capability of the particle. When ω is small/large, the algorithm has a strong local/global search capability. The weight coefficient ω i is applied to accelerate the convergence and prevent a local optimal solution, i.e.: where ω max represents the initial inertia weight, and ω min represents the inertia weight of iteration reaching the maximum number. N i and N max represent the current and maximum numbers of iterations, respectively. In the . . where X = [x 1 , x 2 , …, x n ] T is a n×1 vector, consisting of n unknown variables. For the distribution network with DGs, X represents the required fault components and f(X) = [f 1(x) , f 2(x) , …, f m(x) ] T is an m-dimensional vector function. In order to solve (21), the unconstrained optimization approach is applied by converting it into a nonlinear least squares form as: (22) can be transformed into a 2-norn in the vector space as: (23) can then be converted into an unconstrained minimax optimization model, described as:

Fig. 2 Positive-sequence equivalent circuit of the DG
where f i (X) represents the objective function. (24) can be used as the fitness function which can be solved by the improved PSO algorithm. The fitness of the particle is inversely correlated to the function value.

Fast-iterative calculation method
A symmetrical component method can be used to establish the fitness function for the fault calculation. The fault components can then be determined with the improved PSO algorithm. The specific procedures of the calculation are explained below.

Establish the equivalent sequence circuit
To establish the equivalent sequence circuit, the DGs are replaced by positive sequence current sources, while the conventional synchronous generators are equivalent to the voltage sources in series with the reactance. For the negative and zero sequence networks, the DGs are considered as in an open state.

Establish the voltage equation at each node
According to each sequence network, there are: where subscripts (1), (2) and (0) represent the positive, negative and zero sequence networks, respectively. V i , I i and Z i are respectively the voltage, current and impedance of node i. Z DG is the impedance of the transformer and the DG branch, while I DG is the positive sequence current of the DG.

Establish constraints of the fault boundary
The constraints of the fault boundary can be established according to the fault type, e.g., for a two-phase short-circuit fault, the fault boundary constraints are: where U f (1) and U f (2) represent the fault boundary constraints in the positive and negative sequence networks, respectively.

Define the fitness function
The equivalent model of the DGs, the voltage relationship and the fault boundary constraints as respectively described in (18), (25) and (26) can be combined to establish the objective function f i (X). The fitness function can be constructed according to (24), and all particle fitness values can then be calculated. For the distribution network with the integration of an IIDG and a DFIG as illustrated in Fig. 5, the vector X has a dimension of 12 and consists of 12 variables, including the voltage/current amplitudes (Mag U-IIDG /Mag I-IIDG ) and the phase angles (Phi U-IIDG /Phi I-IIDG ) of the IIDG, the voltage/current amplitudes (Mag U-DFIG /Mag I-IDFIG ) and the phase angles (Phi U-DFIG /Phi I-DFIG ) of the DFIG, the current amplitude (Mag Is ) and the phase angle (Phi Is ) of the system, and the current amplitude (Mag If )) and the phase angle (Phi If ) at the faulty point. The objective functions are defined in Table 2 under different faults where 3P-SC-F, 2P-SC-F, 2P-G-F, and 1P-G-F represent a three-phase short-circuit fault, a two-phase short-circuit fault, a twophase to ground fault, and a single-phase to ground fault, respectively.

Initialize randomly with the improved PSO
In setting the maximum number of iterations, the number of variables and the size of the particle swarm, the speed and position of the particles are randomly initialized in the particle swarm. The constraints of the initial position are:

Update the global optimal solution
The extremum of each individual is the optimal solution obtained for each particle and a global value is identified from these optimal solutions as the global optimal solution for this iteration, which is then compared with the historical global optimization to update its value.

Update inertia weight, speed and position
The inertia weight is calculated with (20), while the velocity and position of the particle are updated and bounded according to (19).

Evaluate iteration
It is evaluated whether the precision e or maximum number of iterations is achieved. If the requirement is met, the method continues to Step 9, otherwise, it goes back to Step 6.

Calculate the fault components
The fault current components of each branch are calculated based on the calculated node voltage and network impedance according to (25). A comprehensive flow chart of the proposed fast-iterative calculation method is depicted in Fig. 3.

Simulation system
The distribution network with an IIDG and a DFIG is established in PSCAD/EMTDC, as shown in Fig. 4. Two different fault conditions at point f, i.e., a three-phase short-circuit fault and a two-phase short-circuit fault, are simulated. The proposed fast-iterative calculation (27) method is used to calculate the fault components from the established equivalent model under symmetrical and asymmetrical faults. The calculation results are compared with those of the simulations for verification. The   main parameters of the simulation system are listed in Table 3 [15,22,30].

Case studies 4.2.1 Symmetrical fault
When a three-phase fault occurs at point f at 1 s, the PCC voltages of the IIDG and DFIG are presented in Figs. 5 and 6, respectively. When the fault occurs, the voltage of the IIDG drops to 0.295 p.u. and the voltage of the DFIG drops to 0.365 p.u. within 0.03 s. The fault current characteristics of the IIDG and DFIG are presented in Figs. 7 and 8, respectively. The current references of the IIDG and DFIG are calculated based on the respective PCC voltages. From (1), (2) and (14), i dref and i qref of the IIDG are 1 p.u. and 0 p.u. before the fault, while i rdref and i rqref of the DFIG are 1.107 p.u. and − 0.44 p.u. before the fault. When the fault occurs, the d-axis current of the IIDG tracks the reference (0 p.u.) within 0.01 s and the q-axis current tracks the reference (1.2 p.u.) within 0.05 s as shown in Fig. 7. In Fig. 8, the d-axis current of the DFIG tracks the reference (0.977 p.u.) within 0.02 s and the q-axis current tracks the reference (0.698 p.u.) within 0.04 s. The maximum overshoot is 13.24% and the average tracking error after reaching the steady state is approximately 2.4%. This indicates that the output currents can rapidly track the references and conform to the LVRT strategies.
The theoretical fault current calculations are obtained by the proposed fast-iterative calculation method, and are compared to the simulation results in Table 4, while the percentage magnitude and phase errors between the theoretical calculations and simulation results are depicted in Fig. 9. These indicate that the results are very close and the maximum error is within 2%. This demonstrates the effectiveness of the equivalent model of the distribution network with multi-type DGs and the proposed fast-iterative calculation method with high accuracy under a symmetrical fault.

Asymmetrical fault
When a two-phase short-circuit fault occurs at point f, the theoretical calculations and simulation results are compared in Table 5. As can be seen, the negative sequence currents of the IIDG and DFIG are mitigated to 0.038 kA and 0.042 kA respectively by the controls. As both negative sequence currents are sufficiently small, the IIDG and DFIG can be equivalent to positive sequence current sources under the asymmetrical fault. As shown in Table 5, the calculation and simulation results of the amplitude and phase of the currents are very close. In addition, the percentage magnitude and phase errors between the theoretical calculations and simulation results as shown in Fig. 10 indicate a maximum error of less than 1%, validating the effectiveness of the equivalent model and the proposed fast-iterative calculation method with high accuracy under an asymmetrical fault. The convergence of the fast-iterative calculation method is evaluated under a three-phase short-circuit fault as an example, and the error precision of the algorithm is presented in Fig. 11. When the number of iterations increases, the error precision of the algorithm decreases, which demonstrates the convergence performance of the proposed fast-iterative calculation method. When the iteration reaches 284 times, the error precision reaches the set value e. The improved PSO algorithm exits the loop and the correct fault components are produced. The proposed fast-iterative calculation method is compared with the conventional approach with the separation of real and imaginary components in Table 6 [15,22,31]. As can be seen there, the fault current calculation results of the two approaches are very close.
The comparison demonstrates that the proposed calculation method of the fault components has high accuracy and can assist with the settings of the protection threshold. Using a computer with an Intel Core i7-6700HQ processor and 16 GB RAM, the calculation time of the improved PSO algorithm is 2.75 s, whereas the calculation time of the traditional algorithm is 42.48 s, which is over 15 times of that of the proposed approach. Hence, the calculation method proposed in this paper is fast and has high efficiency, and can meet the requirements in practical application.

Conclusions
In this paper, a distribution network with multi-type DGs has been investigated, and the fault current characteristics of the IIDG and DFIG have been analyzed. Based on an analytical study, an equivalent model of multi-type DGs has been established. A fast-iterative calculation method for fault current component calculation, which has high calculation efficiency while avoids local convergence, has been proposed. Specifically, the following conclusions are drawn.     It is noted that the PSO algorithm may not be the optimal algorithm and is unlikely to be the only approach that can solve the nonlinear equations presented in this paper. Thus, further investigation will be carried out in future research.