A Generic Power Flow Algorithm for the Microgrid Based on Time Domain Iteration Concept

The power flow calculation is an important analysis tool for the power system. The essence of the traditional power flow algorithm is to solve a set of non-linear power flow equations. If initial conditions are not properly set, it may occur convergence issues. On the other hand, as more and more power electronic devices are connected into the power grid, the traditional algorithm becomes unavailable due to the diverse operation characteristics of these devices. To solve the above issues, this paper proposes a time domain iteration (TDI) based power flow algorithm for the power electronics dominated power system, and takes the microgrid system as an example to analyze. Firstly, the proposed TDI concept is introduced, which uses time domain deduction instead of equation calculation in the power flow calculation, so that convergence issues can be avoided. Secondly, the operation characteristics of power electronic based distributed generation (DG) units are analyzed, with the corresponding DG unit, network and load models established. After this, the iteration process is presented: the DG unit models output voltage or current to the network, and receive feedback voltage or current from the latter to generate the output for the next iteration. The output generation link is independent and edited refer to the DG control strategy. In theory, any strategy can be edited into this link, making the proposed power flow algorithm more flexible. Finally, the algorithm is verified through simulation results from a 38-bus power system containing droop-controlled DG units.


I. INTRODUCTION
With the development of power electronic technology, the power electronic dominated power system has become the future trend. Compared to the traditional electromechanical device, the power electronic device has strong customized control ability, which is widely used in the control of power quality, power generation and equipment automation, etc. The application of power electronic based converters in the source, network and load links is expanding [1]- [3].
In the link of power generation, renewable energy sources (RESs), such as wind, solar, are growing rapidly and have huge potential. Due to the low energy density and fluctuation characteristics, most RESs are interfaced to the power grid through power electronic devices, which provide high energy The associate editor coordinating the review of this manuscript and approving it for publication was Bin Zhou . conversion efficiency and grid-connection performance [4]. In the link of power transmission, the high voltage direct current (HVDC) technology has become the primary choice for long-distance transmission and asynchronous interconnection due to its economic advantages and networking performance. In the link of power consumption, the converter technology has penetrated into industrial, civil and traffic loads. The motor supplied by the inverter is gradually replacing the traditional motor due to the excellent speed regulation performance; high-speed trains, subways, electric cars and DC terminals are newly added converter loads. Therefore, the scale of power electronic devices in the whole source, network and load links has gradually appeared [5]. For the power system, power flow calculation is the basis of network planning and real-time optimization. After many years of research and development, the current power flow algorithm has the advantages of high accuracy, low cost and strong robustness. Conventional algorithms include Newton-Raphson (NR) method [6]- [8], forward-backward sweep (BFS) method [9]- [11]. The essence of the NR method is to solve nonlinear power flow equations, which inevitably encounter multi-solution and convergence issues in the iterative solution process. And the BFS method is constrained by the network structure. Moreover, for future power electronic dominated power systems, the source and load characteristics are quite different from the traditional system, making conventional power flow algorithms limited. Take the microgrid system as an example, in which power is supplied by the inner distributed energy resources (DERs) through power electronic converters [12]- [14]. Normally, the droop control can be applied to distributed generation (DG) units to realize decentralized control [15]- [17]. As the voltage magnitude and frequency of DG units are determined by the real and reactive power, the bus cannot be described as a VF or PQ bus, which limits the application of the conventional power flow method. To calculate the microgrid power flow, many improved algorithms have been proposed.
In [18], a two-step power flow algorithm that accurately and efficiently represents power electronic interface DG units is proposed.
Step 1, based on the existing conventional power flow algorithm, solves for the AC network of a microgrid.
Step 2, utilizes the calculated electrical variables at the point of common coupling (PCC) of each DG unit and solves for the internal DG variables, including limits and constraints. The DG unit connected buses are analogous to slack, PQ and PV buses, so the proposed method is basically a conventional NR method. A unified power flow algorithm is proposed for the voltage-sourced converter (VSC) based DG unit [19]. In the proposed method, additional set of equations are augmented with the conventional NR power-flow equations to describe the control actions of the flexible AC transmission system (FACTS) devices; and the operational constraints/limits of the DG unit and its interface VSC are enforced after each iteration. However, only the VF and PQ control methods are discussed, the droop control method has not been mentioned. In [20], a generalized power flow algorithm based on Newton-trust region method is proposed for islanded microgrids. In this method, the generation bus is represented as droop, PV, or PQ bus, and the power sharing characteristic of droop control is utilized. However, the power sharing is an ideal state, and the droop control method is rarely used alone. For example, to deal with the reactive power sharing issue and control stability issue of droop control in low-voltage microgrids, the virtual impedance method is always used [21]- [23], which leads to the change of the DG voltage. In [24], by modelling the virtual impedance in the power flow formulation, the calculation accuracy is improved, and its basic algorithm is a globally convergent Newton-trust region method. In [25], a novel approach based on modified NR method is proposed for the droop controlled microgrids. And in [26], a generic power flow analysis approach is proposed, which adopts the symmetrical sequence component model rather than phase-coordinate model.
From the above analyses, it can be found that most power flow algorithms for the microgrid are based on NR method. That is because the traditional BFS method and most of its variants are restricted to the radial and weakly meshed networks [27], [28], while the microgrid has more diverse networks. However, the NR based methods need to modify Jacobian matrix in each iteration, taking a lot of computing time. To solve this issue, the fast decoupling NR (FDNR) method is proposed [29], [30]. By assuming that X R in the feeder line, the variable Jacobian matrix is not required. But in microgrids, the resistance component of feeder line cannot be ignored, which will lead to the inaccuracy of the calculation result. Although some modified methods have been proposed for FDNR, the calculation amount increases, or else the sensibility of R/X ratio still exists [31], [32].
To easily calculate the power flow of a power electronics dominated power system, such as the microgrid, this paper proposes a time domain iteration (TDI) based power flow algorithm. In each TDI step, the DG unit and load models output voltage or current to the network model, and receive feedback voltage or current from the latter to generate output variables for the next TDI. The advantages of the proposed TDI algorithm are as follows: 1) Compared with NR based methods, the proposed TDI algorithm does not need to solve equations, and thus avoids the convergence issue fundamentally. 2) Only constant matrixes are used in the algorithm, and there is no variable matrix, such as Jacobian matrix. So, its computing amount is limited. 3) Compared with the FDNR method, it is not sensitive to the R/X ratio, making it more suitable for the low-voltage microgrid. 4) Compared with the BFS method, the proposed TDI method is not limited to the network structure. 5) It contains an independent link for the output variable generation, where load and DG unit characteristics can be edited. So, the TDI algorithm is much flexible. This paper is organized as follows. Some basic concepts are introduced in Section II, including microgrid, TDI, and etc. The detailed modeling process is presented in Section III. As well as the working process of TDI. Algorithm validation and case study are respectively carried out in Section IV and V. Finally, conclusions are drawn in Section VI.

