A Robust Electric Spring Model and Modified Backward Forward Solution Method for Microgrids with Distributed Generation

Abstract

The electric spring (ES) is a contemporary device that has emerged as a viable alternative for solving problems associated with voltage and power stability in distributed generation-based smart grids (SG). In order to study the integration of ESs into the electrical network, the steady-state simulation models have been developed as an essential tool. Typically, these models require an equivalent electrical circuit of the in-test networks, which implies adding restrictions for its implementation in simulation software. These restrictions generate simplified models, limiting their application to specific scenarios, which, in some cases, do not fully apply to the needs of modern power systems. Therefore, a robust steady-state model for the ES is proposed in this work to adequately represent the power exchange of multiples ESs in radial micro-grids (µGs) and renewable energy sources regardless of their physical location and without the need of additional restrictions. For solving and controlling the model simulation, a modified backward–forward sweep method (MBFSM) is implemented. In contrast, the voltage control determines the operating conditions of the ESs from the steady-state solution and the reference voltages established for each ES. The model and algorithms of the solution and the control are validated with dynamic simulations. For the quasi-stationary case with distributed renewable generation, the results show an improvement higher than 95% when the ESs are installed. On the other hand, the MBFSM reduces the number of iterations by 34% on average compared to the BFSM.

1. Introduction

The integration of renewable energy sources, the technological evolution of electric vehicles, and the development of more efficient and higher capacity energy storage systems are some alternatives to restrain CO₂ gas emissions and increase sustainable development [1]. To achieve these objectives, the transition from the actual electrical system to a smart grid (SG) requires the improvement of the reliability, security, and efficiency of the electrical system [2]. Furthermore, the migration to a SG involves the flexible handling of the power intermittences associated with the stochastic nature of renewable energy sources (RESs) [3], implementation of monitoring systems, and control schemes that guarantee bidirectional power flow and optimal operation of the SG [4]. An alternative that the electrical industry has been using to modernize the electrical distribution network is the implementation of electrical microgrids (µGs) [5,6]. The µGs are defined as small local distribution systems that operate with the same philosophy of SGs. They promote the distributed generation (DG) and the integration of RESs such as small hydropower units, wind turbines, and photovoltaic arrays; yet, non-RESs such as fuel cells, gas turbines, and combined cycle units are sometimes considered as well [7,8]. Besides, the µGs can be connected to the electrical grid or operate independently (standalone mode) during external faults or disturbances. This versatility increases local reliability, reduces power losses, and improves power quality [9]. One of the significant challenges of integrating DG into the electrical grid is related to the abrupt power variations caused by the weather dependence of RESs. These unexpected power variations can cause congestion in the electrical network when there is a surplus of power or lack of dispatch, affecting the stability of the electrical network [10]. In this sense, the Electric Spring (ES) is a new device that has emerged as a promising alternative to increase the flexibility of the electrical system since it has the particular feature of demand-side load management. This feature allows the ES to manage energy demand locally, increasing the stability of the electrical network and damping both the power oscillations produced by RESs and the dynamic variations of loads [11].

In recent years, the interest of the scientific community to study the impact of ESs on the electrical network has grown. These studies have been oriented to the design and implementation of collaborative control schemes that improve voltage regulation [12,13,14], keep the frequency within safe operating ranges [15,16], and operate the electrical systems in optimal conditions to reduce power losses [17,18]. These situations are studied using simulation models of radial µGs or radial distribution systems, which are divided into two main groups: time-domain simulation models and steady-state simulation models. Time-domain simulations are implemented in specialized electromagnetic transient software and are appropriate when it is required to analyze the behavior of the network during transients. However, they have as a drawback the requirement of a tremendous computational effort which, under certain conditions, limits the analysis of large-scale electric systems [19,20].

On the other hand, steady-state simulation models are widely used to analyze large-scale electric systems, and they are suitable to evaluate the slow dynamic of electromechanical systems using phasor solution methods. The steady-state analyses for the ES that haven been reported in the literature carry out different studies to understand the behavior and interaction of ES with the electrical network. For instance, the steady-state analysis performed in [21] shows experimentally the operating conditions of the ES that determine the exchange of active and reactive power with the electrical network. In [22], the operating range of the ES is analyzed according to its parameters, and a control scheme is established to minimize the action of the ES for regulating the voltage in the point of common coupling (PCC). The design of a dynamic controller to regulate the active and reactive power exchange between the ES and the electrical network is proposed in [23]; for the controller design, the authors carried out the steady-state analysis of a simplified electrical network and a single ES. From the steady-state model, the geometric relationships between the electrical variables involved were derived. As a result, a control law that depends on the parameters of the electrical network and the electrical variables at the PCC is obtained. In [24], a steady-state model to describe the exchange of active and reactive powers in the ES is also obtained. The simulation model is implemented with a test network and its impact is evaluated when it operates in conjunction with the RES and electric vehicles. In [25], a steady-state model of the ES and its solution and control algorithms for the implementation of study cases with multiple ES are proposed. The presented model decouples the SL into a virtual generator and an equivalent impedance connected to the PCC. The exchange of power with the electrical network for both devices is modeled from their power equations. The virtual generator represented by a voltage source is used to achieve the power balance in the PCC associated with the operation of the SL. This task is performed iteratively using a voltage control algorithm and the BFSM. These algorithms have the particularity that they are closely linked and, therefore, it is not possible to decouple them, which limits the implementation of voltage control schemes. Other limitations are that the study cases presented have been tested with purely resistive critical loads and include restrictions on the voltage of the critical loads to a value of 0.7 p.u. Despite obtaining promising results, the above-mentioned steady-state models were developed to describe the exchange of active and reactive power of the ES with the electrical network by using an equivalent circuit (Thevenin equivalent), which requires the addition of restrictions such as the parameters of the electrical network, the physical location of a single ES, the power factor of the loads, and/or the operation mode of the ES. These restrictions generate simpler models in terms of computational cost or mathematical complexity but limit their applicability to more general and real scenarios as the ones presented in the modern power systems. In this regard, the development of more robust ES models that allow reducing for the number of restrictions, e.g., the implementation of multiple ES distributed along a µG, is still of paramount importance for satisfying the new requirements of modern electrical networks, e.g., the implementation of collaborative control schemes.

In this work, a new steady-state model for the ES is developed. In general, the proposed model allows for considering any location of the ES into the µG without the need to include additional restrictions and describing the power exchange with the µG regardless of (i) its mode of operation, (ii) the distribution line parameters of the µG, and (iii) the presence of external disturbances associated with load variations, injection of active power of REs, and external voltage variations. All these features increase the robustness of the proposed model. On the other hand, the algorithms for the solution and control of the µG with the ESs are also presented. The steady-state solution of the µG is calculated through the proposal of a modified backward–forward sweep method (MBFSM); due to the ES model design and the MBFSM, it was possible to omit the forward propagation stage, making the solution process simpler. Regarding the voltage control, the proposed algorithm operates as a global control to determine the set point of each ES considering its operating limits. The implementation for the ES model, solution method, and control algorithm is validated by several study cases and compared with dynamic simulations, demonstrating the usefulness and effectiveness of the proposal.

2. Theoretical Background for the Steady-State ES Model

The first version of the ES was presented in 2012 [11]; from this version, different topology modifications have been presented over the years. In this sense, the development of steady-state ES models has been also presented for its analysis in power systems [26,27,28]. This section presents the basis for the modelling of ESs and their operating principle in the time-domain and frequency-domain.

2.1. ES Time-Domain Model

Figure 1 shows the ES dynamic model, which includes its controller and power stage. The power stage of the ES incorporates a DC source, V_DC, a DC–AC converter, and a low pass filter integrated by an inductor-capacitor (L_f and C_f). The converter operates as a single-phase full-bridge inverter consisting of four switches (S₁–S₄), which are commutated using the pulse width modulation (PWM) technique. The low pass filter is in charge of filtering high frequency components and providing a sinusoidal output voltage, V_ES, at the fundamental frequency [29]. The control stage contains two controllers that work together to establish the operating conditions of the ES, i.e., the inductive or capacitive operation mode. The non-critical current, I_NC, and the voltage on the PPC, Vs, are fed back to the controller. A reference voltage, Vs_{_}ref, is established, which is compared with Vs to obtain the controller error, e(t). Furthermore, an integral proportional (PI) controller determines the magnitude and operation mode of the ES [29]. The PI controller output is the magnitude of the control voltage, V_{ES_mag}. A phase-locked loop (PLL) is responsible for calculating the phase angle of the control voltage, V_{ES_phase}, which is delayed by 90 degrees with respect to the I_NC for reactive compensation mode. The capacitive or inductive mode is determined by the sign of V_{ES_mag}. If the sign is positive, I_NC leads V_{ES_ref} by 90 degrees, causing the ES to operate in capacitive mode and, consequently, increasing Vs; otherwise, if the sign is negative, the ES operates in inductive mode, decreasing Vs [29].

