Quasi-Static Time-Series Power Flow Solution for Islanded and Unbalanced Three-Phase Microgrids

Recent advances in inverter-based distributed energy resources (DERs) allow microgrids to operate in grid-connected and islanded modes with ease. However, determining a steady-state load flow solution of a real-world unbalanced and islanded three-phase microgrid remains a challenge. Existing methods are either unsuitable, labor-intensive or computationally demanding for a long-term analysis of islanded microgrids. To address these issues, we propose a practical solution framework that employs externally-updated system frequency, power quantities, and droop characteristics governing the voltage and phase angle of DERs in every iteration as inputs to an off-the-shelf open-source multi-phase, multi-wire, unbalanced power flow software tool to solve an islanded-microgrid power flow. Using our framework, users only need to implement a set of droop equations and update system variables. The power flow solution is accomplished entirely by an off-the-shelf power flow tool, e.g., OpenDSS. Our proposed framework can model a microgrid of any arbitrary configuration and operating condition. We demonstrate and validate the proposed method’s efficacy by comparing its results with those modeled using a time-domain simulation with PSCAD/EMTDC. Finally, an example of a quasi-static time-series study is presented to better illustrate the application of such a tool.


I. INTRODUCTION
T HE presence of distributed energy resources (DERs), primarily inverter-interfaced distributed generation (DG) and energy storage systems (ESSs), continues to increase in distribution systems. DERs have the potential to permit portions of a distribution system to operate independently as one or more islanded microgrids. The flexibility introduced by stand-alone microgrids benefits system operation by improving reliability and power quality as well as system resiliency in the face of extreme weather events and natural disasters [1].
To realize the operational benefits, distribution system operators must update current engineering practices to include the operational characteristics of DERs. In particular, DGs and ESSs operating in an islanded microgrid must contribute to voltage and frequency regulation. Most DERs follow a hierarchical control structure divided into primary, secondary, and tertiary levels [2], similar to those of synchronous generators. At least one DER must switch its mode of operation from grid-following to grid-forming (grid-supporting) [2]. If load demand in an isolated microgrid surpasses the generation capacity, the operator must VOLUME 8, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ implement load shedding strategies to guarantee continued reliable operation. The existing protection, monitoring, management, communications, and safety schemes also need revision [3].
The immaturity of current methods for controlling and operating islanded microgrids makes adapting utility operating procedures a challenge for distribution utilities. In the literature, there are several methods for the secondary control level, classified as centralized [4], [5], decentralized [6], [7], and distributed [8], [9]. These solutions differ on the communication available within a given microgrid. Similarly, there are several methods for the tertiary control [10], [11] and load shedding [12], [13]. The methods referenced above were implemented in electromagnetic transient (EMT) simulation programs, where electrical machines, power converters, and their controls were thoroughly modeled. However, such EMT models are complex, extensive, and extremely specific. Each unique simulation requires many hours from high-skilled engineers, incurring high costs for utilities.
For power quality studies, such as voltage regulation analysis, steady-state power flow models of grid components are sufficient to capture the desired voltage and loading profiles. Unfortunately, existing power flow tools are not adequate to model inherent characteristics of islanded microgrids. They require a slack bus and lack the ability to model the contribution of DGs and ESSs to voltage and frequency regulation. Although recently published methods address the challenges associated with microgrid operational requirements, there are still limitations in their implementation. For example, methods based on backward/forward sweep [14], [15] are not efficient for systems with loops, requiring breakpoints. The symmetrical sequence component method proposed in [16] demands additional factorization for sequence-components, which increases the computational burden. The Newton trust-region method presented in [17] requires an optimization method to minimize the mismatches. This increases the complexity of the formulation and the execution time. The modified Newton-Raphson method presented in [18] requires additional equations for each DER to be included in the Jacobian Matrix. Additionally, it is unclear how to implement this method for multi-phase multi-wire systems, as only the single-phase formulation is shown. The power flow method based on neural networks presented in [19] is disadvantageous, requiring training on the method, which increases the complexity of the formulation. The methods mentioned above have not been successfully integrated into existing power flow routines. So, complete power flow and droop control formulations must be implemented by the user. This implementation is complex, especially for multi-phase multi-wire systems, making these methods less attractive for educational or industrial purposes.
To address these issues, this paper proposes a practical power flow solution framework employing a set of droop equations that govern voltage and phase angle of DERs and, along with the system frequency, to directly solve the power flow of a multi-phase multi-wire unbalanced islanded microgrid. Using our framework, users need only implement a set of droop equations. An off-the-shelf, open-source, unbalanced power flow software tool, OpenDSS, completely solves the power flow solution. As such, our proposed framework can model a microgrid of any arbitrary configurations and operating conditions. This paper is organized as follows. Section II describes the element models (e.g., loads, generators, energy storage systems) and presents the proposed formulation for the QSTS power flow. Section III compares the results obtained by the proposed formulation with the results simulated in PSCAD. Section IV demonstrates potential QSTS studies conducted using the proposed method. Section V presents the conclusions for this work.