II. FUNDAMENTAL PRINCIPLE OF THE PROPOSED TDI ALGORITHM A. MICROGRID CONCEPT
The microgrid can be regarded as a small-scale power system. It is a typical system composed of power electronic devices, which provides an efficient integration solution for DERs and energy storage systems (ESSs). In this paper, the microgrid is taken as the study object of the proposed power flow algorithm due to its power electronic dominance. And the research results can be extended to power systems. VOLUME 8, 2020 Figure 1 shows the configuration of a microgrid composed of several DG units and loads. Each DG unit is interfaced to the microgrid with a converter. By controlling the solid-state switch (STS) at the PCC, the microgrid can operates in either grid-connected or islanded mode. In grid-connected mode, as the microgrid voltage is supported by the main grid, the DG unit can adopt either the PQ control method or droop control method to track power references. When the microgrid turns to islanded mode, the load power must be properly shared among a bank of DG units working in parallel. If the energy management is not contained, DG units should output power proportionally to their capacities. To achieve power sharing, there are mainly two control schemes. 1) Master-slave Control Scheme: In this scheme, a certain DG unit is selected as the master unit, which adopts the VF control method to support the microgrid voltage, while the other DG units adopt the PQ control method to balance the load demand. Normally, the selected master unit should have sufficient capacity and energy. A microgrid central controller (MGCC) is required in this control scheme, which works as the power divider. Although accurate power sharing can be achieved, the expansion capability and reliability of system are limited by the capacity of the master unit.
2) Peer-to-peer Control Scheme: In this scheme, microgrid voltage is supported by multi droop-controlled DG units. Due to the distributed control manner, the flexibility and stability of system are improved. Meanwhile, the communication cost is saved. However, the droop control always suffers reactive power sharing issues due to mismatched network parameters.
Moreover, by introducing the MGCC into the peer-to-peer control system, different versions of secondary regulation methods are developed to solve reactive power sharing issues and recovery voltage and frequency [33].

