A Homotopy-Based Approach to Solve the Power Flow Problem in Islanded Microgrid with Droop-Controlled Distributed Generation Units

Abstract

This paper proposes a homotopy-based approach to solve the power flow problem (PFP) in islanded microgrid networks with droop-controlled distributed generation (DG) units. The technique is based on modifying an “easy” problem solution that evolves with the computation of intermediate results to the PFP solution of interest. These intermediate results require the solution of nonlinear equations through Newton–Raphson (NR) method. In favor of convergence, the intermediate solutions are close to each other, strengthening the convergence qualities of the technique for the solution of interest. The DG units are modeled with operational power limits and three types of droop-control strategies, while the loads are both magnitude voltage- and frequency-dependent. To evaluate the method performance, simulations are performed considering the proposed and classical NR methods, both departing from a flat start estimation. Tests are carried out in three test systems. Different load and DG unit scenarios are implemented for a 6-, 38-, and 69-bus test system. A base case is studied for all systems, while for the two larger models, a loading factor is used to simulate the load augmenting up to the maximum value. The results demonstrated that for the largest-size model system, only the homotopy-based approach could solve the PFP for stringent requirements such as the diversification of the load profile and hard loading operation point.

1. Introduction

In the last decades, transmission/distribution companies, governments, and society have been working to overcome power shortages, energy diversification, energy saving, and eliminating environmental pollution due to fossil fuels, which can not be solved only through the traditional electrical grid [1]. Distribution generation (DG) units, such as full cell technology [2], photovoltaic systems, and wind turbines [3], have been added to the existing grids. Places that need a continuous power supply, such as hospitals, educational centers, airports, and ports, may be eligible for using DG. The penetration of DG into the existing distribution grids has led to the concept of microgrids. Microgrids (MG) are formed by the accumulation of DG units, energy storage systems [4], and loads that operate in conjunction with each other to ensure a reliable power supply to the microgrid network [5]. The hybrid energy sources integrated into the MG requiring embedded DC/AC, DC/DC, AC/DC, and/or AC/AC converters, with different management control, as fuzzy logic controllers [6,7] to supply the AC grid. The choice of control strategies can have technical and financial impacts on the network [8]. They can operate in grid-connected mode by connection to the distribution grid through a point of common coupling or independently of the distribution grid in islanded mode. The last operation mode exists when an outage occurs in a previously grid-connected microgrid. Thus, the islanding phenomenon is one of the greatest challenges for microgrids. Associated with challenges are low-inertia [9]; system stability [10]; power sharing among DGs [11]; the maintenance of voltage and frequency [12,13]; and, predominantly, power flow (PF) studies [14,15,16].

Traditional PF solvers, such as one that uses the Newton–Raphson (NR) method, and their decoupled versions have been used in the grid-connected operation mode. In this case, the main grid controls both the voltage and frequency [17] of the microgrid (distribution system). In operation in the islanded mode, the voltage and frequency of the microgrid must be controlled according to load variations through appropriate control techniques. The insertion of DG units has gained considerable interest with the increasing use of renewable energy sources (RES). Associated with the DG units, an essential strategy for controlling the power distributions and frequency is the operation in the droop mode control [18,19]. It should be noted that the computational tools used in PF studies for islanded microgrid systems have the following particularities: an absence of slack bus, frequency as an inherent variable in load flow, and generators with the option to control the generated power and frequency through droop units. In large power systems, the conventional PF is a computational tool for planning, analyzing, controlling, and operating purposes. Generally, PF consists in solving nonlinear equations relating voltage and power injection in each system node but in the synchronous frequency. The solution of these equations represents the steady state conditions. Conventional PF programs cannot solve the problem of isolated microgrids. Thus, several works have been proposed for this aim. In [20], the conventional power flow in conjunction with a large distribution generator was proposed to simulate a slack bus to solve an isolated microgrid problem. Recently, research has been conducted on the developing of power flow algorithms for isolated microgrids [14,21,22,23].

In the literature, various methods and techniques have been proposed for obtaining the PF solution of the islanded microgrids via Newton methods [5,24] and their derivatives [18,21]. According to [25,26], the power flow algorithms for islanded microgrids allow for the following classification:

• Jacobian-based methods—these include NR methods and their variants to solve the nonlinear system related to power flow through Jacobian matrix calculation [18,21,27];
• Backward/forward methods—these use the Kirchhoff law [28,29] and radial topology to attractively solve the power flow equations based on the fixed-point method;
• Gauss-Zbus methods—these use an iterative process to obtain the solution of a power network through nodal matrix analysis [30,31];
• Approximations—these utilize the first order Taylor expression to obtain the Jacobian matrix through nodal matrix analysis using the current injection [32,33].

Among the methods listed previously for obtaining the power flow solution in islanded microgrids, different techniques are available in the literature [5,18] based on NR approaches. In the present day, industry and researchers have been using these approaches with different electrical variables [12,14,34,35]. In [5], an algorithm is presented based on the Newton trust region (NTR) that includes the consideration of the distinguishable characteristics of a microgrid to overcome the limitations of conventional power flow methods. An improved modified NR (IMNR) is proposed in [34]. It proposes to extend the traditional Newton–Raphson technique to the islanded microgrid case with complex loads. A nested-iterative Newton–Raphson (NINR) for islanded microgrids was proposed in [36] to adjust the voltage of a slack bus and grid frequency to force active and reactive power flowing through a slack bus be zero. In turn, in [18,37], a non-linearity in the droop control equation is presented with a modified NR (MNR) method. Recently, linear-based NR approximation [32] and linear power flow formulation [38], both based on NR and depending on the Jacobian matrix for a solution in islanded microgrids, have been proposed.