II. DESCRIPTION OF THE PROPOSED METHOD
The power system element models (e.g., loads, DG, and ESS) and the power flow formulation described in this section are directly associated with OpenDSS elements. The formulation of multi-phase and multi-wire components is done internally in OpenDSS. It represents self-and mutually-coupled impedance of a multi-phase multi-wire system in a Y-admittance matrix. Loads and active elements are converted to equivalent current sources. Steady-state quantities of voltages are obtained by solving a system of linear equations YV = I, where V and I are vectors of steady-state voltage and current phasors. Details of the formulation and solution are presented in [20]. The additional steps added to the OpenDSS power flow solution method were developed using Python via a C-API interface [21].

A. ELEMENT MODELS 1) LOADS
The real and reactive powers of a load are influenced by the voltage magnitude and system frequency. These relationships are specified in (1) and (2) [22]: where α P , α I , α Z , β P , β I , β Z , are the ZIP coefficients for active and reactive powers, respectively, V 0 and V are the load nominal and actual (operational) voltages, P 0 and Q 0 are the nominal active and reactive powers for a given operating point, ω 0 and ω are the nominal and actual system angular frequencies, and K p and K q are sensitivity parameters between the load demand and system frequency. Typically, K p varies between 0 and 3 while K q ranges from −2 to 0 [22].

2) GENERIC GENERATOR MODEL
The ''generator'' element in OpenDSS can model DGs commonly found in microgrids. They include synchronous machines (e.g., diesel generators), induction machines (e.g., doubly-fed induction generators used in wind turbines), and grid-coupled power converters (e.g., full converter wind turbines). In QSTS simulations, this element presents a dispatch shape vector that OpenDSS multiples to calculate the nominal generator power. Users can include shapes that account for output power changes due to wind speed variation. Although photovoltaic (PV) generators can also be modeled using this element, OpenDSS has a unique element for this generator, covered in Section II-A3.
In traditional power flows, generator buses are classified according to their controlled variables: a Vθ bus when voltage magnitude and angle are known, a PQ bus when active and reactive powers are known, and a PV bus when active power and voltage magnitude are known. For islanded microgrids, an additional bus type is proposed, i.e., VF buses [18]. These buses account for DGs capability to contribute to voltage and frequency regulation. Their active and reactive powers are calculated based on their defined voltage magnitude and frequency setpoints.
For a distributed generator classified as a VF bus, the calculation of active and reactive powers uses droop equations. They correspond to the primary control level (local control), as shown in (3) [18]: where P G and Q G are the dispatched active and reactive powers for the DG, P 0 and Q 0 are the power setpoints typically set to zero, shape dis corresponds to the profile of available active power over time (such as a vector for wind speed or irradiance variation), ω ref and V ref are setpoints for the angular frequency and local voltage magnitude, and m ω p and n v q are slopes defined according to the DG sizes. The slopes for the droop equations are calculated as in (4) [17], where ω max and ω min are the maximum and minimum values accepted by the system angular frequency, V max and V min are the maximum and minimum values of voltage magnitudes, and P nom and Q nom are the specifications of each DG. The droop equations presented in (3) are typically applied to a system with a high X/R ratio. In such a system, reactances have a more significant impact than resistances, resulting in a strong coupling of P-ω and Q-V. Thus, any variation of active power has more effect on the system frequency than on the local voltage. Likewise, a variation of reactive power has more effect on the local voltage than on the system frequency. However, low voltage systems generally exhibit low X/R ratios. In these systems, there is a strong coupling in P-V and Q-ω, resulting in the use of the following droop equations shown in (5) [18]: where: Generators that are grid-coupled through power converters have a virtual impedance. Even for systems with a low X/R ratio, these inverters can emulate a synthetic inductive characteristic in their coupling [2]. As long as an adequate virtual impedance is defined, the coupling of P-ω and Q-V can be implemented regardless of the system characteristics [2].

