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

: 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.


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 CO2 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 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.

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 frequencydomain. Figure 1 shows the ES dynamic model, which includes its controller and power stage. The power stage of the ES incorporates a DC source, VDC, a DC-AC converter, and a low pass filter integrated by an inductor-capacitor (Lf and Cf). The converter operates as a single-phase full-bridge inverter consisting of four switches (S1-S4), 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, VES, 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, INC, 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, VES_mag. A phase-locked loop (PLL) is responsible for calculating the phase angle of the control voltage, VES_phase, which is delayed by 90 degrees with respect to the INC for reactive compensation mode. The capacitive or inductive mode is determined by the sign of VES_mag. If the sign is positive, INC leads VES_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].

ES Time-Domain Model
The configuration of the non-critical load, ZNC, 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 VES in both modes of operation has repercussions on the voltage applied to the non-critical load, VNC. 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].

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, G V  and its equivalent impedance, ZL. The ES is modeled as a controlled voltage source ES V  , which is connected in series between the non-critical load, ZNC, and the PCC Vs  . 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 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 Vs , and the non-critical load angle, ϕNC, is given by [31]: By solving Equation (1) for VNC, the following expression is obtained.
Once VNC is obtained, the next step is to calculate the reactive powers of the ES, QES, and the noncritical load, QNC. These are given in Equations (3) and (4) The expression for reactive power of the smart load (QSL) is given below Substituting (2)-(4) in (5), the Equation (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, VES, as well as the non-critical impedance, ZNC, 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.

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 timedomain 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, VG, that represents the utility grid and a set of elements connected through the line impedances (ZL1, ZL2,…., ZLn). The voltages of the µG are numbered sequentially, starting from the closest one to VG until reaching the last node (Vs1, Vs2, ..., Vsn). 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 (Z1, Z2,…, Zn) and non-critical loads (ZNC1, ZNC2,…, ZNCn). 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 (IREN1, IREN2,…, IRENn), where the subscript REN means renewable energies. The ESs are modeled as a controlled voltage sources (VES1, VES2,…, VESn) which modify their voltages based on their control variables (m1, m2, ..., mn).
The steady-state solution of the µG is obtained using the MBFSM algorithm. The algorithm receives as input the operating conditions of the ESs (m1, m2, ..., mn), and the active power injected by each of the RESs (PREN1, PREN2, ..., PRENn). In each iteration, the currents delivered by the current sources are calculated by Equation (7) to inject active power only.
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_ref1, Vs_ref2, ..., Vs_refn).
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 (Lf1-Cf-Lf2), 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.

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, By substituting Equations (9) and (10) in Equation (8), it is obtained that: Equations (14) and (15) Substituting Equations (13) and (16) in Equation (11 Now, defining 1 G V  as follows: equations can be expressed in a compact form by: Additionally, by applying a KVL in the j Vs  node, the equation associated with the SL is given by: 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 , 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 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 VG1 and 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  (19) and (20) is obtained. It can be seen in Equations (21) and (22) that the magnitude of VG1 and 1  depend on the impedances connected to the PCC (ZLj and Zj), the active power, PRENj, the magnitude of the current Ijk, the magnitude of the voltage on the PCC, Vsj, 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 Vsj, the BFSM provides the values of j Vs  and jk I  . The rest of the terms are parameters and operating conditions of the PCC, which are known.
The angle 0  depends on the power factor of ZLj as follows: The ES operating conditions are defined by the parameter mj, and this parameter takes values in the range of [−1, 1] and determines the magnitude of Vsi and the angle θ according to the following relationships.
The parameters a and b are used to simplify Equations (24) and (25). They are given by: 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.
The angle 4  is calculated by Equation (29) 2 2 It can be seen in Figure 5a,b that the geometric relationships take as reference the phasor NCj I  . 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.
In the same way, the phase angle for the voltage Vsi, θsi is calculated taking as reference θSj in the next equation.
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.

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 mj. 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. modified algorithm. 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 (PREN1, PREN2, .., PRENx) and the operation conditions of the ESs (m1, m2, …, my) 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, Vsn. 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, VG. 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 VG calculated in the backward propagation may not correspond to the real voltage of the power supply, VGn, the output error, ErrSal, is calculated as the absolute value of the difference between VG and VGn. 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 = 1x10 -12 was used, which was determined through an exhaustive experimentation. If the condition is not met, a new value of Vsn is calculated for the next iteration. The new value of Vsn 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 VG and VGn. The process is repeated until the maximum number of iterations, Niter, 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. 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, Vmax, then the algorithm stops. Vmax 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 Vmax are identified, and their mβ are not modified in this iteration. For the ESs that did not reach Vmax, 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. 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, VG, five distribution lines (ZL1, ZL2, ZL3, ZL4, and ZL5), three ESs (ES2, ES4, and ES5), two RES (IREN2 and IREN5) with powers (PREN2 and PREN5), three non-critical loads (ZNC2, ZNC4, and ZNC5), and five loads (Z1, Z2, Z3, Z4, and Z5), where Z2, Z4, and Z5 are the critical loads. The five node voltages are indicated as (Vs1, Vs2, Vs3, Vs4, and Vs5). Table 1 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].