The configuration of the non-critical load, Z_NC, and the ES is known as a smart load (SL). The SL is capable of controlling the active and reactive power exchanged with the electrical grid, even if the ES only operates in reactive compensation mode. The modification of the V_ES in both modes of operation has repercussions on the voltage applied to the non-critical load, V_NC. Therefore, the power consumed by the non-critical load can be controlled and the ES simultaneously regulates the voltage at the PCC, managing any intermittence in the RES locally [30].

2.2. ES Frequency-Domain Model

The frequency-domain modeling involves the development of new ES models that adequately and precisely describe the power exchange of the ES with the electrical grid. Unlike time-domain models, in the frequency-domain models, it is not possible to implement controllers that automatically respond to the disturbances and adjust the ES operating point.; therefore, to derivate a steady-state model, the power exchange and the operating principle must be included. These models are developed using phasor representations.

The ES models reported in the literature are obtained from simplified electrical circuits, such as the one shown in Figure 2a [21,24]. The electrical network is represented by a voltage source, ${\overrightarrow{V}}_G$ and its equivalent impedance, Z_L. The ES is modeled as a controlled voltage source ${\overrightarrow{V}}_{ES}$, which is connected in series between the non-critical load, Z_NC, and the PCC $\overrightarrow{V}s$. These models describe the power exchange of the ES with the electrical network, but they have the limitation that they are valid only when the SL is grid-tied at the farthest node of the electrical network [31]. In order to implement this model in simulation software, it is necessary to add additional restrictions that depend on the ES operation mode, the power factor of the loads, and the particular conditions of any study case.

Figure 2b shows the phasor diagram on the PCC of Figure 2a. In the analysis, it is considered that the non-critical load has a lagging power factor, and the ES operates in capacitive operation mode. The geometric relationships between the phasors $\overrightarrow{V}s$, ${\overrightarrow{V}}_{NC}$, ${\overrightarrow{V}}_{ES}$, ${\overrightarrow{I}}_{NC}$, and the non-critical load angle, ϕ_NC, is given by [31]:

$$V{s}^2={\left(V_{NC}-V_{ES}\sin {\varphi }_{NC}\right)}^2+{\left(V_{ES}\cos {\varphi }_{NC}\right)}^2.$$

(1)

By solving Equation (1) for V_NC, the following expression is obtained.

$$V_{NC}=\pm \sqrt{V{s}^2-{\left(V_{ES}\cos {\varphi }_{NC}\right)}^2}+V_{ES}\sin {\varphi }_{NC}.$$

(2)

Once V_NC is obtained, the next step is to calculate the reactive powers of the ES, Q_ES, and the non-critical load, Q_NC. These are given in Equations (3) and (4), respectively

$$Q_{ES}=V_{ES}{I}_{NC}\sin \left(-90°\right)=-V_{ES}{I}_{NC}=-\frac{V_{ES}{V}_{NC}}{Z_{NC}},$$

(3)

$$Q_{NC}=V_{NC}{I}_{NC}\sin \left({\varphi }_{NC}\right)=\frac{V_{NC}^2}{Z_{NC}}\sin \left({\varphi }_{NC}\right).$$

(4)

The expression for reactive power of the smart load (Q_SL) is given below

$$Q_{SL}=Q_{ES}+Q_{NC}.$$

(5)

Substituting (2)–(4) in (5), the Equation (6) is obtained

$$Q_{SL}=\frac{-V_{ES}\left(\pm \sqrt{V{s}^2-{\left(V_{ES}\cos {\varphi }_{NC}\right)}^2}+V_{ES}\sin {\varphi }_{NC}\right)}{Z_{NC}}+\frac{{\left(\pm \sqrt{V{s}^2-{\left(V_{ES}\cos {\varphi }_{NC}\right)}^2}+V_{ES}\sin {\varphi }_{NC}\right)}^2}{Z_{NC}}\sin \left({\varphi }_{NC}\right).$$

(6)

Equation (6) represents the steady-state model associated with the reactive power exchange of the SL with the electrical network. It can be seen that the model depends on the magnitudes of Vs, V_ES, as well as the non-critical impedance, Z_NC, and its power factor. This equation is valid only when the ES operates in capacitive mode and the power factor of the non-critical is lagging. If the operation mode changes or the critical load is modified, the equation would not be applicable and, consequently, it would be necessary to derive a new equation for those particular operating conditions. Therefore, it is necessary to develop new mathematical models of the ES that adequately represent its behavior, including its different operation modes and considering external disturbances. These considerations will allow for studying the impact of multiple distributed ESs in a µG and the inclusion of collaborative controls in DG to improve the reliability of the µGs.

3. Model Obtaining

In order to obtain the proposed steady-state model, the simulation diagram shown in Figure 3 is proposed. It involves DG and the collaborative operation of multiple ESs distributed into the µG. Each ESs is modeled as a voltage source and controlled by its local control (see Figure 3b). This feature provides flexibility to the model since it is possible to analyze the individual operation of the ESs or establish a collaborative control scheme. On the other hand, the DG is modeled using a set of current sources distributed in the µG (see Figure 3c). These current sources inject active power and are used to represent the stochastic behavior of the RES. For validation purposes, the steady-state and time-domain models of the µG are implemented under the same operating conditions.

The steady-state model consists of three main stages, which are the µG, the MBFSM algorithm, and the voltage control algorithm, as shown in Figure 3a. The µG has a radial topology, which is one of the most used in distribution systems, and this is composed of a voltage source, V_G, that represents the utility grid and a set of elements connected through the line impedances (Z_L1, Z_L2, …, Z_Ln). The voltages of the µG are numbered sequentially, starting from the closest one to V_G until reaching the last node (Vs₁, Vs₂, …, Vs_n). The physical location of an element is given by the subscript number of the node to which it is connected. These two actions allow for the homogenizing of the information and for the identifying of the elements that compound the µG. The µG contains critical loads (Z₁, Z₂, …, Z_n) and non-critical loads (Z_NC1, Z_NC2, …, Z_NCn). The critical loads are shunt connected with the SLs, and non-critical loads are series connected with the ESs. The DG is represented employing a set of current sources (I_REN1, I_REN2, …, I_RENn), where the subscript REN means renewable energies. The ESs are modeled as a controlled voltage sources (V_ES1, V_ES2, …, V_ESn) which modify their voltages based on their control variables (m₁, m₂, …, m_n).

The steady-state solution of the µG is obtained using the MBFSM algorithm. The algorithm receives as input the operating conditions of the ESs (m₁, m₂, …, m_n), and the active power injected by each of the RESs (P_REN1, P_REN2, …, P_RENn). In each iteration, the currents delivered by the current sources are calculated by Equation (7) to inject active power only.

$$I_{RENx}=\frac{P_{RENx}}{V{s}_x},$$

(7)

where the subscript x is associated with the physical location of the current source and takes the values x = 1, 2, …, n, when the MBFWM algorithm converges and the steady-state solution is obtained. An additional algorithm is in charge of controlling the voltage of each ES. The voltage control algorithm receives as input the steady-state solution calculated by the MBFSM and adjusts the operating conditions of the ESs until reaching the reference voltages (Vs_ref₁, Vs_ref₂, …, Vs_ref_n).

The validation of the proposed model is carried out by comparing the obtained results with the ones obtained by the solution in the steady-state condition of the time-domain-based model for the µG under the same operating conditions. For the time-domain simulation, the ES is implemented by using its dynamic model (see Figure 3b). At the same time, the current sources associated with the DG are replaced by the dynamic model depicted in Figure 3c. Each current source consists of a DC-AC converter and a low-pass filter (L_f1-C_f-L_f2), operating as a controlled current source. The magnitude and phase for the current signal are calculated from the delivered active power and the voltage in the PCC.

3.1. Steady-State ES Mathematical Model

The dotted rectangle in Figure 4 contains the elements and variables involved in the derivation of the ES mathematical model. This part corresponds to a general section of the circuit shown in Figure 3a. There are three nodes involved in obtaining the model, which are indicated by different subscripts (i, j, and k). The subscript j represents the physical location of the ES, while i and k are used to indicate the neighboring nodes. The current contributions of all the elements connected to the PCC, $\overrightarrow{V}{s}_j$, are considered. It can be seen that the components connected to $\overrightarrow{V}{s}_j$ are those that are usually included in the steady-state models, e.g., the critical loads Z_j and the SL. To increase the robustness of the model, the injection of active power into the PCC, ${\overrightarrow{I}}_{RENj}$, and the effect of external disturbances associated with voltage variations, $\overrightarrow{V}{s}_i$, and current fluctuations ${\overrightarrow{I}}_{jk}$ are included. $\overrightarrow{V}{s}_i$ in the model involves the voltage variations attributed to ${\overrightarrow{V}}_G$ or the power variations at nodes that precede the location of the ES, and the impact on its operation. ${\overrightarrow{I}}_{jk}$ corresponds to the current demanded by either the loads or other elements connected in the µG. On the other hand, when ${\overrightarrow{I}}_{jk}$ is included, the effect of power variations on the rest of the network and the impact that they have on the operation of the ES are also included in the model.