In light of the previous discussion, most of these algorithms are well-suited for the operating conditions of microgrids. However, in some cases, the application of these different algorithms is limited or restricted to various constraints/limitations. In the context presented in this paper, motivate the research by a new power flow algorithm for microgrid use. A technique to solve complex nonlinear systems as the power flow problem is homotopy. From a mathematical point of view, homotopy theory systematically studies situations where maps of a “difficult” problem are handled [39]. Despite this theory originating from abstract aspects of algebraic topology [40], nowadays, this theory has also been used in other areas of mathematics, such as algebraic geometry (for example, homotopy theory), electric circuits [41,42], and electric power networks [43]. In [44], homotopy methods and their numerical implementation for a load flow multi-solution of power systems are investigated. Since then, several applications of power flow studies of electrical networks with DGs [43] have been studied. In [45], a set of sufficient conditions for the uniqueness of the power flow solution is proposed using the uniqueness of a continuation path in a homotopy method. Recently, a flat start guess homotopy-based power flow method was investigated in [27]. In this publication, the homotopy technique was used to solve ill-conditioned power flow problems in grid.

The isolated microgrid for relatively large-scale systems tends to present hard PFP convergence when solved by the traditional Newton–Raphson and departing from a flat start guess. Motivated by this fact and the publication in [27], this paper presents a homotopy-based approach to solve the power flow for isolated microgrid considering droop-controlled DG units and different load characteristics and loading factors. The PF is conceived of so that a flat start guess always initializes the problem formulation. However, a final solution is successfully obtained through the computation of intermediate PF solutions along the homotopy path. Simulations using an isolated microgrid system demonstrate that the traditional NR solver cannot solve the same problem. The problem considers detailed modeling for loads and DGs, as found in several publications [5,18,46]. The loads are represented in different forms through an exponential model considering voltage and frequency dependency. The DG units are droop-control mode dependent and have the traditional inductive, resistive, or complex impedance characteristics. Simulations considering the proposed approach and the traditional NR method with flat start initialization are performed in three test systems with 6, 38, and 69 buses.

The primary contributions of this paper are:

• The presentation of a homotopy-based approach appropriate to deal with the PFP for hard loading scenarios, considering a frequent variation of DG unit limits and diversity in load modeling;
• The use of the NR method to determine the homotopy path solution (intermediate results) considering the limitations of the DG units at each point of the path;
• The exemption from the analytical calculation of the terms of the Jacobian matrix, in preference to the computation through numerical derivatives, thus facilitating the manipulation of frequent calculations when variable limits are exceeded;
• The initialization of the PFP always departing from a flat start guess.

The paper is organized as follows. Section 2 details the primary element models composing the isolated microgrid, i.e., three DG unit droop-control strategy models are addressed, besides a series feeder circuit and a magnitude- and frequency-dependent load model. Additionally, a discussion on the reference bus of the isolated microgrid is presented. Section 3 introduces the power flow problem. It is shown that the power balance equations can be handled as a complex-valued function (power injections), and the Jacobian matrix entries can be computed numerically, which is beneficial for including equipment operational limits and system operational conditions. This section also addresses the proposed homotopy-based technique and outlines the details that are required to compute the points of homotopy path through the NR method. Tests and experiments are assessed in Section 4. In this section, several simulations are presented in three systems with 6, 38, and 69 buses to demonstrate the performance of the proposed technique compared with the classical NR method. Finally, Section 5 presents the conclusions of the paper.

2. Islanding Microgrid Modeling

This section presents the basics of static MG network modeling. The fundamentals are well documented in [5,18,46] and other publications in such a way that additional information can be accessed in those references. The focus is on the power flow in the microgrid, considering studies with different load profiles and DG unit representations. Therefore, a primary aspect to consider is the components’ modeling for use in the power flow problem (PFP).

It is well-known that the treatment given to the isolated microgrid load flow problem differs from that used for the classical grid. Mainly because, in the case of MG, the frequency becomes an additional variable of the problem. In this case, the components also depend on this variable, which should be the same in all parts of the MG. A second aspect that strongly impacts pursuing PFP solutions is that distributed generation units contribute with small generation amounts. Then, imposing limit values for all generation units is mandatory because a possible hard reverse power flow in the network depends on the loads’ on/off status and DGs. Thus, nonlinearity associated with the ceiling values of distributed generators can dramatically affect the search for the power flow problem as performed in cases where limits are little required.

In conventional PFP, DG units are modeled through a PV/PQ type representation. However, when this type of representation is preferred, it is ignored that DG unit responses are sensitive to both voltage and frequency variations. However, these aspects must also be incorporated into the representation of DG in isolated microgrids into a denomination known as droop, besides the generation limits. Therefore, considering all of these aspects, a DG unit needs to be adjusted locally.

Balanced power flow calculations are assumed in such a way that the focus of this paper considers only data and representation for the positive sequence network.

2.1. Distributed Generator Modeling

All DG units represented in this paper are assumed to be controlled through a droop scheme, with each DG having its maximum active and reactive power limited. The DG unit operation droop allows for the proper supply of load in MG. Therefore, droop modeling must be properly represented.

The droop representation must be addressed according to the DG type of output impedance [18], i.e., depending on how the real and imaginary part of the output impedance is considered, a droop type is characterized as resistive, inductive, or a combination of both types. Each DG unit is treated as connected to a PQ bus for the modeling of the PFP.

2.1.1. Inductive-Type Droop DG Unit

In this type of characterization, the imaginary part of the DG output impedance is greater than the resistive. For this situation, the droop equations representing the DG active and reactive power are dependent on frequency and terminal voltage magnitude, respectively. According to [18], for the ith DG, the output powers are calculated as

$$P_{Gi}-P_{Gi}^0=\frac{1}{m_{pi}}({\omega }_0-\omega ),$$

(1)

$$Q_{Gi}-Q_{Gi}^0=\frac{1}{n_{qi}}(V_{refi}-V_i),$$

(2)

where $P_{Gi}$ and $Q_{Gi}$ are the injected active and reactive powers of the ith DG unit into the terminal bus in pu, and $P_{Gi}^0$ and $Q_{Gi}^0$ are the active and reactive power setpoints, which are set zero in this paper. ${\omega }_0$ is the nominal frequency in pu and set in 1.0 pu; $\omega$ is the variable frequency in pu; $V_{refi}$ is the voltage magnitude in the reference bus; $V_i$ is the voltage magnitude in the DG unit terminal bus; and $m_{pi}$ and $n_{qi}$ stand for frequency droop and voltage magnitude droop coefficients, respectively, in pu(rad/s)/pu(power) or simply pu. The latter notation will be adopted in this work.