B. DROOP CONTROL CONCEPT
For a droop-controlled DG unit, the frequency and voltage magnitude are regulated by its output real and reactive powers, making its operating characteristics quite different from that of the traditional generator set. The droop control expressions are as follows: where f 0 and E 0 are the initial values of the frequency and voltage magnitude, respectively; P i and Q i are the measured real and reactive powers of DG unit i after the first-order lowpass filtering (LPF), respectively; D pi and D qi are the real and reactive power droop slopes, respectively. Normally, they are associated with the capacity of the DG unit: the larger the capacity, the smaller the droop slopes. So that when all the DG units operate under the same frequency and voltage magnitude, large capacity DG units will output more real and reactive powers according to (1) and (2). But in practice, only the DG frequency can be unified, the voltage magnitude can hardly be unified due to the mismatched feeder lines.

C. PROPOSED TDI CONCEPT
In this paper, the system model for power flow calculation is divided into network model, DG unit model and load model. The network model contains feeders and RL loads, which have constant impedances. According to the DG control objective, those DG units are substituted by controlled voltage or current sources. The DG unit adopting VF control or droop control is substituted by controlled voltage source (CVS), and the DG unit adopting PQ control is substituted by controlled current source (CCS). The PQ load is regarded as a specific current source, which absorbs power from the grid, just like the ESS working in charging state. The structure of the concept system model is illustrated in Figure 2, where ω com is the rotating frequency of the power system, and n is the current number of iterations. Each DG unit has its own control reference frame, and the rotating frequency is ω i . The CVS and CCS output voltageU o and currentİ o to the network model, respectively. Since the network is passive, it can immediately feed current or voltage back according to the inputs. According to this concept system model, the operation of the microgrid can be analyzed.
For a single DG unit, its operation process can be described as follows. At t = t 1 , the DG unit samples necessary electrical parameters from the network and generates voltage or current references for its local control. At t = t 2 , when the voltage or current tracking is completed, the DG unit provides set voltage or current to the network. Then the network changes its state accordingly and remains unchanged after t 13 . When next t 1 comes, the DG unit begins sample again. Assuming that all the DG units operate synchronously and the interval between t 1 and t 3 is extremely short, the microgrid operation can be treat as an iteration process as shown in Figure 3. By deducing the operation process of the microgrid, the power flow can be calculated. Different from the iterative algorithm in conventional methods, the TDI algorithm which mimics the real-time operation of the microgrid, doesn't need solve power equations. So, the convergence problem will not be encountered.

III. POWER SYSTEM MODEL BASED ON TDI A. MODELING OF THE DG UNIT
According to the control methods mentioned in Section II, the corresponding DG unit model is established. Firstly, the voltage and current of the CVS in the common reference frame can be described asU where U oi and δ ovi are the magnitude and angle of the DG voltage in the common reference frame, respectively; I bi and δ bci are the magnitude and angle of the feedback current in the common reference frame, respectively. The voltage and current of the CCS in the common reference frame can be described asU where I oj and δ ocj are the magnitude and angle of the DG current in the common reference frame, respectively; U bj and δ bvj are the magnitude and angle of the feedback voltage in the common reference frame, respectively.
As each DG unit controls its voltage and current according to the own reference frame, which has angle difference from the outer common reference frame as shown in Figure 4, the expression need be converted when voltage and current variable move from one reference frame to another. Assume that the DG unit is controlled in DQ reference frame, when the local variable in individual reference frame, such as the output voltage, goes into the common reference frame, following equations can be used whereḞ i is the voltage or current variable in the common reference frame; F di and F qi are the dq axis variables in the individual reference frame; δ i is angle difference between the two reference frames. When the grid side variable in the common reference frame, such as the feedback current, goes into the individual reference frame, following equations can be used Then, the DG units with VF, PQ and droop control method are modeled, respectively.

1) VF CONTROL METHOD
According to the characteristics of the VF control method, the voltage of the CVS in individual reference frame is where U odi and U oqi are the dq axis output voltages of CVS i, respectively; U * odi and U * oqi are the corresponding dq axis voltage references, respectively. Since the CVS voltage is controlled locally, it need be converted using (7) and (8) when goes into the common reference frame. VOLUME 8, 2020 2) PQ CONTROL METHOD According to the characteristics of the PQ control method, the current of the CVS in common reference frame can be calculated as where P * j and Q * j are the real and reactive power references of CCS j, respectively.