A steady-state analysis is performed on the elements found within the dotted rectangle in Figure 4 to obtain the ESs mathematical model. By applying Kirchhoff’s Voltage Law (KVL) between $\overrightarrow{V}{s}_i$ and $\overrightarrow{V}{s}_j$, the following equation is obtained.

$$\overrightarrow{V}{s}_i={\overrightarrow{I}}_{ij}{Z}_{Lj}+\overrightarrow{V}{s}_j$$

(8)

where ${\overrightarrow{I}}_{ij}$ is the current that flows through the distribution line impedance, Z_Lj. Applying a Kirchhoff’s Current Law (KLC) at the node $\overrightarrow{V}{s}_j$, the following expression is obtained,

$${\overrightarrow{I}}_{ij}={\overrightarrow{I}}_j+{\overrightarrow{I}}_{N{C}_j}-{\overrightarrow{I}}_{RENj}+{\overrightarrow{I}}_{jk},$$

(9)

where ${\overrightarrow{I}}_j$, ${\overrightarrow{I}}_{NCj}$, …, and ${\overrightarrow{I}}_{jk}$ correspond to the currents in the critical load, non-critical load, the current delivered by the RESs, and the current demanded by other elements of the µG. The current ${\overrightarrow{I}}_j$ is expressed in terms of the voltage $\overrightarrow{V}{s}_j$ and the impedance Z_j as follows:

$${\overrightarrow{I}}_j=\frac{\overrightarrow{V}{s}_j}{Z_j}.$$

(10)

By substituting Equations (9) and (10) in Equation (8), it is obtained that:

$$\overrightarrow{V}{s}_i=\overrightarrow{V}{s}_j\left[\frac{Z_{Lj}}{Z_j}+1\right]+Z_{Lj}{\overrightarrow{I}}_{NCj}-Z_{Lj}{\overrightarrow{I}}_{RENj}+Z_{Lj}{\overrightarrow{I}}_{jk}.$$

(11)

Expressing the active power delivered for the renewable energy source, P_RENj, in terms of $\overrightarrow{V}{s}_j$ and ${\left({\overrightarrow{I}}_{RENj}\right)}^{*}$ in polar form, it results that:

$$P_{RENj}=\overrightarrow{V}{s}_j{\left({\overrightarrow{I}}_{RENj}\right)}^{*}=\left(V{s}_j\overline{)\theta s_j}\right)\left(I_{RENj}\overline{)-\theta s_j}\right),$$

(12)

where θs_j is the phase angle of $\overrightarrow{V}{s}_j$ and ${\left({\overrightarrow{I}}_{RENj}\right)}^{*}$ is the conjugated of the current ${\overrightarrow{I}}_{RENj}$.

From Equation (12), ${\overrightarrow{I}}_{RENj}$ can be represented in terms of the phasor $\overrightarrow{V}{s}_j$ as given in Equation (13)

$${\overrightarrow{I}}_{RENj}=I_{RENj}\overline{)\theta s_j}=\frac{P_{REN}}{V{s}_j^2}\overrightarrow{V}{s}_j.$$

(13)

Equations (14) and (15) give the $\overrightarrow{V}{s}_j$ and ${\overrightarrow{I}}_{jk}$ phasors represented in polar form, respectively.

$$\overrightarrow{V}{s}_j=V{s}_j\overline{){\theta }_{sj}},$$

(14)

$${\overrightarrow{I}}_{jk}=I_{jk}\overline{){\theta }_{jk}},$$

(15)

where θ_jk is the phase angle of the current ${\overrightarrow{I}}_{jk}$.

Defining α = θs_j − θ_jk, then ${\overrightarrow{I}}_{jk}$ can be expressed in terms of $\overrightarrow{V}{s}_j$ as shown in Equation (16)

$${\overrightarrow{I}}_{jk}=\left(\frac{I_{jk}}{V{s}_j}\overline{){\theta }_{jk}-\theta s_j}\right)\left(V{s}_j\overline{)\theta s_j}\right)=\left(\frac{I_{jk}}{V{s}_j}\overline{)-}\alpha \right)\overrightarrow{V}{s}_j.$$

(16)

Substituting Equations (13) and (16) in Equation (11), it results:

$$\overrightarrow{V}{s}_i=\overrightarrow{V}{s}_j\left[1+\frac{Z_{Lj}}{Z_j}+\frac{Z_{Lj}{I}_{jk}}{V{s}_j}\overline{)-}\alpha -\frac{Z_{Lj}{P}_{REN}}{V{s}_j^2}\right]+Z_{Lj}{\overrightarrow{I}}_{NCj}.$$

(17)

Now, defining ${\overrightarrow{V}}_{G1}$ as follows:

$${\overrightarrow{V}}_{G1}=\overrightarrow{V}{s}_j\left[1+\frac{Z_{Lj}}{Z_j}+\frac{Z_{Lj}{I}_{jk}}{V{s}_j}\overline{)-}\alpha -\frac{Z_{Lj}{P}_{REN}}{V{s}_j^2}\right],$$

(18)

equations can be expressed in a compact form by:

$$\overrightarrow{V}{s}_i={\overrightarrow{V}}_{G1}+Z_{Lj}{\overrightarrow{I}}_{NCj}.$$

(19)

Additionally, by applying a KVL in the $\overrightarrow{V}{s}_j$ node, the equation associated with the SL is given by:

$$\overrightarrow{V}{s}_j={\overrightarrow{V}}_{ESj}+{\overrightarrow{I}}_{NCj}{Z}_{NCj}.$$

(20)

Unlike other works, the proposed model and the calculation of geometric relationships among the electrical variables of the ES consider the load contributions in the PCC and the current contributions associated with external disturbances, as was shown in the approach to derive Equations (19) and (20). This model, unlike the one reported in [23], is not obtained from an equivalent circuit of the electrical network, which increases its robustness. The geometric relationships of the model described in Equations (19) and (20) are illustrated in Figure 5a,b. The phasors represent the Equation (19) in black color. This equation considers the voltage drop in Z_Lj produced by the current ${\overrightarrow{I}}_{NCj}$. This voltage drop is divided into two components related to the resistance, R_Lj, and reactance, X_Lj of Z_Lj. The other term involves ${\overrightarrow{V}}_{G1}$, whose magnitude and phase are influenced by different parameters, as shown in Equation (18).

On the other hand, Equation (20) is represented for the blue color phasors in Figure 5a,b. One term represents the voltage drop across the non-critical load ${\overrightarrow{V}}_{NCj}$, having two components associated with the resistance, R_NCj, and reactance, X_NCj of Z_NCj, and the other term is linked with the ES operation mode. In inductive operation mode ${\overrightarrow{V}}_{ESj}$ is lagging ${\overrightarrow{I}}_{NCj}$ by 90^o degrees, as shown in Figure 5a. In capacitive operation mode, ${\overrightarrow{V}}_{ESj}$ it is leading ${\overrightarrow{I}}_{NCj}$ by 90^o, as shown in Figure 5b. The geometric relationships involve four angles. The angle ${\phi }_1$ is measured from $\overrightarrow{V}{s}_j$ to ${\overrightarrow{V}}_{G1}$, and the angle ${\phi }_0$ is associated with the power factor of Z_Lj. The angle ${\phi }_4$ is measured from ${\overrightarrow{V}}_{G1}$ to $\overrightarrow{V}{s}_i$, and finally, $\theta /2$ is the angle between ${\overrightarrow{I}}_{NCj}$ and $\overrightarrow{V}{s}_j$.

When performing the analysis of the geometric relationships of the phasor diagrams in Figure 5,b, the mathematical model described by the next Equations (21)–(29) is obtained. This model differs from that reported in [23] in the calculation of V_G1 and ${\phi }_1$, while the rest of the geometric relationships remain unchanged. This fact is closely linked with the procedure to derive Equations (13) and (16), which allowed finding a connection between ${\overrightarrow{V}}_{G1}$ and $\overrightarrow{V}{s}_j$. As a result, a model like the one expressed in Equations (19) and (20) is obtained.