2.1.2. Resistive-Type Droop DG Unit

The distribution network can be predominantly resistive. This can be verified because of the absence of the coupling inductor or the presence of highly resistive feeders in the distribution network. For this situation, the droop type is called “inverse droop”, and now active and reactive powers are associated with voltage magnitude and frequency, respectively [18]. Power Equations (3) and (4) below characterize this type of DG unit:

$$P_{Gi}-P_{Gi}^0=\frac{1}{n_{qi}}(V_{refi}-V_i),$$

(3)

$$Q_{Gi}-Q_{Gi}^0=\frac{1}{m_{pi}}(\omega -{\omega }_0).$$

(4)

2.1.3. Complex-Type Droop DG Unit

This situation occurs when no assumption is assumed concerning the DG output impedance, and it can be considered complex. In this case, the active and reactive power of the DG unit are not decoupled as presented in the previous models. Therefore, following the model of DG unit output power suggested in [18], the expressions adopted for the power are:

$$P_{Gi}-P_{Gi}^0=\frac{1}{2}\left[\frac{1}{m_{pi}}({\omega }_0-\omega )+\frac{1}{n_{qi}}(V_{refi}-V_i)\right],$$

(5)

$$Q_{Gi}-Q_{Gi}^0=\frac{1}{2}\left[\frac{1}{n_{qi}}(V_{refi}-V_i)-\frac{1}{m_{pi}}({\omega }_0-\omega )\right].$$

(6)

For all DG unit power outputs in (1)–(6), the outputs are constrained according to the inequalities $0\le P_{Gi}\le P_{Gmaxi}$ and $Q_{Gmini}\le Q_{Gi}\le Q_{Gmaxi}$, in which $P_{Gmaxi}$ is the maximum active power set for the ith DG unit; $Q_{Gmini}$ and $Q_{Gmaxi}$ are the minimum and maximum reactive power set for the ith DG, respectively.

2.2. Feeders Modeling

The system frequency $\omega$ is a variable in islanded MG. Then, the impedance considered in the modeling of line and feeders in the conventional PFP must be modified according to the frequency. The changes need to be implemented only in the reactance because assuming that the interest is in the computation of frequency near the nominal, the resistance of the cables is kept constant. Additionally, the capacitance of the cables is neglected. Then, for a feeder interconnecting buses $\#i$ and $\#j$, the impedance is $Z_{ij}=R_{ij}+j\omega L_{ij}$, in which $R_{ij}$ and $L_{ij}$ are the resistance and inductance of the feeder, respectively, and it is assumed that $L_{ij}=X_{ij}^0/{\omega }_0$, where $X_{ij}^0$ is given for the nominal frequency ${\omega }_0$.

2.3. Load Modeling

Active and reactive power demanded in the distribution system characterizing the loads can be represented as a function of the bus voltage magnitude $V_i$ and angular frequency deviation. A well-accepted model is given by the exponential form [5] is:

$$P_{Li}=P_{0i}{\left[\frac{V_i}{V_{0i}}\right]}^{{\alpha }_i}(1+K_{pi}\Delta \omega ),$$

(7)

$$Q_{Li}=Q_{0i}{\left[\frac{V_i}{V_{0i}}\right]}^{{\beta }_i}(1+K_{qi}\Delta \omega ),$$

(8)

where $\Delta \omega =\omega -{\omega }_0$ is the angular frequency deviation in pu; $P_{0i}$ and $Q_{0i}$ are the active and reactive power measured for the voltage $V_{0i}$, respectively; ${\alpha }_i$ and ${\beta }_i$ are the active and reactive power exponents, respectively; and $K_{pi}$ and $K_{qi}$ are frequency load coefficients in pu.

2.4. Reference Bus and System Losses

The modeling of the isolated MG requires a reference bus. It differs from the conventional PFP, in which a slack bus is assigned. In the case of the slack bus, the voltage is wholly defined in this bus, and it supplies the power to balance loads and losses of the whole network. In the case of an isolated MG, only the angular reference is set, and the voltage magnitude of the bus needs to be computed as a state variable. However, the reference bus alone cannot supply the power to balance the losses in the system. On the other hand, all DG units and PV buses connected to the MG must contribute to this balance. In this paper, we assume the connection of only DG units in the isolated MG, in such a way that the total power generated by DGs are:

$$P_G^{total}=\sum _{k=1}^{N_{DG}}{P}_{Gk},$$

(9)

$$Q_G^{total}=\sum _{k=1}^{N_{DG}}{Q}_{Gk},$$

(10)

in which $N_{DG}$ is the number of DG units, and $P_{Gk}$ is the injected active power into the bus $\#k$ according to a type of DG unit defined in (1), (3), or (5); similarly, $Q_{Gk}$ is the injected reactive defined in (2), (4), or (6). All individual powers are limited to meeting the physical requirements of each DG unit.

3. The Power Flow Problem

The standard power flow problem for a grid network of $N_b$ buses is formulated to determine the nodal voltages at the network [47]. However, its formulation for an isolated MG must also include the frequency as a state.

3.1. Power Balance Equations

In the isolated MG, the problem is similar to the standard PFP except that now the voltage magnitude in the reference bus and the frequency $\omega$ must be considered, in which $N_{DG}$ is the number of DG units and $P_{Gk}$ is the injected active power into the bus $\#k$ according to a type of DG unit defined in (1), (3), or (5); similarly, $Q_{Gk}$ is the injected reactive defined in (2), (4), or (6). All individual powers are limited to meeting the physical requirements of each DG unit. Mathematically, the power balance equations (PBEs) to solve the PFP for an isolated MG can be presented from the current Kirchoff law in each bus (11), combined with the power injection for each bus (12).

$${\overline{I}}_k\left(\omega \right)={\overline{\mathbf{Y}}}_{bk}\left(\omega \right)\overline{\mathbf{V}}\left(\omega \right),k=1,2,\dots ,N_b,$$

(11)

$$0={\overline{V}}_k\left(\omega \right){\overline{I}}_k^{*}\left(\omega \right)-{\overline{S}}_k(V_k,\omega ),k=1,2,\dots ,N_b.$$

(12)