3) DROOP CONTROL METHOD
According to the working process of the droop control method, the real and reactive powers of the DG unit is firstly calculated.
Then the expected voltage of the droop control method is calculated as (19) and where R vi and X vi are the virtual resistance and reactance of CVS i; U vdi and U vqi are the dq axis voltage drops on the virtual impedance, respectively; i bdi and i bqi are the dq axis currents calculated fromİ bi .
The phase angle of the individual reference frame is where τ s is a time constant for the angle calculation. With the δ i obtained from (22), the output voltage in individual reference frame can goes into the common reference frame by using (7) and (8).
From the above DG unit models, it can be seen that the tracking process of the output voltage or current is neglected. And the operation process of microgrid is discretized. As the power flow analysis aims at obtaining the steady state values of the power flow, this simplification will not affect the final calculation results.

B. MODELING OF THE NETWORK
The network model to be established is sketched in Figure 5, containing feeders and RL loads. Assume that there are a 1 CVSs, a 2 CCSs, b 1 network nodes, b 2 network feeders, c 1 RL loads and c 2 PQ loads in the system. The detailed modeling process is shown as follows:

1) FEEDER MODEL
Assume that current flows from nodes j to k, the voltage drops on feeder i between nodes j and k can be expressed aṡ whereU ni andU nk are the voltage of node j and k;İ fi is the current of feeder i; and Z fi is the impedance of feeder i. If node j is connected with a CVS through feeder i, equation (23) can be rewritten aṡ (24) By combining all the feeder voltage together, the whole feeder model can be derived as (25) where

B fn U n + B fo U o = A f I f
In (25), B fn (b 2 ×b 1 ) and B fo (b 2 ×a 1 ) are mapping matrices which reflect the connecting relationship between the node and feeder. If the DG unit directly connects to the network node without any connected inductance, a relatively small impedance can be introduced between the DG and node.

2) RL & PQ LOAD MODEL
Both the PQ load and RL load are considered here. The load voltage of node i can be expressed aṡ (26) whereİ li is the current of load i; and Z li is the impedance of load i. If no RL load is connected to node i, a sufficiently large ground resistance can be introduced. Hence, the whole load model can be given by U n = A l I l (27) where . . .
If PQ load is connected to node k, it can be treated as a special CCS, which absorbs current form the network. And its current can be expressed aṡ (28) where I ik and δ ick are the magnitude and angle of the load input current, respectively.

3) NETWORK MODEL
According to the network connection, the relationship among the line current, load current and CCS current can be expressed as (29) where

I l = B lf I f + B lo I o + B li I i
. . .
In (28), the matrices B lf , B lo and B li are size of b 1 ×b 2 , b 1 × a 2 and b 1 ×c 2 , respectively, which map the connection among the RL load, PQ load, CCS and feeder. By transforming (25), (27), and (29), the complete network model can be derived as (30) where

I f = A nv U o + A nc I o + A ni I i
With the calculation result of the feeder current, the load current and node voltage of the network can be obtained from (27) and (29), respectively.

IV. ALGORITHM FLOW OF THE PROPOSED METHOD
The complete algorithm flow of the proposed TDI method is illustrated in Figure 6.
Before the proposed TDI starts, some necessary data need to be prepared, such as the DG control parameter, constant matrices and the maximum iteration number. In each iteration, the algorithm performs the following steps: Step 1: Update the state variables of the network model, and get necessary feedback variables for the calculation of input variables; If n = 1, the current variables will be zero, and the voltage variables will be the rated value.
Step 2: Calculate the input variables of the network model, including U o and I o (or I i ), using the DG unit model or PQ load model.
Step 3: According to (29), the state variable I f is calculated firstly; A nv , A nc and A ni are all constant matrices, which will not change unless the network changes; Then, calculate the other state variables U n and I l for the next iteration.
Step 4: When the iteration number reaches its maximum or the difference between two iteration results is small enough, the iteration will be ended.
With the obtained node voltage and feeder current, power flow analysis can be well performed. The power flow from node j to node k can be calculated as And the line loss of feeder i can be calculated as where R fi is the resistance of feeder i. The accuracy and effectiveness of the proposed method will be verified in the next section.