It can be seen in Equations (21) and (22) that the magnitude of V_G1 and ${\phi }_1$ depend on the impedances connected to the PCC (Z_Lj and Z_j), the active power, P_RENj, the magnitude of the current I_jk, the magnitude of the voltage on the PCC, V_sj, and the angle α. Variations of any of these parameters will influence the calculation of their geometric relationships. Therefore, the modeling and implementation of the µG in simulation is not a trivial task. However, the BFSM provides an advantage in the implementation of the ES mathematical model. As previously mentioned, this solution technique starts at the node farthest from the µG and, employing the backward propagation, the voltages and currents are calculated until reaching the power supply. This advantage provided by BFSM allows the ES model to be implemented in any µG location without the need for additional restrictions. For instance, to calculate the geometric relationships of an ES connected to the PCC at node Vs_j, the BFSM provides the values of $\overrightarrow{V}{s}_j$ and ${\overrightarrow{I}}_{jk}$. The rest of the terms are parameters and operating conditions of the PCC, which are known.

$$V_{G1}=\left|1+\frac{Z_{Lj}}{Z_j}+\frac{Z_{Lj}{I}_{jk}}{V{s}_j}\overline{)-}\alpha -\frac{Z_{Lj}{P}_{RENj}}{V{s}_j^2}\right|V{s}_j.$$

(21)

$${\phi }_1=\arg \left(1+\frac{Z_{Lj}}{Z_j}+\frac{Z_{Lj}{I}_{jk}}{V{s}_j}\overline{)-}\alpha -\frac{Z_{Lj}{P}_{RENj}}{V{s}_j^2}\right).$$

(22)

The angle ${\phi }_0$ depends on the power factor of Z_Lj as follows:

$${\phi }_0=\arg \left(Z_{Lj}\right).$$

(23)

The ES operating conditions are defined by the parameter m_j, and this parameter takes values in the range of [−1, 1] and determines the magnitude of Vs_i and the angle θ according to the following relationships.

$$V{s}_i=\sqrt{m_j\sqrt{a^2+b^2}+V_{G1}^2+b}.$$

(24)

$$\theta ={\sin }^{-1}\left(m_j\right)-{\tan }^{-1}\left(\frac{b}{a}\right).$$

(25)

The parameters a and b are used to simplify Equations (24) and (25). They are given by:

$$a=\frac{V_{G1}V{s}_j}{R_{NCj}}\sqrt{R_{Lj}^2+X_{Lj}^2}\cos \left({\phi }_0+{\phi }_1\right).$$

(26)

$$b=\frac{V{s}_j^2\left(R_{Lj}^2+X_{Lj}^2\right)}{2R_{NCj}^2}+\frac{V_{G1}V{s}_j}{R_{NCj}}\sqrt{R_{Lj}^2+X_{Lj}^2}\sin \left({\phi }_0+{\phi }_1\right).$$

(27)

The magnitude of the ES voltage is calculated by Equation (28). If the non-critical load has a unity power factor, the second term of the equation is omitted. If it has a lagging power factor, the sign is negative, and if it is leading, the sign is positive.

$$V_{ESj}=V{s}_j\sin \left(\frac{\theta }{2}\right)\pm I_3{X}_{NCj}.$$

(28)

The angle ${\phi }_4$ is calculated by Equation (29)

$${\phi }_4={\sin }^{-1}\left(\frac{V{s}_j\sqrt{R_j^2+X_j^2}}{R_{NCj}V{s}_i}\cos \left(\frac{\theta }{2}\right)\cos \left(\frac{\theta }{2}+{\phi }_0+{\phi }_1\right)\right).$$

(29)

It can be seen in Figure 5a,b that the geometric relationships take as reference the phasor ${\overrightarrow{I}}_{NCj}$. The reference was established to simplify the calculation of geometric relationships. However, it is necessary to calculate the exact value of the phase angle for the ES voltage, θ_ESj. This value is obtained by the Equation (30) and is calculated from θ_Sj. The value of θ_Sj is provided by the MBFSM in the backward propagation process. If the angle θ takes a positive value, the ES operates in inductive mode; while for a negative value, the operating mode is capacitive.

$${\theta }_{ESj}=\theta s_j-\frac{\theta }{2}+sign\left(\theta \right)\frac{\pi }{2}.$$

(30)

In the same way, the phase angle for the voltage Vs_i, θs_i is calculated taking as reference θ_Sj in the next equation.

$$\theta s_i=\theta s_j+{\phi }_1+{\phi }_4.$$

(31)

Finally, to calculate the current and continue with the backward propagation process, Equation (9) is used. The geometric relationships described in Equations (21)–(31), along with the MBFSM solution technique, allow analyzing the joint operation of RES and ESs distributed in the µG. MBFSM is an iterative process that calculates the steady-state solution of the µG until it converges or reaches a maximum number of iterations. The following subsection describes in detail the implementation of the solution algorithm and the proposed modification to the traditional BFSM.

3.2. Modified Backward–Forward Sweep Method

The BFSM has been extensively used for the steady-state solution of radial electric networks due to its accuracy and fast convergence. There are three variants in the method that depends directly on the type of electrical quantities used: (a) the current summation method, (b) the power summation method, and (c) the admittance summation method [32]. The BFSM procedure is depicted in Figure 6a. It starts with a flat solution in which the voltage at all nodes takes the same value. The backward sweep begins at the farthest node of the radial grid and moves toward the source node computing the currents or power flow in each branch. The forward sweep computes the voltage drop in the electrical grid from the currents or power flows calculated in the backward sweep. The nodal voltages are updated from the voltage source to the last node of the electric grid. During the forward sweep process, the currents or powers calculated in the backward sweep remains constant. The values calculated in the forward sweep correspond to the new voltage values used in the backward sweep for the next iteration. The convergence is reached when the mismatch voltage is less than the specified tolerance.

In this work, it is proposed to modify the traditional BFSM algorithm eliminating the forward propagation process. In general, two factors benefit the proposal and make its implementation viable. The first one is that in the backward propagation, the algorithm simultaneously calculates the currents and voltages at each node. The second one is associated with the proposed ES model, which was designed to operate in conjunction with the solution method. As was seen in the previous section, the model calculates the magnitudes and phase angles at both the PCC and the node that precedes it from the ES operating condition defined by the parameter m_j. Therefore, it is possible to simultaneously calculate the currents and voltages at the nodes where the ESs are physically located, avoiding the need to implement the forward propagation step.

Figure 6b shows the flowchart for the MBFSM algorithm proposed in this work. This algorithm computes the solution for the proposed model (see Figure 3a). The algorithm begins with the reading of the parameters (lines, load, and ESs) that conform to a µG. Then, the active power delivered by the current sources (P_REN1, P_REN2, …, P_RENx) and the operation conditions of the ESs (m₁, m₂, …, m_y) are established. It is worth noting that different subscripts are used in order to generate a complete scenario; for instance, the subscripts x and y are used to indicate that the number of ESs and RES located in the µG may be different and not necessarily be equal to the number of nodes in the µG established by the subscript n. The initial condition of the farthest voltage in the µG is assigned, Vs_n. It can be seen from the flowchart in Figure 6b that forward propagation is omitted because, in the backward propagation, the currents and node voltages are simultaneously calculated until reaching the power supply, V_G. If a node, where an ES is located, is reached during the backward propagation process, the model presented in Equations (21)–(31) is employed to calculate the geometric relationships of the electrical quantities in the PCC. As the voltage V_G calculated in the backward propagation may not correspond to the real voltage of the power supply, V_Gn, the output error, ErrSal, is calculated as the absolute value of the difference between V_G and V_Gn. If ErrSal is less than a tolerance, tol, the algorithm stops. The value assigned to tol directly impacts the accuracy of the solution and execution time. The tuning of this parameter must be carefully done as it depends on the dimensions of the µG and the number of ESs and RESs. In this work, a value of tol = 1 × 10⁻¹² was used, which was determined through an exhaustive experimentation. If the condition is not met, a new value of Vs_n is calculated for the next iteration. The new value of Vs_n is calculated by dividing the current value by an adjustment gain, K. The gain K is calculated by the ratio of the absolute value of V_G and V_Gn. The process is repeated until the maximum number of iterations, N_iter, is reached, or the algorithm converges.

The MBFSM accurately provides the steady-state solution of the µG. This solution is calculated from the operating conditions of the ESs and the active power injected by each of the RES; therefore, it is possible to study the impact that DG has on µG. Another area of opportunity for the proposal is provided by the ability to establish the operating conditions of the ESs. These features allow analyzing the individual or joint operation of the ESs and developing global control schemes that face the challenges of modern electric power networks. In this sense, the following subsection explains in detail the voltage control scheme proposed in this work, which is in charge of regulating the operating conditions of multiple ESs distributed in the µG.