3) PHOTOVOLTAIC GENERATORS
As mentioned in Section II-A2, PV generators are modeled using a unique element named ''PVSystem'' in OpenDSS. This element is modeled similarly to a grid-feeding converter, i.e., without a droop control. Equations (3) or (5) are included in this element to account for contribution to frequency and voltage control. Consequently, grid-feeding converters with droop control (i.e., grid-supporting) are classified in the proposed formulation as VF buses.
Similarly, grid-forming converters with droop control (grid-supporting) are classified as VF buses. Thus, (3) or (5) also apply to these elements. One grid-forming converter must be defined as the Vθ bus, to solve the power flow, while maintaining the voltage and frequency of the system according to its setpoints. If a PV generator is set as the Vθ bus, its droop equations shown in (3) and (5) need to be reorganized as follows, for P-ω and Q-V coupling, and as follows, for P-V and Q-ω coupling [18].

4) ENERGY STORAGE SYSTEMS
Energy storage systems are modeled in OpenDSS as a generic energy storage device called ''Storage.'' The operation modes of this element include idling, charging, and discharging. The state of charge (SoC) is automatically updated in each step of a QSTS simulation. Similar to the ''PVSystem'' element, the ''Storage'' element interfaces the system through a grid-feeding converter, i.e., without a droop control. The contribution to voltage and frequency regulation at the primary level (i.e., grid-supporting) is shown in 9 and 10, for coupling of P-ω and Q-V or as shown below, for coupling of P-V and Q-ω [23]. These equations are valid for grid-supporting converters classified as VF buses. The same formulation can be considered for discharging and charging modes. The droop equations for the ESS consider the SoC in their formulation, similar to the shape vector shape dis used in generators. If the SoC is zero, i.e., the ESS is empty, and this element cannot contribute to voltage or frequency regulation. Similar to PV generators (Section II-A3), if an ESS configured as a grid-forming converter is selected as the Vθ bus, its droop equations must be rearranged as follows, for coupling of P-ω and Q-V or as shown below, for coupling of P-V and Q-ω [23].

B. QUASI-STATIC TIME-SERIES POWER FLOW FORMULATION
There are three main differences between traditional power flow methods and islanded microgrid power flow methods.
In islanded microgrids: • The system frequency and voltage magnitude of the Vθ bus are no longer constant, i.e., they vary according to the system active and reactive power mismatches; • DGs and ESSs have less nominal capacity than bulk generators, so that the system losses cannot be accumulated only in the Vθ bus (commonly called the slack bus); • Multiple DGs and ESSs are considered for frequency and voltage control, classifying them as VF or Vθ buses. DGs and ESSs active and reactive power injections are updated for each iteration according to system frequency and voltage changes.
A small islanded microgrid with three buses and one ESS, PV generator, and load each, shown in Fig. 1, is used to demonstrate the step-by-step implementation of the proposed method. The power sign convention adopted herein defines that power injected into a node is positive, while power consumed is negative.
The detailed description of the proposed method is enumerated as follows:

1) INITIAL CALCULATIONS
The first step is to calculate the droop slopes (m p and n q ) for all elements that contribute to voltage and frequency control (VF and Vθ buses). The parameters ω max , ω min , V max , and V min need to be defined in this step. Equations in Section II-A, describe how (4) is applied for coupling of P-ω and Q-V and (6) is applied for coupling of P-V and Q-ω. For the example system, the ESS and the PV generator contribute to voltage and frequency control and their slopes are calculated using (4).