In (11) and (12), ${\overline{I}}_k\left(\omega \right)\in \mathbb{C}$ is the current injection in bus $\#k$; $N_b$ is the number of bus in the isolated MG; ${\overline{\mathbf{Y}}}_{bk}\in {\mathbb{C}}^{1\times N_b}$ is the kth row of the matrix of nodal admittance (${\mathbf{Y}}_{BUS}$); and $\overline{\mathbf{V}}\left(\omega \right)\in {\mathbb{C}}^{N_b}$ is the vector of nodal voltages, whose entry is ${\overline{V}}_k\in \mathbb{C}$.

Note that (12) comprises a set of $N_b$ complex-valued equations. Therefore, it consists of $2N_b$ real-valued equations. Additionally, the voltage phasor in the bus $\#k$ is ${\overline{V}}_k=V_k\angle {\theta }_k$, for which $V_k\in \mathbb{R}$ is the voltage magnitude and ${\theta }_k\in \mathbb{R}$ is the voltage phase angle.

From (12), we can form the $2N_b$ real-valued equations:

$$g_{rk}\left(\mathbf{x}\right)=real\left[{\overline{V}}_k\left(\omega \right){\overline{I}}_k^{*}\left(\omega \right)-{\overline{S}}_k(V_k,\omega )\right],{\mathbf{g}}_r\left(\mathbf{x}\right)={\left[g_{r1}{g}_{r2}\dots g_{r{N}_b}\right]}^T,$$

(13)

$$g_{mk}\left(\mathbf{x}\right)=imag\left[{\overline{V}}_k\left(\omega \right){\overline{I}}_k^{*}\left(\omega \right)-{\overline{S}}_k(V_k,\omega )\right],{\mathbf{g}}_m\left(\mathbf{x}\right)={\left[g_{m1}{g}_{m2}\dots g_{m{N}_b}\right]}^T.$$

(14)

Finally, a $2N_b$ real-valued equation set is formed as

$$\mathbf{g}\left(\mathbf{x}\right)={\left[{\mathbf{g}}_r^T\left(\mathbf{x}\right){\mathbf{g}}_m^T\left(\mathbf{x}\right)\right]}^T=0,\mathbf{x}={\left[{\theta }^T{\mathbf{V}}^T{V}_{ref}\omega \right]}^T\in {\mathbb{R}}^n,\mathbf{g}\left(\mathbf{x}\right)\in {\mathbb{R}}^n,$$

(15)

where $n=2N_b$; $\theta$ is the vector with the phase angles in PQ buses (loads and DG units); $\mathbf{V}$ is the voltage magnitude vector for PQ buses; $V_{ref}$ is the voltage magnitude in the reference bus; and the phase in the reference bus is assumed zero and is not included in the vector $\theta$ in (15).

In (12), the complex power injection is defined according to the active and reactive powers computed for the load or DG unit. For a DG in the bus $\#k$, ${\overline{S}}_k(V_k,\omega )=P_{Gk}+j{Q}_{Gk}$ (see expressions (1)–(6)). However, for a load ${\overline{S}}_k(V_k,\omega )=-P_{Lk}-j{Q}_{Lk}$ (see expressions (7)–(8)).

In summary, the PFP for isolated MG can be resolved following the modeling for the components and determining the solution of the nonlinear equations system defined by (15) and organized as

$$\mathbf{g}\left(\mathbf{x}\right)=0,\mathbf{x}={\left[{\theta }^T{\mathbf{V}}^T{V}_{ref}\omega \right]}^T\in {\mathbb{R}}^n,\mathbf{g}\left(\mathbf{x}\right)\in {\mathbb{R}}^n.$$

(16)

3.2. Solving the Nonlinear Equations System by the Newton–Raphson Method

The nonlinear equations system in (15) can be resolved by an iterative method. In the sequence, we describe how the PFP is solved by employing the classical NR method.

In [18], the problem was solved by the NR method. In that work, the elements of the Jacobian matrix, $\mathbf{J}\left({\mathbf{x}}^{\left(i\right)}\right)=\frac{\partial \mathbf{g}\left(\mathbf{x}\right)}{\partial \mathbf{x}}$, in $\mathbf{x}={\mathbf{x}}^{\left(i\right)}$, are analytically determined. We propose to solve the procedure carrying out the computations of the elements numerically. This procedure is effective when several limit violations are present along the iterative process, where frequent analytic formulae switching to limit values is verified.

For the implementation of the numerical computation of the Jacobian matrix elements, we have assumed a perturbation $h={10}^{-6}$ for each variable and then used numerical differentiation [48] to compute the derivative of $\mathbf{g}\left(\mathbf{x}\right)$ in (15) in relation to the perturbed variable numerically. In case of the derivative for the variable $x_k^{\left(i\right)}$, the perturbed variable has the value $x_k^{\left(i\right)}+h$ and the kth column of the Jacobian matrix is calculated as $\mathbf{J}(:,k)=\frac{\mathbf{g}({\mathbf{x}}^{\left(i\right)}+h)-\mathbf{g}\left({\mathbf{x}}^{\left(i\right)}\right)}{h}$. Therefore, the use of explicit analytic derivatives is given up.

The iterative solution of the NR method [47] for the ith iteration is obtained as

$${\mathbf{x}}^{(i+1)}={\mathbf{x}}^{\left(i\right)}+\Delta {\mathbf{x}}^{\left(i\right)},$$

(17)

where, considering a given initial estimate, ${\mathbf{x}}^{\left(0\right)}$, the increment $\Delta {\mathbf{x}}^{\left(i\right)}$ is computed for the mismatch $\mathbf{g}\left({\mathbf{x}}^{\left(i\right)}\right)$ as

$$\Delta {\mathbf{x}}^{\left(i\right)}=-{\left[\mathbf{J}\left({\mathbf{x}}^{\left(i\right)}\right)\right]}^{-1}\mathbf{g}\left({\mathbf{x}}^{\left(i\right)}\right).$$

(18)

The update ${\mathbf{x}}^{(i+1)}$ is computed, while a given tolerance for the mismatch is not reached, i.e., $\left|\right|\mathbf{g}\left({\mathbf{x}}^{(i+1)}\right){\left|\right|}_{\infty }<ϵ$. After, DGs’ powers are calculated and checked to determine whether the limits are satisfied. Figure 1 exhibits a flowchart for the NR solver algorithm, and Figure 2 shows the flowchart for the control of limits of DGs.