3.3. Voltage Control Algorithm Applied to Multiple ESs

Figure 7 depicts the flowchart for the voltage control algorithm of a group of ESs distributed in the electric grid. This algorithm controls the interaction of all the elements in Figure 3a. The algorithm starts by initializing the operation conditions of all the ESs (m₁, m₂, …, m_y) and their respective reference voltage (Vs_{_ref1}, Vs_{_ref2},…, Vs__refy). The magnitude of the active powers supplied for the current sources are established, which are kept constant until the algorithm converges (P_REN1, P_REN2, …, P_RENx). The algorithm computes the steady-state solution of the µG in each iteration using the MBFSM. Only the voltages of the nodes where the ESs are located are extracted from the steady-state solution (Vs₁, Vs₂, …, Vs_y).

The root mean square error, Err, is calculated from the Vs_β and Vs_ref_β, where β =1,2, …, y. If Err is less than the tolerance, tolc, or the voltage magnitude of all ESs reaches its maximum value, V_max, then the algorithm stops. V_max is directly related to the magnitude of the voltage source on the DC bus. In one case of the two previous conditions is not fulfilled, it is necessary to modify the operating conditions of each one of the ESs. The following procedure applies to each situation individually. The ESs that reached their V_max are identified, and their m_β are not modified in this iteration. For the ESs that did not reach V_max, the sign is obtained from the subtraction (Vs_β − Vs_ref_β). If the sign is positive, it means that the voltage Vs_β is higher than Vs_ref_β; therefore, it is necessary to consume a higher amount of reactive power to decrease the voltage at node Vs_β, which is achieved by increasing m_β, having as a limit a maximum value of 1. Otherwise, m_β must be reduced to a minimum value of −1. This action is associated with an increment in the reactive power injected by the ES. The process is repeated until the solution converges or the maximum number of iterations has been reached.

Once the µG voltage control and the solution algorithms have been implemented, it is necessary to validate the ES model under different operating conditions. In the next section, the obtained results for the proposal are presented.

4. Study Cases and Results

Figure 8 shows the µG implemented to carry out the model validation. It can be seen that it has the topology shown in Figure 3. The assignment of subscripts in the elements corresponds to the criteria described in Section 3. The µG is composed of a power supply, V_G, five distribution lines (Z_L1, Z_L2, Z_L3, Z_L4, and Z_L5), three ESs (ES2, ES4, and ES5), two RES (I_REN2 and I_REN5) with powers (P_REN2 and P_REN5), three non-critical loads (Z_NC2, Z_NC4, and Z_NC5), and five loads (Z₁, Z₂, Z₃, Z₄, and Z₅), where Z₂, Z₄, and Z₅ are the critical loads. The five node voltages are indicated as (Vs₁, Vs₂, Vs₃, Vs₄, and Vs₅).

Table 1. Parameters of the single-phase µG.

Distribution Lines			Loads			Non-Critical Loads		DC Bus for ESs
Element	Resistance (Ω)	Reactance (Ω)	Element	Resistance (Ω)	Reactance (Ω)	Element	Resistance (Ω)	Element	Voltage (V)
ZL1	0.1	0.4599	Z1	0	−70	ZNC2	50	VDC2	100
ZL2	0.1	0.9048	Z2	53	0	ZNC4	50	VDC4	100
ZL3	0.1	0.9048	Z3	53	50	ZNC5	50	VDC5	100
ZL4	0.1	0.4599	Z4	0	−70
ZL5	0.1	0.4599	Z5	53	50

shows the numerical values for the elements that integrate the µG, and these parameters are similar to the ones used in other works to study the impact of the ESs in the µGs [33].

For the implementation of the dynamic simulation model, each ES has its control and has no communication with other ESs. The control has a PLL to operate in reactive power compensation mode and a PI controller to regulate the magnitude and operation mode of the ESs (see Figure 3b). For each I_REN, a PLL is used to maintain the injected current in phase with the voltage to deliver active power only (see Figure 3c). The controllers’ gains were adjusted experimentally using the best transient response as the fit criterion. The obtained gains are shown in

Table 2. Gains of the PI controllers and phase-locked loops (PLLs).

	Sampling Time	Proportional Gain (Kp)	Integral Gain (Ki)	Integral Gain (Kd)
PI Controllers	0.0167 s	10	30
PLL	0.001 s	180	3200	1

.

In the next subsections, the study cases carried out to demonstrate the usefulness and effectiveness of the proposed method are shown. They are all developed using MATLAB software. It is worth noting that the setpoints adjustment of the ESs and the variation in the power injection of the RES are the scenarios that mostly occur in the operation of a µG with DG and distributed ESs; therefore, their study is fundamental (study cases 1 and 2). For validation purposes, the µG time-domain model is implemented, which includes the dynamics of the power and control stages of the ESs. These dynamic models represent an adequate point of comparison since they have been exhaustively studied in different publications reported in the literature [19,20,26,29]. The results obtained when the dynamic model reaches its steady state are the criteria used to carry out the validation; therefore, aspects related to transitory response are not considered in the discussion of the results. Finally, the performance of the model and control algorithms for the implementation of a quasi-stationary study case is presented (study case 3). It provides a handy tool to explore operating scenarios where the transient behavior does not play an important role. As a result, more extended operating scenarios can be implemented with the active participation of more devices and with less computational effort compared to dynamic simulation models.

In order to quantitively evaluate the performance of the MBFSM against the BFSM, a study case with different operating conditions is implemented. The µG of Figure 8 and the parameters of Table 1 are used. A set of n = 100 random values given by x_j are generated for the powers injected by the RES (P_REN2 and P_REN5). The active power is in the range of 0 to 1000 W, which corresponds to the maximum and minimum values used in the next study cases of Section 4.1 and Section 4.2. These arrays are mathematically represented by Equation (32). In the same way, 100 random values are generated for the operating conditions of the ESs. (m₂, m₄ y m₅). These values are contained in the range [−1, 1], which allows the ESs to operate in inductive or capacitive compensation mode. The mathematical representation of these arrays is shown in Equation (33).

$$P_{RENi}^n=\left\{\left[x_1;x_2;\dots ;x_n\right]\right\}\hspace{1em}\hspace{1em}\begin{array}{l}{x}_j\in \mathsf{\Re },\forall j=1,2,\dots ,n\\ x_j\in \left[0,1000\right]\mathrm{W}\\ \forall i=2,5.\end{array}$$

(32)

$$m_i^n=\left\{\left[x_1;x_2;\dots ;x_n\right]\right\}\hspace{1em}\hspace{1em}\begin{array}{l}{x}_j\in \mathsf{\Re },\forall j=1,2,\dots ,n\\ x_j\in \left[-1,1\right]\\ \forall i=2,4,5\end{array}$$

(33)

The convergence criterion used in the BFSM is defined by Equations (34)–(36) [34]. The voltage errors (ΔV_i) for all nodes in the electrical network are calculated from the voltages in the current iteration (k) and the previous iteration (k − 1), as shown in Equation (34). A convergence error (err) is defined. If Equations (35) and (36) are satisfied for all nodes, the method stops, and the steady-state solution is obtained.

$$\Delta V_i=V{s}_i^k-V{s}_i^{k-1}\hspace{1em}\forall i=1,2,\dots ,n$$

(34)

$$\left|real\left(\Delta V_i\right)\right|<err\hspace{1em}\forall i=1,2,\dots ,n$$

(35)

$$\left|imag\left(\Delta V_i\right)\right|<err\hspace{1em}\forall i=1,2,\dots ,n$$

(36)

The convergence criterion implemented for the MBFSM differs from that used in the BFSM. This criterion has been extensively discussed in Section 3.2 and is represented by:

$$\left|V_G-V_{Gn}\right|<err$$

(37)

Table 3. Comparison of the computational effort of solution algorithms.

	BFSM Iterations			MBFSM Iterations			Reduction of Average Iterations (%)
err	Max	Min	Average	Max	Min	Average
1 × 10−1	5	3	4	3	1	3	25.0
1 × 10−2	7	3	5	4	2	3	40.0
1 × 10−3	9	3	6	5	2	4	33.3
1 × 10−4	10	4	8	6	3	5	37.5
1 × 10−5	12	5	9	7	3	6	33.3
1 × 10−6	14	5	10	7	4	6	40.0
1 × 10−12	21	9	17	12	7	11	35.2
						Average	34.9

contains the results obtained for different convergence errors. The smallest convergence error is 1 × 10⁻¹, and the maximum one is 1 × 10⁻¹². For each method, the maximum, minimum, and average iterations are extracted, considering the 100 samples from the test bench. It can be seen that for all the convergence errors analyzed, the MBFSM obtains the solution in a smaller number of iterations, achieving a reduction of 34.9% in all the cases, which demonstrates the improvement in the computational effort when implementing the MBFSM.