2) DEFINITION OF Vθ BUS
One grid-forming converter must be defined as a Vθ bus to enable the power flow solution. The bus with the lowest normalized (per unit) value for n q in X > R systems or m p in R > X systems is selected as the Vθ bus. This bus has the highest power capability to control voltages, which speeds up the convergence. For the system presented in Fig. 1, bus 3 is defined as the Vθ bus.
Similar to traditional power flow solution methods, the Vθ bus refers to the bus where the voltage magnitude and phase angle are assumed known. For islanded microgrids, the voltage magnitude at this bus is not fixed; it varies according to the element droop equations. For the system presented in Fig. 1, the voltage magnitude at bus 3 is determined using (13), the modified droop equations for an ESS shown in (11) and (12): where Q ph i mis and P ph i mis are the sum of load demand, system losses, and injected powers for the ith iteration and ph phase, excluding those elements at the Vθ bus, shown in (14). The power balance in the Vθ bus behaves similar to the slack bus in traditional power flows. If it has sufficient power/energy available, the Vθ bus will provide the Q In (13), the voltage magnitude is defined by averaging the calculated voltage in (11) or (12) and its voltage for the past iteration (V ph i−1 V θ ). It improves the convergence of the method, which is detailed in step 5 below. Also, in (13), for R > X, the SoC of ESS influences the voltage determination of the Vθ bus. As SoC reduces, the variation of the voltage bus (V

3) SYSTEM FREQUENCY UPDATE
The system frequency (ω) is determined by the power contributions of all VF and Vθ elements, according to (15) [18]: where P sys and Q sys are the system active and reactive power given by, The P sys (or Q sys ) varies for each iteration and the system frequency is updated in each iteration. For the system shown in Fig. 1, both the ESS and the PV generator are considered for frequency control. As shown in (15), the larger the power injection by a generator or an ESS is, the larger the variation of the system frequency is. The system frequency is calculated on per-phase basis to better account for the unbalanced loads in the microgrid, emulating power converts capable of injecting unbalanced powers. This calculation only represents a mathematical manipulation. The physical value of the angular frequency can be obtained by averaging the frequencies of the three phases.

4) UPDATE POWER INJECTION, CONSUMPTION, AND LOSSES
Active and reactive powers from all system elements are updated according to the equations formulated in Section II-A. The power consumption of frequency dependent loads is updated. ESSs and DGs that contribute to frequency and voltage control (Vθ and VF buses) have their active and reactive powers updated, respecting their operational limits, such as nominal power and energy capacity limit. The other elements (i.e., lines, transformers, voltages regulators, etc) are internally updated by OpenDSS using the actual state variables.
OpenDSS calculates losses during each iteration. This computation is important because losses are considered in (14) and consequently (16). The system frequency and the Vθ bus voltage are also influenced by the system losses. Therefore, all elements that contribute to voltage and frequency control receive a proportional distribution of the system losses. Removing the losses term from (14), would aggregate the system losses to the Vθ bus, similar to traditional power flows.

5) VERIFICATION OF CONVERGENCE
After calculating the system frequency, adjusting the voltage of the Vθ bus, and updating the active and reactive powers for loads and elements located at VF buses, an iteration of the power flow is executed by solving the linear system I = YV. This process is done internally in OpenDSS. Then, steps from 2 to 5 repeat until convergence.
Conceptually, the power flow converges when Q mis , calculated in (14), (or P mis for systems with R > X) and Q V θ , calculated in (3) or (9), (or P V θ calculated in (5) or (10) for systems with R > X) are equal. In the system from Fig. 1, Q mis needs to equal Q V θ as calculated in (9) to converge, because the element in the V θ bus is an ESS. When both variables have the same value, the system power mismatch is fully supplied by the element at the V θ bus. As a consequence, the values for node voltages have a marginal update when this condition is reached. The frequency variation for each iteration is also verified. The power flow converges when the maximum changes of the voltage magnitudes for all buses and the system frequency are lower than a specified tolerance.
As explained in step 2, the voltage of the V θ bus must be modified in each iteration to match Q mis and Q V θ . In (13), the calculated voltage for the V θ bus is averaged with its voltage for the past iteration. This reduces the voltage steps for the initial iterations. As a result, for the initial iterations, variations on reactive power (X > R) or active power (R > X) at VF elements are also reduced, expediting convergence. Fig. 2 shows an example of the convergence process for the system illustrated in Fig. 1. Without voltage averaging, VOLUME 8, 2021  the reactive power in the small system oscillates indefinitely. With the voltage averaging, the system quickly converges. Fig. 3 summerizes the steps required to run an islanded microgrid power flow using OpenDSS. The internal steps for solving the linear equations and the process to advance to the next time step required in a QSTS simulation are included, even though they are internal to OpenDSS. They are included for clarity of the QSTS power flow process.