Study Cases and Results
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 IREN, 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.
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 timedomain 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 xj are generated for the powers injected by the RES (PREN2 and PREN5). 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 Sections 4.1 and 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. (m2, m4 y m5). 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) The convergence criterion used in the BFSM is defined by Equations (34)-(36) [34]. The voltage errors (ΔVi) 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.
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: Table 3 contains the results obtained for different convergence errors. The smallest convergence error is 1e-1, and the maximum one is 1e-12. 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.

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 PREN1 = PREN2 = 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 Vs2, Vs4, and Vs5, 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 Vs2, Vs4, and Vs5 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 VES2 = 40.29 V and absorbs a reactive power QES2 = 90.23 VAR, (ii) ES4 also operates in inductive mode but with VES4 = 22.17 V and QES4 = 51.83 VAR, and (iii) ES5 operates in capacitive mode, delivering a reactive power QES5 = −106.5 VAR and a voltage VES5 = 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_ref2 = 121 V, and the voltages of the other two ESs remain constant at Vs_ref4 = 119 V and Vs_ref5 = 119 V. The new value of Vs_ref2 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 Vs5 reaches 119 V, while the other voltages are controlled at the values, Vs4 = 119.2 V and Vs2 = 120.5 V. The modification in Vs_ref2 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 VES2 = 70.71 V and delivers a reactive power QES2 = −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 VES4 = 70.71 V, but it operates in inductive mode, consuming a reactive power QES4 = 135.7 VAR. The ES5 is kept in capacitive mode with a voltage VES5 = 29.7 V and a reactive power QES5 = −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. 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.

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 PREN2 and PREN5 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 PREN2 takes in each interval is 500 W, 500 W, and 1000 W, while PREN5 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_ref2 = Vs_ref4 = Vs_ref5 = 120V, which is determined based on the nominal voltage of the µG power supply. Figure 10a shows the voltages at the nodes Vs2, Vs4, and Vs5. It can be seen that Vs2 and Vs4 are controlled at 120 V established by their reference voltages, regardless of the disturbances produced by the RES. On the other hand, Vs5 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 In the last interval, the total active power increases to 2000 W, distributed equally between PREN2 and PREN5, 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 results obtained with the steady-state model are equal to the ones obtained with the timedomain 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.

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. PREN2 and PREN5 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 PREN5 is a sudden power drop of 300 W that occurs in the interval from 10 to 11 h. The disturbance of PREN2 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 PREN2 and PREN5. 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. This fact is due to that it did not reach its maximum allowable voltage for these operating conditions. As the power of PREN2 and PREN5 increase, the power losses in the distribution lines decrease, causing Vs2 and Vs5 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 Vs2, Vs4, and Vs5 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 PREN5. The two ESs most affected by this disturbance are ES4 and ES5 due to the proximity with PREN5. 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 Vs1 and Vs3, 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 Vs1 are shown in Figure 12a, the voltages of Vs2 in Figure 12b, the voltages of Vs3 in Figure  12c, the voltages of Vs4 in Figure 12d, and finally the voltages of Vs5 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. 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 Vs3. Then, the gradients of Vs4 and Vs5 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 Vs2 and Vs5, 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 Vs5 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 Vs4 since, throughout the study case, the ES regulates the voltage to 120V established by its reference voltage. For its part in Vs2 and Vs5, Δ1 = 0 because it never exceeds 120V. 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 Vs1 and Vs3, 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 Vs1 are Δ1 = −0.6423 V and Δ2 = 0.0427 V. The voltage in Vs3 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.

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 PREN2 and PREN5, 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 PREN2 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 PREN5, its behavior corresponds to a partially cloudy day with power disturbances of greater severity and duration than those shown by PREN2.
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 (VES2 = VES4 = 70.71 V) and in capacitive mode with reactive powers (QES2 = QES4 = −136.6 VAR), as shown in Figure 13c,d, respectively. On the other hand, Vs4 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.

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 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 shortterm 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.

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. Funding: This research was funded by the "Consejo Nacional de Ciencia y Tecnología (CONACYT)" under the scholarship 728379.

Conflicts of Interest:
The authors declare no conflicts of interest.