4.1. Study Case: Constant Active Power and Changes in Reference Voltages

In this study case, the active power injected by the RES is kept constant, and the reference voltages of the ESs are modified. In this way, the source of the disturbances is fully identified, and the effect that it has on the µG is analyzed. The study case is divided into three time-intervals to show the different operating conditions of the ESs. The duration time of the intervals was determined based on the transient response of the µG. The time to reach the steady state depends on the severity of the disturbance and the collective action of the devices that conform to the µG. Therefore, this time can vary from one interval to another. In order to unify the criteria, it was adjusted to a value of 75 s, which was obtained experimentally. The active power injected by the RESs remains constant at the values P_REN1 = P_REN2 = 500 W. These values represent approximately 70% of the energy consumed in the µG under nominal conditions. Therefore, it is possible to reduce the power losses and adjust the reference voltages of the ESs to values close to the nominal voltage of the µG.

Figure 9 shows the behavior of the electrical quantities for the nodes Vs₂, Vs₄, and Vs₅, which correspond to the location of the ESs. The results obtained with the proposed steady-state model are identified with the subscript model. On the other hand, the results of the dynamic simulation are presented using root mean squarer values (RMS) and instantaneous values. When it corresponds to an RMS value, the dynamic subscript is added, and when the electric variable corresponds to an instantaneous value, it is indicated with the subscript dynamic_ins.

In the first interval, the reference voltages of the ESs are set to 119 V. It can be seen in Figure 9a that all the ESs reach the established reference voltage (red line). To achieve the objective of controlling the voltages on the Vs₂, Vs₄, and Vs₅ nodes to a value of 119 V, the ESs modify their voltage, reactive power, and operation mode (see Figure 9b,c). It can be seen that all the electrical variables in the dynamic simulation start with zero initial conditions, and after a transient, it reaches the steady state. It is essential to note the differences in the transient duration times for each voltage even though they include a controller with the same gains. These differences are presented because all the elements are interconnected, and the action of one ES affects the others. However, the ES local controller allows the gradual adjustment until the steady state is reached.

Regarding the steady-state model, this effect is already considered in the modeling and is adjusted using the solution and control algorithms. Results calculated with the steady-state model are represented as a constant value throughout the interval (black line). When comparing both results, it is verified that the results obtained by the steady-state model and the dynamic model are the same. From Figure 9b,c in the first interval, it can be seen that (i) ES2 operates in inductive mode with a voltage V_ES2 = 40.29 V and absorbs a reactive power Q_ES2 = 90.23 VAR, (ii) ES4 also operates in inductive mode but with V_ES4 = 22.17 V and Q_ES4 = 51.83 VAR, and (iii) ES5 operates in capacitive mode, delivering a reactive power Q_ES5 = −106.5 VAR and a voltage V_ES5 = 49.12 V.

The second interval shows the performance of the control algorithm when the ESs operate under saturation conditions. This effect occurs when the limits of voltage in the DC bus for the ES are exceeded. Regarding dynamic simulation, this limit is established by the element that powered the DC bus. If the required control action is higher than its limit, the ES voltage is adjusted to the maximum value and remains constant. If this restriction is not included in the steady-state model, the calculated ESs voltage will be wrong, and the collaborative operation of the ESs will not be appropriately modeled. In this regard, the reference voltage of the ES2 changes to Vs_ref₂ = 121 V, and the voltages of the other two ESs remain constant at Vs_ref₄ = 119 V and Vs_ref₅ = 119 V. The new value of Vs_ref₂ affects the entire µG causing changes in the voltages and modes of operation of the ESs since they try to reach their established reference voltage. Only Vs₅ reaches 119 V, while the other voltages are controlled at the values, Vs₄ = 119.2 V and Vs₂ = 120.5 V. The modification in Vs_ref₂ causes the change in the operation of the ES from inductive mode to a capacitive mode in order to increase the voltage in its PCC. However, it reaches its maximum voltage V_ES2 = 70.71 V and delivers a reactive power Q_ES2 = −137.9 VAR. The maximum voltage of the ESs is restricted by the voltage source that powered the DC bus. In this study case, it has a value of 100 V, and it is indicated in Figure 9b. The ES4 also reaches its maximum voltage V_ES4 = 70.71 V, but it operates in inductive mode, consuming a reactive power Q_ES4 = 135.7 VAR. The ES5 is kept in capacitive mode with a voltage V_ES5 = 29.7 V and a reactive power Q_ES5 = −68.45 VAR.

As was seen in the second interval, the modification in the reference voltage of one ES directly affects the other ESs. It is possible that in some cases, the ES reaches its operating limits, which is not adequate since this condition will not allow the appropriate handling of the power intermittencies in the µG. Therefore, it is essential to show that an increment or decrement of voltage is a collective task of the ESs, which is directly related to the assignment of the reference voltages of the ESs. At this point, the collaborative action is fundamental since if the reference voltages are arbitrarily adjusted, the ESs will have a conflict to establish their set points. This fact sometimes puts the ESs in saturation mode, resulting in an inadequate energy management which will directly impact on the µG’s voltage regulation and, consequently, its reliability.

Finally, in the third interval, the reference voltages are adjusted to the values Vs_ref₂ = 122 V, Vs_ref₄ = 122.5 V, and Vs_ref₅ = 122 V. In this interval, none of the ESs reaches their maximum voltage, allowing the voltage control at the nodes Vs₂, Vs₄, and Vs₅ with their reference voltages, respectively. The ES2 remains in capacitive mode with voltage and reactive power V_ES2 = 63.51 V and Q_ES2 = −132.3 VAR, respectively. The ES4 operates in capacitive mode, and its voltage and reactive power are V_ES4 = 50.31 V and Q_ES4 = −112.4 VAR, respectively. The ES5 operates in inductive mode with a voltage and reactive power V_ES5 = 1.695V and Q_ES5 = 4.112 VAR, respectively.

The results presented in this study case were allowed to validate the ES model and the proposed algorithms subjected to modifications in their reference voltages. The steady-state comparison of both simulation models showed the same results in the three simulation intervals; however, the computational effort for the dynamic simulation model is higher because it is necessary to compute the transient response to get the steady-state. The proposed model and the algorithms for its solution and control described the collective operation of multiple ESs in µGs adequately. The tests carried out properly showed the action in reactive power compensation mode and the saturation effect on the ESs.

4.2. Study Case: Constant Reference Voltages and Changes in ACTIVE POwer

In this case study, the reference voltages of the ESs are kept constant throughout the simulation, while the powers of P_REN2 and P_REN5 are modified. The purpose of this study case is to validate the model when it is subject to power disturbances, which is one of the most critical requirements in µGs with DG. The study case is divided into three intervals of 75 s. The values that P_REN2 takes in each interval is 500 W, 500 W, and 1000 W, while P_REN5 takes the values 0 W, 500 W, and 1000W. These values represent 35%, 70%, and 140% of the power consumed by the loads in nominal operating conditions. Therefore, the study case will allow observing the power handling of the ESs when there is a lack of dispatch (intervals 1 and 2) or surpluses of power (interval 3). The obtained results allowed demonstrating the proper operation of the proposal to manage power intermittencies and support voltage regulation locally simultaneously. The reference voltage of the ESs is adjusted to a fixed value Vs_ref₂ = Vs_ref₄ = Vs_ref₅ = 120V, which is determined based on the nominal voltage of the µG power supply.

Figure 10a shows the voltages at the nodes Vs₂, Vs₄, and Vs₅. It can be seen that Vs₂ and Vs₄ are controlled at 120 V established by their reference voltages, regardless of the disturbances produced by the RES. On the other hand, Vs₅ in the first interval reaches a value of 119.7 V (maximum allowed voltage), and in the following two ranges, it reaches the reference voltage of 120 V. The voltages obtained in the dynamic simulation show a transient at the beginning of each time interval, and later the steady state is reached.

The power disturbances produced in each interval cause the adjustment of the voltages and the operation modes of the ESs since the voltages in each PCC have to be controlled (see Figure 10b,c). During the first interval, the ESs operate in capacitive mode since the active power provided by P_REN2 = 500 W is not enough to compensate for the voltage drops in the µG. Consequently, the ESs operate in capacitive mode to increase the voltage in their respective PCC. In this interval, the voltages of the ESs are V_ES2 = 23.04 V, V_ES4 = 4.541 V, and V_ES5 = 70.71 V, and the reactive powers are Q_ES2 = −54.27 VAR, Q_ES4 = −10.89 VAR, and Q_ES5 = −136.6 VAR. On the other hand, in the second time interval, the active power increases to the values P_REN2 = 500 W and P_REN5 = 500 W. The increment in the delivered active power reduces the voltage drops in the µG, causing a reduction in the amount of reactive power delivered by the ESs. The ES2 and ES5 keep operating in capacitive mode, but now the voltages are modified to the values V_ES2 = 2.236 V and V_ES5 = 50.36 V, while the reactive powers are adjusted to Q_ES2 = −5.58 VAR and Q_ES5 = −109.7 VAR. Regarding the ES4, it changed its operation mode to inductive mode with a voltage V_ES4 = 22.53 V and a reactive power Q_ES4 = 53.11 VAR.