In controlling the DG unit limits, active and reactive power control is carried out in different forms. In case the active power of the ℓth DG is violated, the power of the DG is fixed in its maximum value, i.e., $P_{G\ell }=P_{Gmax\ell }$, and an inductive power factor equal 0.9 is imposed. However, in case a negative value for $P_{G\ell }$ is computed, the value is set for $P_{G\ell }=0$, and the reactive power is $Q_{G\ell }$ in the situation that the reactive power is kept within limits. When the reactive power is violated, $Q_{G\ell }$ is set in its limit.

When the DG unit violates any limit, the generated power is converted into a constant negative load.

The NR method faces difficulties depending on the network loading level and the requirement to control the limits of the DG units in the isolated MG. For this reason, in some situations, the method may fail. Consequently, searching for a robust method for these situations in an isolated network should be encouraged. To improve the performance of the NR solver, a homotopy based-approach using the NR technique is proposed.

3.3. Homotopy-Based Approach for Solving the PFP

Homotopy-based techniques have robust characteristics for solving numerically nonlinear algebraic equation systems as (15) considering limit violations of DGs [49]. The method allows for the transformation of a “difficult” problem into an “easy” one through the manipulation of a parameter t in the parametric space $t\in [0,1]$. The primary objective is to calculate the roots of the nonlinear equations system $\mathbf{g}\left(\mathbf{x}\right)=0$. However, the user needs to introduce an auxiliary function ${\mathbf{g}}_0\left(\mathbf{x}\right)=0$ to form a mixed function $\mathbf{G}(\mathbf{x},t)$. The function ${\mathbf{g}}_0\left(\mathbf{x}\right)=0$ must be chosen in such a way that for $t=0$, the solution is easily found. On the other hand, for $t=1$, the solution of $\mathbf{g}\left(\mathbf{x}\right)=0$ is determined. A good strategy to select ${\mathbf{g}}_0\left(\mathbf{x}\right)=0$ is to establish a practical starting point for this function. The point should stay around the solution, which is searched for the increments of the homotopy parameter until reaching $t=1$.

The modified problem, including t, $\mathbf{g}\left(\mathbf{x}\right)=0$, and ${\mathbf{g}}_0\left(\mathbf{x}\right)=0$, can be presented as the nonlinear set of equations:

$$\mathbf{G}(\mathbf{x},t)=t\mathbf{g}\left(\mathbf{x}\right)+(1-t){\mathbf{g}}_0\left(\mathbf{x}\right)=\mathbf{0},$$

(19)

where $t\in \mathbb{R}$ is the homotopy parameter, and the solution set ${\mathbf{x}}_{*}$ for all $t\in [0,1]$ forms the homotopy trajectory (path).

Now, the focus is on determining the root of $\mathbf{G}(\mathbf{x},t)$ in (19). The definition of ${\mathbf{g}}_0\left(\mathbf{x}\right)$ is crucial to varying the parameter t and calculating the solution of interest in $t=1$. In [27], a homotopy approach was proposed to solve an ill-conditioned problem for a large-scale grid. The method is modified and extended in this paper to consider the isolated MG with constraining in DG units. As in [27], the interest is also to initialize the homotopy problem from a flat start estimation, i.e., all voltages in 1.0 pu (including the reference) and nodal phase angles zero for all buses.

The idea in [27] for $t=0$ is to set ${\mathbf{g}}_0\left(\mathbf{x}\right)$ in such a way that, given an initial estimate ${\mathbf{x}}^{\left(0\right)}$ the solution, ${\mathbf{x}}_{*}^{\left(0\right)}$, for the nonlinear equations system $\mathbf{G}(\mathbf{x},0)=0$ is also ${\mathbf{x}}_{*}^{\left(0\right)}={\mathbf{x}}^{\left(0\right)}$. Therefore, we conclude that ${\mathbf{g}}_0\left(\mathbf{x}\right)=\mathbf{x}-{\mathbf{x}}^{\left(0\right)}$. This means that in $t=0$, we solve a power flow where each bus has its load and a fictitious generation to supply this local load. Hence, no power flow leaves the bus or arrives from other adjacent buses. Similarly, for a DG unit bus, a fictitious load consumes the power generated locally. However, for $t>0$, increasing at discrete step $\Delta t$, such that $t_{i+1}=t_i+\Delta t$, until $t_{i+1}=1$, the fictitious loads and generations are gradually removed according to (19). Obviously, for $t_{i+1}=1$, $\mathbf{G}(\mathbf{x},1)=\mathbf{g}\left(\mathbf{x}\right)=0$ and the PFP solution is determined.

Some remarks concerning the homotopy process can be stated as follows:

• A nonlinear system $\mathbf{G}(\mathbf{x},t_{i+1})=0$ must be solved for each parameter $t_{i+1}$, but the only one whose results are of interest is the one for $t_{i+1}=1$;
• The resolution of $\mathbf{G}(\mathbf{x},t_{i+1})=0$ can also be performed by the standard NR method according to the procedure presented in Figure 1;
• The initial estimate used to solve $\mathbf{G}(\mathbf{x},t_{i+1})=0$ is the solution ${\mathbf{x}}_{*}$ found for $\mathbf{G}(\mathbf{x},t_i)=0$ in the step $t_i$;
• The solution is closer when $\Delta t$ is small, facilitating the convergence of the NR solver;
• As at each step $t_{i+1}$, a NR solver is used and provides a solution for the homotopy path (solutions of $\mathbf{G}(\mathbf{x},t_{i+1})=0$). In the case of DG units, all their limits are verified at the point $t_{i+1}$ before proceeding to the next point. This procedure differs from the case without homotopy (only the classical NR) when the violation condition is only verified for the solution of $\mathbf{g}\left(\mathbf{x}\right)=0$.

Considering the aspects listed for the homotopy-based approach, its application to solve the PFP for isolating MG exploits other qualities not found in the standard NR solver, mainly because it allows one to depart from a desired initial estimate and “guide” the solutions along the homotopy path considering solutions that are initializations for the next points of the path. Additionally, devices submitted to the control of limits of variables are monitored gradually along the homotopy path.