III. VALIDATION OF THE PROPOSED METHOD USING PSCAD
The IEEE 34-node test feeder, shown in Fig. 4, is used to validate the proposed method by comparing the power flow results with those solved using PSCAD. The original system is reduced by removing 14 buses and the voltage regulators but maintaining the original loading. The purpose of the reduction is to simplify the PSCAD system model. Four PV generators are added to the system and the one at bus 812 is assigned as the Vθ bus. The relationship adopted for the droop equations is P-ω and Q-V (Section II-A). Table 1 compares the node voltages obtained using OpenDSS (the proposed method) and PSCAD for all buses for the reduced IEEE 34-node test feeder. Positive errors indicate that OpenDSS values are higher than PSCAD ones.
Both voltage magnitudes and angles match well, with an absolute maximum error of 0.2769% and 0.2 • , respectively. The active and reactive powers injected by the PV generators had an absolute maximum error of 0.2545% and 1.2346%, respectively, as shown in Table 2. The system frequency error was −0.035%.
The accuracy of the proposed method presented in this section matches the precision of the existing methods in literature [14]- [18]. The main advantages of the proposed method are its straightforward implementation in an existing power flow software (i.e., OpenDSS) and its capability to handle multi-phase, multi-wire systems without increasing the formulation complexity. These characteristics are lacking in the existing methods presented in the literature.

IV. EXAMPLE OF A QSTS STUDY USING THE PROPOSED ISLANDED MICROGRID POWER FLOW
This section demonstrates the proposed QSTS power flow solution method on using the complete IEEE 34-node test feeder in both a grid connected and islanded state.
Originally, the loads in this system do not vary with time. For a QSTS analysis, a residential load shape with hourly resolution from [24] is applied to all loads. An irradiance curve corresponding to a clear-sky day is also included. The hourly resolution is adopted for the following QSTS simulations and ω 0 = 376.99 rad/s, ω max = 1.01ω 0 , ω min = 0.99ω 0 , V max = 1.05 pu, and V min = 0.95 pu. Finally, six PV generators and three ESS were included in this system, with the specifications detailed in Table 3. Bus 860 was chosen as the Vθ bus and has a PV generator.
The settings of the three ESS included in the test system are defined based on the system operation mode. For the grid-connected mode, all ESSs are programmed as follows: set to idle from 0h to 9h, charge at 20% of their nominal power from 11h to 15h, and discharge at 10% of their nominal power for the rest of the time. For the islanded microgrid mode, all ESSs are configured to charge at 20% of their nominal power rating during the highest period of PV generation (from 10h to 15h) and discharge according to their droop equations for the rest of the day.

A. SYSTEM FREQUENCY AND VOLTAGE RESULTS
The variation of the system frequency is shown in Fig. 5. The behavior of the system frequency depends on the total system loading. Due to the droop control characteristics, the higher the loading, the lower the system frequency is. The limits VOLUME 8, 2021  defined for the system frequency can be adjusted by changing the parameters ω max and ω min .
The grid-connected and islanded microgrid voltage profiles are presented in Fig. 6. The system operates from 0h to 5h and 7h to 17h. However, there is not sufficient power or energy capacity to sustain all loads at 6h and from 18h to 0h. A load shedding could maintain continuous operation of the system. The proposed power flow method can be utilized to design load shedding strategies for multiple systems in a fast and efficient manner. The voltage profile shows that the ESSs fail to ensure that the islanded system remains within regulatory limits. These findings reinforce how useful the proposed method can be to update existing standards and best practices for islanded microgrids.