V. ALGORITHM FLOW OF THE PROPOSED METHOD
To test the proposed power flow calculation method, a test 6-bus low voltage microgrid system is established in Matlab simulation. The single line diagram of the test system is shown in Figure 7.
The corresponding network and DG unit parameters are listed in Table 1. The three DG units are identical, and each VOLUME 8, 2020  of them adopts the droop control method. Thus, the load power can be properly shared by the DG units without using communication. The simulated power and voltage of the test system are shown in Table 2. Meanwhile, two different case studies are carried out on the test system. In the first case, the test system has been solved using a conventional power flow method [33], where DG 1 is modeled as a slack bus, DG 2 and DG 3 are modeled as PV buses. In the second case study, the test system is solved using the proposed algorithm that accurately model the DG units operating in droop control. The calculation results are also shown in Table 2.
In the simulation and calculation, the reference frame of DG 1 is set as the common reference frame. It can be seen that the two calculation results are different, and only the result from the proposed method agrees with the simulation result. In calculation result 1, although the voltage magnitude and real power results agree with the simulation results, the angle and reactive power results have large deviations. The conventional power flow algorithm cannot solve the droop controlled islanded microgrid. The above calculation results validate the limitations of the conventional algorithm and the effectiveness of the proposed algorithm in solving the power electronics dominated power system.
The iteration process of the proposed method is shown in Figure 8. The result shows that the proposed method can get correct results after a few iterations. By increasing the time constant τ s , the number of iterations can be reduced, but too large time constant will make the results diverge.

VI. CASE STUDIES
In this section, a 38-bus microgrid system has been chosen to test the effectiveness of the proposed algorithm. The power flow algorithm is implemented in Matlab environment. Figure 9 shows the 38-bus balanced studied system used in [24], [34]. The feeder parameters, load nominal power, type, and the real and reactive power exponents are given in [34].  In this case study, all the loads are modeled as PQ loads. The DG unit parameters are selected to be the same as in Section V.  Table 3 shows the voltage profile, power flow and the line loss in per unit when all the DG units operate in droop mode. The reference frame of DG 1 is selected as the common reference frame. The proposed algorithm has converged after 5 iterations. The results show that the 5 DG units well share the load real power, but the load reactive power is not properly shared. This agrees with the characteristics of the droop control. The voltage profile for study case 1 is shown in Figure 10. All the node voltages are between 0.95 and 0.99. By increasing the initial voltage in voltage droop control, the average voltage of the system can be raised.

B. CASE 2: ISLANDED MICRORGRID IN MIXED MODE
In this study case, DG 1 turns to the VF mode, while DG 2-5 turn to the PQ mode (P * = 0.833S base , Q * = 0.6S base ). The voltage profile for study case 2 is shown in Figure 10.
Then, DG 1 and 2 both turn to the droop mode, while the others remain the PQ mode. It can be seen that the more the droop-controlled DG units, the closer the node voltage will be. In the three cases, the line loss of case 1 is the lowest.  Table 4 shows the voltage profile, power flow and the line loss in per unit when the microgrid operates in grid-connected mode. In this study case, the 5 DG units operate in the droop mode. It can be seen that all the DG units output their rated real power. However, due to the different voltage deviations, their reactive power outputs are quite different, and some of them even absorb reactive power. The results also show that the load real power is provided by the DG units, while the reactive power is provided by the main grid. Since the 5 DG units have output enough real power, the extra real power is transmitted to the main grid. The voltage profile for study case 3 is shown in Figure 11. The average voltage in grid-connected mode is higher than that in islanded mode. Most node voltages are between 0.99 and 1.02.

D. CASE 4: GRID-CONNECTED MICRORGRID IN PQ MODE
In this study case, the test microgrid still operates in gridconnected mode, but all the DG units turn to the PQ mode. Firstly, the DG units output rated real and reactive powers; Then, the output real and reactive powers are reduced to 85%; Finally, the output powers are reduced to 70%. The voltage profile for study case 4 is shown in Figure 12. As the result shows, as the output power reducing, the voltage of the test system is also reduced. If the output powers of DG units are too low, the voltage of the microgrid will drop out of normal range (0.95-1.05 p.u.). So, the DG units can help the grid to support the terminal voltage.

VII. CONCLUSION
In this paper, a generalized power flow algorithm based on TDI is proposed for power electronics dominated power systems. The proposed algorithm deduces the operation of the microgrid instead of solving power flow equations. The simulation results have well verified the effectiveness of the proposed algorithm. Compare with the traditional NR based methods, the proposed method needn't change matrices during iteration process, thus reducing the calculation amount. The method is not restricted to the R/X ratio and the network type, so its application is extensive. Finally, several case studies for both islanded and grid-connected microgrid test systems have been carried out to show the reliability and effectiveness of the proposed TDI power flow algorithm. The results show good convergence characteristics in all operating conditions. In the future, the enhanced TDI algorithm used for the unbalanced or nonlinear load conditions will be studied.