In the next section, experiments are carried out considering some isolated MG and different characteristics of loads and loading cases to demonstrate the numerical performance of the method presented in this paper.

4. Experiments and Results

This section presents simulations and results for demonstrating the efficiency of the homotopy-based approach to PFP in isolated MG. All experiments are performed using a modified version of the MATPOWER tool to incorporate input data for microgrid and load varying with voltage magnitude and frequency. The MATPOWER tool was prepared to run in the OCTAVE platform release 6.4.0. Three isolated MG having only DG units as energy sources supplying the network were studied. The first system is the 6-bus MG studied in [5,18]. The second MG is the 38-bus system presented in [18,46]. The third and largest MG is the 69-bus MG described in [50] and adapted in this work to meet other requirements of loads.

All the initialization of the methods is of the type flat start, i.e., voltage magnitude 1.0 pu and phase angle zero. Step $\Delta t=0.25$ is used when the homotopy-based approach is employed.

The simulation was carried out using a personal computer, Intel(R) Core(TM) i7-8565U 8th Gen, 256 GB SSD CPU, 1.80 GHz, 16 GB RAM, 64-bit OS.

4.1. 6-Bus System Simulations

This subsection describes simulations and results for the 6-bus system. The objective is to show the results determined with the NR method and the homotopy-based approach for a tutorial isolated MG presented in [5,18]. The one-line diagram reproduced from [18] is depicted in Figure 3.

Table 1. 6-bus impedance parameters.

Branches	1–2	1–4	2–3	2–5	3–6
$R_{km}$ ($\mathsf{\Omega }$)	0.300	0.300	0.150	0.200	0.050
$L_{km}$ (mH)	0.318	0.350	1.843	0.250	0.050

shows the parameters $R_{km}$ and $L_{km}$ of the network series branches $k-m$. They were converted to pu for a three-phase power $S_b=10$ kVA and phase-phase voltage base $V_b=220$ V. The loads $L_1$ and $L_3$ were calculated for the nominal voltage, ${\omega }_0=377$ rad/s and impedances $Z_1=6.95+j{\omega }_{0}0.0122\mathsf{\Omega }$ and $Z_3=5.014+j{\omega }_{0}0.0094\mathsf{\Omega }$. Assuming these values and constant impedance, the powers were computed giving $P_{10}=0.4842$ pu, $Q_{10}=0.3204$ pu, $P_{30}=0.6436$ pu, and $Q_{30}=0.4549$ pu. Then, these power values were adopted for simulations considering different parameters $\alpha$ and $\beta$ for load computation as a function of the voltage magnitude and frequency. We kept the parameters $K_{pi}$ and $K_{qi}$ constant for any bus #i, equal to 1 and −1, respectively.

The DG unit basic data for the network simulations were obtained from [5]. The droop data parameters $m_p$ and $n_q$ for the DGs in that reference are given in absolute units, in (rad/s)/W and V/VAr, respectively. We have used $m_p=9.45\times {10}^{-5}$ (rad/s)/W and $n_q=0.0016$ (V/VAr) for all DGs. In our simulations, the data were used in pu. Then, the values $m_p\left(pu\right)=m_p{S}_b/{\omega }_0=0.0025$ pu and $n_q\left(pu\right)=n_q{S}_b/V_b=0.073$ pu. Simulations were performed for all droop types, the NR method, and the homotopy approach. The homotopy was implemented by using the NR method at each point $t_i$ and step $\Delta t=0.25$. For all cases, the convergence was successfully obtained for the tolerance $ϵ={10}^{-8}$ pu for no more than 5 iterations. The DGs were assumed with no power limit.

Table 2. Voltages for the 6-bus system considering different load characteristics and determined by the NR method and homotopy approach.

	$\mathit{\alpha }=0$,  $\mathit{\beta }=0$		$\mathit{\alpha }=1$,  $\mathit{\beta }=1$		$\mathit{\alpha }=2$,  $\mathit{\beta }=2$		$\mathit{\alpha }=0$,  $\mathit{\beta }=2$
Bus	$\mathbf{V}$(pu)	$\mathbf{\theta }$(degree)	$\mathbf{V}$(pu)	$\mathbf{\theta }$(degree)	$\mathbf{V}$(pu)	$\mathbf{\theta }$(degree)	$\mathbf{V}$(pu)	$\mathbf{\theta }$(degree)
#1	0.9564	0.0000	0.9583	0.0000	0.9600	0.0000	0.9582	0.0000
#2	0.9702	−0.5602	0.9714	−0.5401	0.9725	−0.5207	0.9716	−0.5049
#3	0.9609	−2.8716	0.9624	−2.7636	0.9638	−2.6711	0.9631	−2.8156
#4	0.9860	−0.0881	0.9866	−0.0819	0.9872	−0.0737	0.9872	−0.0277
#5	0.9892	−0.4783	0.9896	−0.4629	0.9900	−0.4455	0.9903	−0.3876
#6	0.9668	−3.0697	0.9681	−2.9539	0.9692	−2.8540	0.9689	−2.9953

shows the simulation results for the NR and homotopy approaches. The results are similar, considering the given tolerance. The second and third columns in the table exhibit the voltage magnitude and phase angle for $\alpha =0$ and $\beta =0$. For this result, $\omega =0.9990$ pu was obtained. The fourth to the ninth columns show similar results but for other load characteristics. The fourth and fifth columns exhibit the information for $\alpha =1$ and $\beta =1$, and the value determined for $\omega$ was 0.9991 pu. The sixth and seventh columns depict the results for $\alpha =2$ and $\beta =2$, for which the frequency $\omega =0.9991$ pu was calculated. Finally, the eighth and ninth columns show the results for $\alpha =0$ and $\beta =2$, for which the value $\omega =0.9990$ pu was computed.

From Table 2, it is observed that the NR method and the homotopy approach provided the same results, even for different types of loads. Still, these experiments’ respective values for magnitude voltages and phase angles are close, considering the distinct characteristics of loads.

The tutorial demonstrates the agreement between the NR method and the homotopy-based approach performance.