In the last interval, the total active power increases to 2000 W, distributed equally between P_REN2 and P_REN5, so this causes an increment in the node voltages and, consequently, the reactive power exchanged between the ESs and µG is modified. The ES2 and ES4 change their operating mode to inductive mode, and the ES5 decreases the amount of reactive power delivered. The new voltages at which the ESs operate are V_ES2 = 26.24 V, V_ES4 = 13.34 V, and V_ES5 = 7.318 V, and the reactive powers are Q_ES2 = 61.44 V, Q_ES4 = 15.51 VAR, and Q_ES5 = −17.53 VAR.

The results obtained with the steady-state model are equal to the ones obtained with the time-domain model, demonstrating that the model is valid to adequately and accurately represent the behavior of the multiple ESs under power disturbances. As identical results are obtained, the proposed model can be used and implemented to study large-scale systems with more extended operating scenarios and fewer computational efforts compared with the dynamic simulation models. An important aspect to highlight is the modification of the operating mode in the ESs when the active power increases. This fact is closely linked with his excellent ability to handle power locally. When the active power generated is lesser than the one consumed, the ESs operate in capacitive mode to support the voltage. On the contrary, when the generated power is higher than the one consumed, the ESs absorb reactive power and, consequently, the voltage is decreased.

4.3. Quasi-Stationary Simulation Model with Distributed Renewable Generation

This section implements a quasi-stationary study case. The study case allows for observing the correct performance of the proposal to implement more extended operating scenarios and with different operating conditions of the ESs and RESs distributed in the µG. A quasi-stationary study case with a variable power injection from the RESs with a duration time of 14 h and step increments of one minute is performed. P_REN2 and P_REN5 are modeled using a gaussian photovoltaic generation curve with a maximum power of 1000 W that occurs at 12 h in both cases. Two disturbances are implemented simulating the effect of shading on the photovoltaic panels. The disturbance affecting P_REN5 is a sudden power drop of 300 W that occurs in the interval from 10 to 11 h. The disturbance of P_REN2 is of magnitude 500 W, and it occurs in the interval from 12 to 13 h. In both disturbances, a time of 1 h was used to make more evident and clearly show the behavior of the ESs in the management of sudden power changes. Figure 11a shows P_REN2 and P_REN5. In Figure 11b, the behavior of the five voltages of the µG is shown. Figure 11c,d correspond to the voltages and reactive powers of the three ESs installed.

The power delivered by the RES before the 5 h and after the 18 h is zero, causing the voltages to remain constant with values Vs₁ = 119.95 V, Vs₂ = 119.88 V, Vs₃ = 119.6818, Vs₄ = 120 V, and Vs₅ = 119.72. This is obtained because the ES2 and ES5 reach their maximum voltage limits, i.e., V_ES2 = V_ES5 = 70.71 V operating in capacitive mode with a reactive power of Q_ES2 = Q_ES4 = −136.6 VARs. The Vs₄ voltage is the only one that is maintained at 120 V with a V_ES4 = 8543 V and a reactive power Q_ES4 = −20.45 VAR. This fact is due to that it did not reach its maximum allowable voltage for these operating conditions. As the power of P_REN2 and P_REN5 increase, the power losses in the distribution lines decrease, causing Vs₂ and Vs₅ to reach the 120 V reference at t = 7 h and t = 9 h, respectively, and, simultaneously, both ESs leave the operation in maximum reactive power compensation. It can be seen from 9 to 15 h that the voltages Vs₂, Vs₄, and Vs₅ remain at 120 V, which is possible because the ESs can mitigate the disturbances produced by the power injection. The first disturbance occurs in the interval from 10 h to 11 h and only affects P_REN5. The two ESs most affected by this disturbance are ES4 and ES5 due to the proximity with P_REN5. To maintain the voltage at 120 V, the ES4 increases the amount of reactive power consumed, while the ES5 simultaneously increases the amount of reactive power delivered. The second disturbance occurs in the time interval from 12 to 13 h, which is completely mitigated by the ES2. To achieve this, the ES2 decreased the amount of reactive power consumed; in a short time, it operated in capacitive mode. Regarding the voltages Vs₁ and Vs₃, they are more affected by the disturbances, which are directly reflected in a greater variation in their voltage.

Figure 12 illustrates the behavior for the voltages of the µG with the operation of the ESs (solid line) and without them (dotted line) under the same considerations of active power injection. The voltages of Vs₁ are shown in Figure 12a, the voltages of Vs₂ in Figure 12b, the voltages of Vs₃ in Figure 12c, the voltages of Vs₄ in Figure 12d, and finally the voltages of Vs₅ in Figure 12e. It is evident that when the ESs are not used, there are more considerable voltage variations, and their behavior is closely related to the power injection curve. On the other hand, by including the ESs, the voltage regulation is improved substantially. To quantify the voltage variations, the gradients Δ1 and Δ2 are calculated. For the calculation of the gradients, the value of 120 V is taken as reference, which corresponds to the nominal voltage of the power supply. The gradient Δ1 is obtained from the difference between the reference voltage and the maximum voltage, while Δ2 is calculated with the difference between the reference voltage and the minimum voltage. The gradients are graphically indicated in the voltages of Figure 12a–e, and their numerical values are concentrated in

Table 4. Voltage gradients of the µG.

	Without Control		With Control		Improvement on the Voltage Eduction (%)
Voltages	Δ1	Δ2	Δ1	Δ2	Δ1	Δ2
Vs1	−0.8604 V	1.3681 V	−0.6423 V	0.0427 V	NA	NA
Vs2	−0.6533 V	3.8775 V	0 V	0.1098 V	100	97.1683
Vs3	−0.1660 V	5.3455 V	0.2848 V	0.4661 V	NA	NA
Vs4	−0.6166 V	5.7038 V	0 V	0 V	100	100
Vs5	−0.5413 V	6.5444 V	0 V	0.2758 V	100	95.7857

NA: Not applicable because there is no ES at the node.

. The shaded area contains the voltages of the µG when the ESs are included. These voltages are zoomed for a better appreciation.

A negative gradient means that the maximum voltage exceeds the reference voltage, while a positive gradient indicates that the minimum voltage is below the reference. For a better context of the results, the gradients calculated without the ESs are first discussed. Column 1 in Table 4 shows the values of Δ1, which are negative and maintain a trend of increasing its value until reaching Δ1 = −0.1660 V, which corresponds to Vs₃. Then, the gradients of Vs₄ and Vs₅ decrease to the values Δ1 = −0.6166 V and Δ1 = −0.5413 V, respectively. This alteration in the trend is closely linked with the injection of active power in Vs₂ and Vs₅, which decreases the voltage drop in the µG and, consequently, increases the voltage in the farthest nodes. Regarding Δ2 values, they are all positive, showing a continuous magnitude increment as it moves away from the power supply. This behavior is originated because the minima voltages occur when the active power provided by the RES is zero, i.e., in t = 5 h and t = 19 h, generating a gradual voltage drop and being Vs₅ the most affected with a Δ2 =6.5444 V.

Columns 3 and 4 of Table 4 contain the gradients of the voltages when the ESs operate in the µG. The values of Δ1 and Δ2 are zero in Vs₄ since, throughout the study case, the ES regulates the voltage to 120 V established by its reference voltage. For its part in Vs₂ and Vs₅, Δ1 = 0 because it never exceeds 120 V. However, in the intervals where the ESs enter in saturation and, therefore, cannot control the voltage, a lower voltage was obtained, i.e., Δ2 = 0.1098 V and Δ2 = 0.2758 V, respectively. Because there are no ESs in nodes Vs₁ and Vs₃, a higher variation in voltage is obtained; in all the remaining cases, an improvement on the voltage reduction higher than 95% is obtained. The values obtained for Vs₁ are Δ1 = −0.6423 V and Δ2 = 0.0427 V. The voltage in Vs₃ never exceeds 120 V, causing both gradients to be positive with values Δ1 = 0.2848 V and Δ2 = 0.4661 V.