B. RESULTS OF POWER SHARING AMONG GENERATORS
Power-sharing among elements is categorized by asset class (i.e., PV and ESS classes) as shown in Fig. 7. As expected, the active power contribution is dominated by the ESSs during the period without PV generation (0h to 5h), which gradually decreases as the PV generation increases. The period from 11h to 15h is dominated by PV generators because the ESSs are set to charge. PV generators ability to inject reactive power, when not injecting active power, allows them to dominate reactive power sharing. During the ESS charging period, the reactive power is entirely provided by the PV generators, since the ESSs charges with a unity power factor.
The normalized active power for each asset is shown in Fig. 8. The normalization is calculated using the maximum power point tracking (MPPT) for PV generators and the nominal active power for ESSs. This figure shows that the loading of the PV generators is proportionally the same except for the PV on bus 822. System power imbalances during periods with low PV generation caused this outlier (e.g., from 7h to 9h). However, it has a marginal effect on the PV loading because when the available active power for this asset becomes relevant (i.e., between 10h to 17h) its loading reaches values equivalent to the other generators.
For ESSs, the active power loading depends on the energy capacity of each element, as stated in (9) and (10). Although the ESS installed at bus 828 and 846 have different settings, their normalized load sizes are similar. Their ratio of energy capacity and nominal active power is equal. Analogously, the ESSs installed at buses 800 and 828 share similar specifications, but since their energy capacities are distinct their normalized active powers are different.

C. NUMBER OF ITERATIONS, PROCESSING TIME, AND SCALABILITY
The number of iterations required to converge each power flow and the respective processing time is shown in Fig. 9.   OpenDSS has a tolerance of 1 × 10 −4 . The simulations were executed on a Ryzen 7 XT CPU 3.9 GHz, with 32 GB of RAM. The first power flow instance (i.e., 0h) requires the highest number of iterations because it provides the first solution for the islanded microgrid. Subsequent instances, require fewer iterations as they are initialized with the solution of the previous instance. Fewer iterations are required when the ESSs are charging (between 11h and 15h), as there are fewer elements contributing to voltage and frequency regulation. It takes only 18s of processing time to execute this daily power flow. Therefore, the proposed method is viable in extensive studies, such as stochastic assessments or the evaluation of several systems.
To demonstrate the reproducible and scalable nature of this power flow formulation for islanded microgrids, a realworld 11.9-kV feeder from a Brazilian utility [25] operating in islanded mode is analyzed using the proposed method. This feeder has 1,815 buses, 646 lines, and the peak load is around 5 MW. The X/R ratio of the lines vary from 0.32 to 2.89. The topology of the system is shown in Fig. 10. The same PV generators and ESS from the 34-node test feeder simulation were added, with their ratings doubled to supply the load. Fig. 11 shows the number of iterations and processing time of a QSTS daily simulation of this feeder.
The total simulation took only 24.5s. This indicates that the processing time of our method marginally escalates with the system size, as the number of nodes increased by 3,226% and the processing time increased only 36%. It represents another advantage compared to the methods proposed in the literature.

V. CONCLUSION
This work proposes a quasi-static time-series power flow for islanded microgrid systems. The proposed formulation utilizes the solution mechanism from OpenDSS (i.e., a Yadmittance matrix formulation), making it adequate for solving multi-phase, multi-wire, and unbalanced systems. A major advantage of the proposed method is its simplicity in adapting the traditional power flow from OpenDSS into a suitable formulation for solving the power flow of an islanded microgrid.
The proposed method is validated through time-domain simulations using PSCAD/EMTDC. The accuracy obtained on these simulations is equivalent to the existing methods presented in the literature. In a quantitative analysis, the maximum error obtained for voltage magnitude and angle (i.e., power flow state variables) was 0.0027 pu and 0.2 • , respectively. A real world simulated example of a QSTS power flow study on an islanded microgrid was presented for a better understanding of the applications and opportunities for this tool. These studies also demonstrated the practicality of implementing the proposed method of the QSTS power flows for islanded microgrids.