4.2. 38-Bus System

A larger model based on a 38-bus network was studied in this subsection. The 38-bus system data are based on the network configuration presented in [18], whose one-line diagram is exhibited in Figure 4. The microgrid has 5 DGs of the $P-\omega$ and $Q-V$ droop type. All loads in the system are characterized according to [46] of three types as industrial (I), commercial (C), and residential (R) (see the indication I, C, and R in the one-line diagram). However, in [46], the frequency load dependency is neglected. This dependency is considered in [5] and is adopted in [18] and in this paper.

We emphasize that the feeder connections found in [18,46] differ concerning the configuration detected for the two branches. First, the branch $18-37$ is present in the network in [46], but it is modified to $22–37$ in [18]. Second, the branch $9–35$ in [46] is not present in [18]. Instead, the branch $29–35$ has equal parameters to the branch $9–35$ found in [46]. In view of the explanation, all branch data (line impedance) can be recovered from [46] for a pu base 12.66 kV and 1000 kVA.

The load data can be obtained from [46], including the load type characterized by the exponential parameters $\alpha$ and $\beta$. For load type I, $\alpha =0.18$ and $\beta =6.00$; for load type C, $\alpha =1.51$ and $\beta =3.40$; and for load type R, $\alpha =0.92$ and $\beta =4.04$. The frequency-load dependency parameters were adopted as $K_{pi}=1$ and $K_{qi}=-1$ for the active and reactive part, respectively, such as in [5,18].

DG characterization is considered according to [5,46]. All DGs have inductive droop type and magnitude reference voltage 1.01 pu. The remaining parameters are exhibited in

Table 3. 38-bus system DG parameters.

Bus	${\mathit{m}}_{\mathit{p}}$$\times {10}^{-3}$ (pu)	${\mathit{n}}_{\mathit{q}}$ (pu)	${\mathit{P}}_{\mathbf{max}}$ (MW)	${\mathit{Q}}_{\mathbf{max}}$ (MVar)
#34	5.102	0.020	2	0.9
#35	1.502	0.033	2	0.6
#36	4.506	0.020	2	0.9
#37	2.253	0.050	2	0.3
#38	2.253	0.050	2	0.3

.

Simulations were performed using the NR method proposed in [18] and the homotopy-based approach. The results are in agreement with those found in [18], where the frequency converges to 0.9981 pu, and the reactive power limit of the DG in the bus #38 exceeded its value.

Other experiments were carried out for different loading factors, i.e., modified loads by a factor $\lambda$, such that $P_{Li}=\lambda {\widehat{P}}_{Li}$ and $Q_{Li}=\lambda {\widehat{Q}}_{Li}$, where ${\widehat{P}}_{Li}$ and ${\widehat{Q}}_{Li}$ are the loads for the original case studied in [18]. Progressive values of $\lambda$ were used until the reactive power limits were surpassed or divergence was detected. The last value of $\lambda$ for which we still have the power inside the reactive limits occurred for $\lambda =1.45$.

Table 4. Voltages and power in the DG buses for the 38-bus system.

	$\mathit{\lambda }=1.00$				$\mathit{\lambda }=1.45$
Bus	$\mathbf{V}$(pu)	$\mathbf{\theta }$(degree)	${\mathbf{P}}_{\mathbf{G}}$(pu)	${\mathbf{Q}}_{\mathbf{G}}$(pu)	$\mathbf{V}$(pu)	$\mathbf{\theta }$(degree)	${\mathbf{P}}_{\mathbf{G}}$(pu)	${\mathbf{Q}}_{\mathbf{G}}$(pu)
#34	0.9695	−0.7746	0.3668	0.6773	0.9825	−1.1838	0.5296	0.9000
#35	0.9993	1.2888	1.2461	0.3199	0.9904	1.4674	1.7415	0.5885
#36	0.9971	−1.2086	0.4154	0.6443	0.9839	−1.8753	0.5900	0.9000
#37	0.9973	0.6148	0.8307	0.2536	0.9832	0.9954	1.1992	0.3000
#38	0.9848	0.1868	0.8326	0.3000	0.9650	0.4039	1.1992	0.3000

shows the DG unit generated powers ($P_G$ and $Q_G$) and voltages (magnitude V and phase angle $\theta$) for the loading factors $\lambda =1.00$ and $\lambda =1.45$.

In Table 4, for $\lambda =1.00$, the calculated frequency is $\omega =0.9981$ pu, while for $\lambda =1.45$, the frequency is $\omega =0.9974$ pu. Only the DG unit in the bus #35 for the latter loading kept the power values within its limits.

When $\lambda =1.46$, all DG power reactive limits were surpassed despite the frequency being kept approximately in $\omega =0.9974$ pu.

Again, as it was also verified for the 6-bus system, both the NR method and homotopy-based techniques presented a similar performance in relation to the determination of results. Even with the hard loading and monitoring of the DG unit power limits, both methods converged and presented the same results.

4.3. 69-Bus System

In this section, the objective is to apply the NR method and homotopy-based approach for a higher-order model system with different types of loads and DG units, including loads with year-season profiles. Additionally, it is verified how the two techniques behave with a progressive loading factor.

A modified 69-bus isolated microgrid was generated based on the original radial network presented in [50]. The one-line diagram for the modified network that includes 10 DGs is shown in Figure 5. All data of branches are given in [50] and available in [51] as case69 for pu base 10 MVA and 12.7 kV. The load data (type R, I, or C) were introduced here and characterized by the parameters in

Table 5. Load characteristics informed from tests on the 69-bus Ontario Hydro system [52].

	Summer Day		Summer Night		Winter Day		Winter Night
Load Type	$\mathbf{\alpha }$	$\mathbf{\beta }$	$\mathbf{\alpha }$	$\mathbf{\beta }$	$\mathbf{\alpha }$	$\mathbf{\beta }$	$\mathbf{\alpha }$	$\mathbf{\beta }$
Constant	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
Residential (R)	0.72	2.96	0.92	4.04	1.04	4.19	1.30	4.38
Commercial (C)	1.25	3.50	0.99	3.95	1.50	3.15	1.51	3.40
Industrial (I)	0.18	6.00	0.18	6.00	0.18	6.00	0.18	6.00