The results presented in this study case show the robustness and flexibility of both the proposed model and the solution and control algorithms to simulate a variable power injection. As can be seen, they are a great tool to analyze the operation of multiple ESs together with RES, including the analysis of µGs or large-scale distribution systems with DG. This characteristic also allows the development of collaborative control algorithms for managing a higher density of RES and alleviates the problems associated with intermittent energy sources without affecting the reliability of the µG.

4.4. Quasi-Stationary Simulation Model with Real Power Profiles Data

The injection of active power from more real conditions is an aspect of great importance to evaluate simulation models of electrical networks with DG. In this sense, the model was tested by using the daily and annual power profiles from the National Renewable Energy Laboratory (NREL) database [35]. The database contains the power profiles of approximately 6000 distributed photovoltaic solar plants in the USA. The power profiles have the information of the days included in a year. Each profile provides the active power injected in 5 min time intervals. The profiles allow evaluating the simulation models with DG under different operating conditions, i.e., sunny, partially cloudy, or cloudy days. Figure 13 depicts the results obtained from the simulation model and the real data. Figure 13a shows the active power profiles used by P_REN2 and P_REN5, respectively. These profiles correspond to the days 50 and 210 of a 50 MVA plant located in Arizona, USA. For the study case, the profiles were adjusted to a nominal power of 1 kVA to operate within the range of the power consumed by the loads in the µG. The power profile of P_REN2 describes the behavior of a sunny day in the interval time from 8 to 13 h, and with minimal power disturbances. Afterward, the power profile undergoes power variations of different magnitude and duration. Regarding P_REN5, its behavior corresponds to a partially cloudy day with power disturbances of greater severity and duration than those shown by P_REN2.

The results presented in Figure 13 retain the same format shown in Figure 11. It can be seen that the behavior of the electrical variables is similar to the ones obtained in the study case of Section 4.3, with differences associated with the new power profiles. Figure 13b shows the behavior of the µGs voltages, which remain constant during the null periods of power injection (before 8 h and after 18 h). In these time intervals, the ES2 and ES4 operate at their maximum voltage (V_ES2 = V_ES4 = 70.71 V) and in capacitive mode with reactive powers (Q_ES2 = Q_ES4 = −136.6 VAR), as shown in Figure 13c,d, respectively. On the other hand, Vs₄ remains constant throughout the simulation at a value of 120 V, which is closely linked to the fact that the ES4 never enters on saturation. Therefore, it can handle power locally and regulate the voltage in the PCC simultaneously. Active power injection comprises the time interval from 8 to 18 h. In this time interval, the ES2 and ES5 come out of saturation, allowing them to regulate the voltage in their PCC to a value of 120 V. The results presented demonstrate a proper behavior of the model and algorithms with real power profiles, where the ESs manage the power variations through three actions, i.e., the adjustment of their voltage magnitude, operation mode, and the amount of reactive power injected or consumed.

5. Comparative Analysis of Analytical Models for the ES Applied in Electrical Systems

In general, the analytical models for the ES used in electrical network simulations are divided into two groups, which are dynamic models and steady-state models.

Table 5. Comparative analysis of the mathematical models for the ES applied in electrical networks.

	Mathematical Model
	Dynamic	Steady-State
		Power Equations	Phasor
References	[15,26,36,37,38]	[21,24,25,39,40,41,42,43]	[22,23]	This work
Type of electric system	Distribution power systems Grid-connected µGs Standalone µGs	Distribution power systems Grid-connected µGs. Standalone µGs	• Grid-connected µGs	Grid-connected µGs Standalone µGs
ESs	Single Multiple	Single Multiple	• Single	Single Multiple
RES	Yes	Yes	No	Yes
Require equivalent power system circuits	No	Yes	Yes	No
Computational effort	High	Low	Low	Low
Additional constraints	No	Yes	No	No
Attributes	The control is implemented independently of the power stage [26] Controllers can be independent by ES and operate locally [38] It is suitable for power quality problems (i.e., transients, flicker, and harmonics) [36,37] It allows studying the transient stability in electrical systems. [15]	They are not capable of reproducing transitory phenomena [39] Different equations are formulated depending on the mode of operation and type of loads [21] Operating regions are determined to model the behavior of multiple Es adequately [41,42,43] The model needs additional restrictions such as topology, parameters, the physical location of the Es, operating limits [24,25,40]	They are not capable of reproducing transitory phenomena [39] The geometric relationships of the electrical quantities are used to obtain the model [22,23] Control laws are based on geometric relationships [22,23]
			• The calculation of geometric relationships has been limited to simplified electrical networks	Calculation of geometric relationships from a more realistic electrical network Active power injection associated with RES is considered A robust model that considers the effect of external disturbances

shows the essential characteristics, advantages, and disadvantages of both types of mathematical models. Dynamic models are widely used for studies of transient stability [15], power quality aspects related to short-term variations in electrical quantities [36], and harmonics [37]. These models allow the evaluation of overshoots and response times for the in-test system to assess their performance for different tasks, e.g., regulation purposes or disturbance mitigation, to avoid more severe network problems. Simulation models reported in the literature include the operation of a single ES up to several ESs, and in some cases, the models incorporate RES. These studies have been implemented in distribution power systems, grid-connected µGs, or standalone µGs. The power and control stages in these models are independent, which facilitates the implementation of simulation models with multiple ESs [38]. The power stage includes the physical restrictions of the ES, and the control stage is responsible for operating the ES under its operating principle [26]. However, dynamic simulation models require a high computational cost associated with both the robustness of the model and the small integration step in the solution process.

On the other hand, steady-state simulation models are obtained from phasor solution methods. These models are not suitable for reproducing the transitory phenomena of electrical networks. Still, they are characterized by having a low computational effort compared to dynamic models [39], which is a suitable feature to analyze large systems. The study of a single ES or multiple ESs in conjunction with RES has been carried out in distribution power systems, grid-connected µGs, and standalone µGs. These models are divided into two main groups.

• The first group is represented by power equations that describe the power exchange of the ES in the PCC; however, they have the constraint of being obtained from a simplified electrical network, indicating that their complexity in the study of multiple ESs will be higher. To adequately describe the behavior of the ESs, the models include additional restrictions related to the characteristics of the network, the physical location of the ES, and their operating limits [24,25,40]. Other alternatives carried out increase the complexity of the model since they include a more significant number of equations [21] or the design of external algorithms that restrict the operation of the ESs within established regions [41,42,43].
• The second group is called phasor models, and within this group is the proposed work. These models are characterized by providing a detailed representation of the geometric relationships of the electrical variables. Although some works have shown its use in grid-connected µGs by considering only one ES [22,23], the needs of modern power systems require the development of more robust models. In this regard, the proposed work outperforms previous works since it can deal with both grid-connected μGs and standalone μGs by considering single and multiple ESs. On the other hand, unlike reported works that consider simplified electrical networks, the proposed model is obtained from a more realistic electrical network by considering the injection of active power into the PCC and the effect of external disturbances on the behavior of the ES. These considerations increase its robustness, allowing for the obtaining of a generalized model suitable for the operation of multiple ESs without the need for additional restrictions. Consequently, a handy tool is obtained to study the requirements of modern energy networks.

6. Conclusions

In this work, a robust steady-state model for the ES is developed. All the considerations involved in its derivation allowed for the reaching of a generalized representation, making suitable its implementation in any physical location of the µG without additional restrictions. This feature expands its applicability in schemes that involve multiple ESs distributed in a µG. The ES model is solved by the MBFSM that allows for the obtaining of the steady-state solution of a µG with DG from the operating conditions of the ESs and the active power delivered by the RESs. Due to the ES model design and the MBFSM, it was possible to omit the forward propagation stage. This omitted step allows simplifying the solution algorithm and reducing the computational effort. Additionally, a control algorithm is implemented to control the µG voltage from the operating conditions of the ESs. The algorithm demonstrates an excellent performance for the control of multiple ESs operating within their limits or under saturation conditions (maximum available DC voltage). The obtained results allow the validation of both the proposed model and the proposed algorithms for the solution and control tasks since the steady-state results obtained by the proposed model are equal to the ones obtained with the dynamic simulation model when the steady-state is achieved. Results also validate of the proposal effectiveness to face active power variations, changes in the reference voltages of the ESs, and different values of loads at each node. The obtained results show an improvement higher than 95% for voltage control when the ESs are installed in the quasi-stationary simulation model; also, the analysis of quasi-stationary study cases with real active power profiles shows the versatility of the proposal for analyzing real scenarios. On the other hand, the MBFSM reduces the number of iterations by 34% on average compared to the BFSM.

The approach presented in this work forms the basis for future research on the operation of multiple ESs in conjunction with DG, allowing the implementation of collaborative control schemes based on intelligent algorithms that guarantee reliability in modern electrical networks. Additionally, it will be possible to analyze operating scenarios with uncertain information through the implementation of estimators and observers.