. The 10 DG units were also proposed in this paper and are connected at buses #2, #5, #15, #38, #42, #47, #56, #57, #60, and #63. The DG units are assumed to be of the inductive droop type, with each having $m_p=0.011$ pu and $n_q=0.02$ pu, but with different limits for active and reactive power.

Table 6. Maximum limits of active and reactive power of DG units for the 69-bus system.

	#2	#5	#15	#38	#42	#47	#56	#57	#60	#63
P${}_{max}$ (MW)	0.50	0.5	0.5	0.5	0.5	0.5	1.5	1.5	1.5	1.5
Q${}_{max}$ (MVar)	0.25	0.2	0.2	0.2	0.2	0.2	0.9	0.8	0.8	0.8

provides the upper limits of DGs. In the original configuration, all loads are of a constant power type, i.e., with $\alpha =0$ and $\beta =0$. The load profile was modified in this paper to reflect a more realistic characterization of a distribution load characteristic. Then, it was also characterized according to the proposed model presented in [52], as residential (R), commercial (C), and industrial (I) but with parameters changing according to the season’s winter or summer, both for day and night. Table 5 exhibits the parameters associated with these characteristics. The load frequency parameters are equal for all loads, being $K_{pi}=1.0$ pu and $K_{qi}=-1.0$ pu, for any bus #i.

Simulations were performed with different scenarios using the NR method and the homotopy-based approach. Experiments were carried out for loading factors changing from $\lambda =1$ (base case) until divergence was detected. The homotopy approach was performed with step $\Delta t=0.25$.

Table 7. Convergence frequency profile for the modified 69-bus system using the Newton–Raphson method and homotopy-based approach with different load characteristics according to Table 5 and several loading factors.

		Summer Day	Summer Night	Winter Day	Winter Night
Method	$\mathit{\lambda }$	$\mathit{\omega }$
Newton–Raphson	1.00	0.9996	0.9996	0.9996	0.9996
	1.70	0.9988	0.9988	0.9988	0.9988
	1.72	-	0.9993	0.9991	0.9989
	1.73	-	-	-	-
	1.75	-	-	-	-
	1.80	-	-	-	-
Homotopy	1.00	0.9996	0.9996	0.9996	0.9996
	1.70	0.9988	0.9988	0.9988	0.9988
	1.72	0.9990	0.9990	0.9990	0.9990
	1.73	0.9990	0.9990	0.9990	0.9990
	1.75	0.9990	0.9990	0.9990	0.9990
	1.80	-	0.9994	0.9992	0.9990

exhibits the calculated frequency for different load scenarios and loading factors. It was verified that the two techniques applied to the solution of the PFP for the 69-bus system behaved differently for high loads. The NR method could no longer reach a solution for a given load ($\lambda >1.72$). In contrast, the homotopy method continued to show convergence up to $\lambda =1.80$ (divergence occurred for all scenarios with $\lambda \ge 1.81$).

Figure 6 depicts magnitude voltage plots according to the bus for the loads considering the year season and day/night profile for the loading factors for the base case, i.e., ${\lambda }_1=1.00$ and a high loading factor ${\lambda }_2=1.75$. The results obtained with the NR method and homotopy-based approach coincide with the base case. However, only the homotopy-based result is available for a higher loading factor since the NR methods diverged for this condition. A significant difference between the magnitude of voltages from the ninth bus to the twenty-seventh bus is verified by visualizing the curves. This justifies the difficulty for the NR method in converging in high-loading scenarios because the technique uses flat start initialization.

On the other hand, despite using the NR method and departing from the flat start guess, the homotopy-based approach uses intermediate results as initialization for the next NR solver in the homotopy path. This process continues until $t_{i+1}=1$ (the last point of the homotopy path and the one of interest for the PFP solution). Therefore, the initialization of the NR method along the homotopy path provides good conditions for the convergence of the method.

Note that, according to Table 7, for $\lambda >1.75$, the summer day scenario presents a divergence of the power flow. Then, the plots in Figure 6 show the base case plots and those for which the maximum loading factor $\lambda$ is verified for all scenarios. When $\lambda =1.80$, three scenarios still present convergence of the power flow; however, for $\lambda =1.81$, the Jacobian matrix is singular, and divergence occurs for all scenarios.

The CPU time was measured only for the largest system, considering it is more significant among the studied systems. The monitoring was carried out for the base case ($\lambda =1.00$) and a hard loading ($\lambda =1.75$) of the 69-bus system. For the base case, the NR method needed 1.43 s to obtain the PFP solution, while the homotopy method required 2.67 s. The higher time for the homotopy-base approach is justified because it needs to solve more nonlinear systems along the homotopy path to achieve the solution of interest. However, for the hard loading, just the homotopy method had the time measured as it was the only one to converge, requiring 4.55 s. The importance of this result is emphasized, considering that the NR method alone was unable to achieve convergence for the same initialization conditions of the homotopy-based approach.

5. Conclusions

This paper presented an efficient technique to solve the power flow problem in an isolated microgrid network. The approach differs from the conventional power flow problem since instead of using a slack bus, a reference bus is adopted in the isolated microgrid. Additionally, in the latter, frequency is a state that must be determined.

A homotopy-based approach to determine the power flow solution was proposed. The technique considers the isolated microgrid based on droop-controlled DG units and different load characteristics. The PF was always evaluated, departing from a flat start guess and computing intermediate solutions to form the homotopy path until the solution of interest was reached. To compare the performance of the method, the NR solver, also with initialization of type flat start, was used to solve the same problem. Simulations were assessed in 6-, 38-, and 69-bus systems with several DG units, loads with different features, and loading factors.

Experiments for the two smaller systems evidenced the similar performance between the traditional NR method and the homotopy-based approach. However, only the latter technique could reach the solution for the largest system with stringent requirements, such as high loading. This fact is justified since the homotopy approach, despite departing from the flat start guess, uses intermediate results, which are also obtained via the NR method, but uses the solution as the initialization for the next point of the homotopy path until the computation of the last point.

The authors intend to extend the methodology to unbalanced three-phase isolated microgrid studies in future works.