New Decentralized Control of Mesh AC Microgrids: Study, Stability, and Robustness Analysis

: In this paper, we investigated the power sharing issues in mesh islanded microgrids that contain several distributed generators (DGs) and loads connected to different points of common coupling (PCC). Firstly, an improved decentralized droop control algorithm is proposed to achieve the active and reactive power sharing of different DGs in reconﬁgurable mesh islanded microgrids. Accurate power sharing was obtained even though line parameters or the mesh microgrid conﬁguration were unknown. Secondly a state-space model of the whole mesh microgrid was developed, considering several generators with their decentralized controllers, line feeders, and dynamic loads. This model was used to design parameters of droop controllers, to study the asymptotic stability and the robustness properties of the system. All strategies and analyses were validated by simulation based on the generic microgrid detailed in the standard IEEE 9bus test feeder. summarize, the main points are, the a new strategy sharing Unlike in the solutions mostly used by authors for power sharing in this modified strategy is decentralized, which it does to the the classical the by and The modified for and in and efficient synchronization strategy affect the The second is the and of a model the and This used


Introduction
Microgrids are able to integrate different distributed generator (DG) systems converting different types of renewable energy and to supply different types of loads. This gives a certain level of independence, allowing them to be connected or disconnected from the main grid. One of the challenges is to synchronize and connect all the distributed generators to an islanded microgrid, while providing the "plug and play" functionality and respecting the active and reactive power sharing between the different distributed generation units (DGs) [1,2]. Currently, the most used methods for power sharing and synchronization in literature are based on the droop control strategy [3][4][5][6][7][8]. However, most of the microgrids considered in these research works have only a single point of common coupling (PCC), which are connected to all DGs through converters and the loads like in Figure 1. In order to level up the independence of microgrids from the main grid by increasing the penetration rate of distributed generators in microgrids, the choice of mesh multi-PCC microgrids seems to be a good solution. However, the intermittency of renewable energy may cause more instability problems in mesh microgrids due to their higher level of complexity compared to microgrids with single points of common coupling. and different states of charge. The validated state-space model made the manipulation of the system so much easier and allowed many more studies to be carried out on the meshed microgrid. For example, one of many scenarios where the state-space model can be useful is when the CPL-injected power in the microgrid increases, and the state model can easily predict the right value of the capacitor that should be added to the PCC where the CPL is also connected. This means that many future studies can be conducted based on this developed and validated state-space model.

Power Sharing Using Droop Control Strategies
The traditional droop strategy imposing the electrical pulsation and the RMS voltage of a distributed generator via Equations (1) and (2) is efficient in microgrids with a single PCC shown in Figure 1, only if the effect of power line impedance is ignored [4]. In mesh multi-PCC reconfigurable microgrids, with many distributed generation sources (DGs) and loads randomly connected to different PCCs, classical synchronization methods and power sharing strategies based on droop control and used in mono-PCC microgrids are less efficient. Indeed, most droop methods assume that the transmission lines are purely inductive or resistive in nature which leads to a linear droop characteristic where active power sharing depends on either frequency or voltage magnitude. However, in microgrids with reconfigurable structure, the sources as well as the networks may become redundant. The changes in the network impedance whenever a branch or source is disconnected in a microgrid means that the active and reactive power cannot be totally decoupled. These changes affect the power sharing among the sources and reduce the stability of the microgrid for any disturbance [9,10].
Much research has focused on solutions that concern the active and reactive power sharing in droop-controlled mesh microgrids [11,12] either by using virtual impedance correction loop and a convergence acceleration strategy to compensate the offset created by the power line impedances and improve the reactive power sharing such as in ref [11], or by adopting a voltage compensation strategy to keep the bus voltage stable at the rated value and improve the reactive power sharing such as in ref [12]. However, only a few studies have focused on active and reactive power sharing in reconfigurable mesh microgrids such as in ref [10], where the authors propose the use of L1 adaptive methods for stable operation of a microgrid with wide range of R/X ratios.
Another problem is the possible unstable behavior of the microgrid caused by the interaction between the DGs and loads as well as the changes in the network configuration. The system stability depends on the type of loads being connected to the microgrid, especially the ones supplied through tightly regulated power converters. These loads behave as constant power loads (CPL) and may cause system instability [13][14][15]. In the literature, the main studies are based on linearization techniques [16]. The non-linear models of the considered power systems are linearized around an operating point and then studied using linear analysis tools. Moreover, linearization tools only predict the stability of the system for small perturbations [17,18] and cannot guarantee the system stability for large perturbations.
In this paper, a mesh microgrid constituting several PCCs connected to two DGs (DG1 with a nominal power of three megawatts and DG2 with a nominal power of two megawatts) through two transformers (6 kV/20 kV) as well as three loads (LOAD1 with a nominal power of one megawatt, LOAD2 with a nominal power of 1.2 megawatts, and LOAD3 with a nominal power of 1.5 megawatts) are considered. The different PCCs are interconnected with power supply lines, modeled by resistance, inductance and capacitance (RLC) circuits inspired by an IEEE 9-bus test feeder ( Figure 2). All the microgrid data are presented in Tables 1 and 2.
With: m = ∆ , n = ∆ and J = K ∫ In the absence of information on E the coefficient J is set to 0 and we return to the classical droop.   In order to ensure an accurate power sharing and to provide the "plug and play" function adapted to the considered multi-PCC mesh microgrid, a new droop control and synchronization strategies is applied, such as the one demonstrated in ref [19]. This new droop control in islanded mode consists of removing the decoupling between the Sustainability 2021, 13, 2243 4 of 25 active and reactive power caused by the line impedances and the change in the microgrid configuration. The first DG connected to the microgrid imposes its voltage and frequency. For other DGs, an adapted synchronization strategy is proposed and applied before their interconnection to the microgrid. The simulation results using a model developed in the Simscape feature of MATLAB/Simulink environment enabled the validation of the efficiency of this new droop control.
To study the stability of the microgrid, a state model of the entire meshed microgrid was developed considering its power lines and loads including a CPL as well as its generators, with their decentralized controllers integrating modified droop algorithms. The comparison of the results obtained from the stablished state model and the developed model in Simscape environment confirm the validity of the proposed microgrid state model. This validated state model was then used to study the microgrid stability by calculating its Jacobean matrix and its eigenvalues at each operating point. Based on these results, the root locus, and eigenvalues sensitivity analysis, depending on the parameters of the established model, were performed. The latter made it possible to find the origin of different frequency components by observing the dynamic behavior of different state variables. The analysis of the system stability was then performed by observing the evolution of the system eigenvalues with respect to droop parameters and different load power variations. It will be shown that the system stability may be affected due to unadapted droop control parameters and/or CPL levels.
The established state-model was also used to verify the robustness of the considered mesh microgrid, including the new droop-control that ensured active and reactive power sharing to the modifications caused by the connection of DGs and large load variations (changes in network configuration).
To summarize, the main points of this paper are, firstly, the development of a new droop control strategy for power sharing especially reactive power in mesh islanded microgrids. Unlike in the solutions mostly used by authors for power sharing in reconfigurable mesh microgrids, this modified droop control strategy is decentralized, which means that it does not require any sort of communication between the DGs connected to the microgrids. Thanks to the non-linear term added to the classical droop control, the power coupling phenomenon created by the line impedances and the change in the microgrid configuration is removed. The modified droop control strategy was proven efficient for active and reactive power sharing in different microgrid structures. Additionally, a fast and efficient synchronization strategy did not affect the power sharing. The second major point is the development and validation of a state-space model describing the mesh microgrid and its control. This state-space model was used to study the microgrid stability especially in the presence of active loads such as CPLs. The state model was also used to test the robustness of the modified droop control under different microgrid configurations and different states of charge. The validated state-space model made the manipulation of the system so much easier and allowed many more studies to be carried out on the meshed microgrid. For example, one of many scenarios where the state-space model can be useful is when the CPL-injected power in the microgrid increases, and the state model can easily predict the right value of the capacitor that should be added to the PCC where the CPL is also connected. This means that many future studies can be conducted based on this developed and validated state-space model.

Power Sharing Using Droop Control Strategies
The traditional droop strategy imposing the electrical pulsation and the RMS voltage of a distributed generator via Equations (1) and (2) is efficient in microgrids with a single PCC shown in Figure 1, only if the effect of power line impedance is ignored [4].
where P i and Q i are the measured values of the active and reactive power of the ith DG, P in and Q in are their nominal values, ω n and E n are the nominal values of the pulsation and voltage of the ith DG, ∆ ω and ∆ E are the permissible variations of pulsation and voltage, and m i and n i are the droop control coefficients. The active and reactive power sharing method is based on a droop control algorithm that sets the frequency and voltage amplitude at the associated PCC for each DG according to Equations (1) and (2). This droop control does not ensure efficient reactive power sharing even in single PCC micro-grids due to the line impedances. A droop control strategy that has been proven effective for active and reactive power sharing in mono-PCC microgrids (3) and (4) was applied to the complex mesh microgrid in Figure 2, and considering one of the PCC voltages as reference potential of the studied microgrid. The simulation results are presented in Figures 3 and 4.
Sustainability 2021, 13, x FOR PEER REVIEW 6 of 27 (a) (b) Figure 3. Evolution of the distributed generator (DG) active (a) and reactive (b) powers (first DG supplies the microgrid up to 5 s, the second DG is connected to the microgrid at 5 s).

Synchronization Strategy in Mesh Islanded Microgrid
Due to the complexity of the mesh microgrids and the intermittency of renewable energies, DGs frequently connect to and disconnect from the microgrid, so a fast and efficient synchronization method is required. To achieve the synchronization of the ith DG to the ith PCC before their interconnection, the voltage amplitude , the pulsation ω and the phase θ of the ith PCC must be approximately equal to those of the ith DG ( , ω , θ ) [2]. To achieve this objective, the errors between the amplitudes, pulsations, and phases of both sides (the ith DG and the ith PCC) are forced to zero by adding pure integral controllers to the ith DG droop control, as shown in Equations (9) and (10) [19]. It should be noted that the binary coefficient in these equations is equal to one only during the synchronization interval, and equal to zero otherwise. For the pulsation and phase synchronization, the pure integrators are added to the pulsation droop control, Equation (9), and for the voltage synchronization, the pure integrator is added to the voltage droop control, Equation (10).  The proposed angular droop aims to indirectly control the voltage at the PCC and its phase to be equal to the nominal values (i.e., E in and δ in ). The added integrators can minimize the static error between the feedback signal and the corresponding nominal values. By choosing identical K a and K e for each generator, an accurate active and reactive power sharing is achieved which no longer depends on the system impedance and is unaffected by the digital errors and disturbances [9].
In practice, real microgrids can have several PCCs, interconnected by several supply lines with non-negligible impedances. In order to study such a microgrid, a slightly modified IEEE 9bus test feeder is considered, composed of two DGs and three loads, interconnected by RLC lines (Figure 2). It can be noticed that the application of the droop method, based on Equations (3) and (4), on the multi-PCC microgrid of Figure 2 leads to a perfect active power sharing when the second DG is connected at 5 s ( Figure 3a). However, the reactive power sharing when the second DG is connected is not achieved (Figure 3b). It can also be remarked that unacceptable disturbing active and reactive power peaks occur due to the lack of synchronization (Figure 3a,b at 5 s).
In mesh microgrids, each line, connecting the ith PCC to the jth, has a non-negligible inductance (λ i,j ) and resistance (ρ i,j ), leading to a line voltage drop between these two PCCs that creates a coupling between the active power (P i,j ) and reactive power (Q i,j ) exchanged between ith and jth PCCs, according to Equation (5): where ϕ is the phase shift between the phase voltage E j and the phase current I i,j . To achieve an efficient reactive power sharing in this type of mesh microgrids, the voltage equation in (2) is modified by adding a decoupling term (see Equation (7)) removing the coupling phenomenon between active and reactive power [12]. This coefficient, called J i , is estimated using a PI controller, forcing the suppression of the error ε i defined in Equation (8). In steady state, when the error ε i tends to zero, the reactive power sharing is well ensured between the DGs. This approach allows the primary control of the voltage of the ith DG, regardless of the operating point of the loads. It should be noted that E re f in Equation (8) is the measured value of one of the PCC voltages of the studied mesh microgrid acting as the pilot node and considered as its reference potential. The rated voltage of this pilot node is called E n .
With : In the absence of information on E re f the coefficient J is set to 0 and we return to the classical droop.

Synchronization Strategy in Mesh Islanded Microgrid
Due to the complexity of the mesh microgrids and the intermittency of renewable energies, DGs frequently connect to and disconnect from the microgrid, so a fast and efficient synchronization method is required. To achieve the synchronization of the ith DG to the ith PCC before their interconnection, the voltage amplitude E pcci , the pulsation ω pcci and the phase θ pcci of the ith PCC must be approximately equal to those of the ith DG (E DGi , ω DGi , θ DGi ) [2]. To achieve this objective, the errors between the amplitudes, pulsations, and phases of both sides (the ith DG and the ith PCC) are forced to zero by adding pure integral controllers to the ith DG droop control, as shown in Equations (9) and (10) [19]. It should be noted that the binary coefficient B i in these equations is equal to one only during the synchronization interval, and equal to zero otherwise. For the pulsation and phase synchronization, the pure integrators are added to the pulsation droop control, Equation (9), and for the voltage synchronization, the pure integrator is added to the voltage droop control, Equation (10).
To prove the efficiency of the proposed synchronization and power sharing strategies, the mesh microgrid in Figure 2 is modeled using the Simscape toolbox of Matlab/Simulink. Source 1 and Source 2 of Figure 2 are modeled by two controllable voltage sources shown in Figure 4, connected to two different PCCs and controlled by E abc re f i , which was generated using the modified droop strategy described by Equations (9) and (10), a power calculation bloc, a power filter, park transform and inverse park transform, as explained in the equivalent synoptic diagram describing a droop-controlled DGi in Figure 5. The main microgrid parameters are listed in Table 1. The powers of sources and loads are listed in Table 2.     Based on this Simscape model, the efficiency of the proposed synchronization and power sharing strategies assumed in Equations (9) and (10) is evaluated. Figure 6a and 6b show the evolution of the active and reactive powers delivered by DG1 and DG2. In the beginning, the first DG (DG1) imposes the frequency of the microgrid as well as the voltages of each PCC; up to t = 5 s, DG1 supplies the loads connected to the microgrid. The second generator (DG2) is synchronized during the interval 1-5, and then it is connected to the microgrid at t = 5 s. The power-sharing of the active and reactive powers are ensured in steady state without being affected by the synchronization procedure. In addition, compared to the previous results shown in Figure 3a,b, the power peaks appearing after the connection of DG2 to the microgrid are considerably attenuated. It should be noted that these performances are maintained with a higher number of DGs connected to the microgrid, even though the results are not presented in this paper.
second generator (DG2) is synchronized during the interval 1-5, and then it is connected to the microgrid at t = 5 s. The power-sharing of the active and reactive powers are ensured in steady state without being affected by the synchronization procedure. In addition, compared to the previous results shown in Figure 3a,b, the power peaks appearing after the connection of DG2 to the microgrid are considerably attenuated. It should be noted that these performances are maintained with a higher number of DGs connected to the microgrid, even though the results are not presented in this paper.
(a) (b) Figure 6. Evolution of the DG active (a) and reactive (b) powers (first DG supplies the microgrid up to 5 s, the second DG is connected to the microgrid at 5 s).
In order to highlight the effect of the loss of the information on , a simulation test was made using the developed model in Simscape environment, and the results are shown in Figure 5. It consists of connecting both DGs to the microgrid at 0 s using the control strategy defined in Equations (6) and (7). The information on is lost at 4 s and becomes available again at 8 s, and in this scenario Figure 7a,b show the evolution of the active and reactive powers, respectively. The active power sharing is not affected by the loss of information on because it only depends on the frequency, being the same in all the microgrid PCCs. However, the reactive power sharing is totally lost when the information on is absent, but it is easily regained when the information on becomes available again. In order to highlight the effect of the loss of the information on E re f , a simulation test was made using the developed model in Simscape environment, and the results are shown in Figure 5. It consists of connecting both DGs to the microgrid at 0 s using the control strategy defined in Equations (6) and (7). The information on E re f is lost at 4 s and becomes available again at 8 s, and in this scenario Figure 7a,b show the evolution of the active and reactive powers, respectively. The active power sharing is not affected by the loss of information on E ref because it only depends on the frequency, being the same in all the microgrid PCCs. However, the reactive power sharing is totally lost when the information on E ref is absent, but it is easily regained when the information on E ref becomes available again. Even though the considered microgrid including its loads and DGs are well controlled, the synthesis of the droop control parameters is not based on a methodical approach; they are tuned by trial-and-error method using Simscape.
This method cannot guarantee the system stability when the parameters of the load change or additional dynamical constant power loads (CPL) are connected to the microgrid. For example, when a CPL is added to the microgrid it is not evident to foresee the Even though the considered microgrid including its loads and DGs are well controlled, the synthesis of the droop control parameters is not based on a methodical approach; they are tuned by trial-and-error method using Simscape.
This method cannot guarantee the system stability when the parameters of the load change or additional dynamical constant power loads (CPL) are connected to the microgrid. For example, when a CPL is added to the microgrid it is not evident to foresee the system stability with the variation of its absorbed active power. Figure 8 shows that, as previously, the first DG1 imposes the frequency of the microgrid as well as each PCC voltage of up to t = 5 s and supplies the R-L loads connected to the microgrid. The second generator (DG2) is synchronized and connected to the microgrid at t = 5 s. A dynamic CPL is connected parallel to load 1 (see Figure 2) at t = 7 s, absorbing an active power of 120 kW and zero reactive power. The active power absorbed by this CPL is increased to 500 kW at t = 8 s. It can be remarked that the connection of the CPL absorbing 120 kW does not impact the stability, but when its absorbed power increases to 500 kW the system becomes unstable. Hence, for each modification, the trial-and-error method should be applied again in order to find the eventual new values of the microgrid control parameters which may ensure its stability. In order to overcome this constraint, one solution consists of establishing the state-model that represents the microgrid including its control and studying the microgrid stability by calculating its Jacobean matrix and its eigenvalues at each operating point, as is detailed in the next section.

System Modeling for its Stability Analysis
Due to the interaction between the DGs and the loads, the stability of the microgrids is strongly influenced. Therefore, to study the stability and robustness of the droop-controlled microgrid, a complete mathematical dynamic model of the whole system is established in this section.

System Modeling for Its Stability Analysis
Due to the interaction between the DGs and the loads, the stability of the microgrids is strongly influenced. Therefore, to study the stability and robustness of the droop-controlled microgrid, a complete mathematical dynamic model of the whole system is established in this section.

State-Space Model of a Distributed Generator (DG)
The equations describing the behavior of the ith DG are presented on the equivalent synoptic diagram of Figure 9. It is composed of four cascaded blocs for which the relationship between their outputs and inputs are presented in the following.

System Modeling for its Stability Analysis
Due to the interaction between the DGs and the loads, the stability of the microgrids is strongly influenced. Therefore, to study the stability and robustness of the droop-controlled microgrid, a complete mathematical dynamic model of the whole system is established in this section.

State-Space Model of a Distributed Generator (DG)
The equations describing the behavior of the ith DG are presented on the equivalent synoptic diagram of Figure 9. It is composed of four cascaded blocs for which the relationship between their outputs and inputs are presented in the following. In the first bloc (Figure 9), the ith DG measured voltages ( , ) and currents ( , ) are used to calculate the instantaneous active and reactive powers using classical equations (11) and (12). In the first bloc (Figure 9), the ith DG measured voltages (E i d , E i q ) and currents (i di , i qi ) are used to calculate the instantaneous active and reactive powers using classical Equations (11) and (12).
In the second bloc, the instantaneous active and reactive powers (P i and Q i ) are filtered using a first-order low-pass filter to obtain the average (filtered) values P f i and Q f i using Equations (13) and (14). It should be emphasized that the cut-off frequency ω f of the filter is related to the dynamics of the droop control-loop.
The third bloc applies the modified droop control in Equations (6)- (8), using P f i and Q f i to achieve the active and reactive power sharing.
The fourth bloc models the delay imposed by the voltage source inverter (VSI), controlling the ith DG output voltages by means of a second-order filter. This filter has a faster dynamic compared to the external droop control loop ω c ω f . The relationship between the output voltage of the ith DG and its voltage reference is expressed by the transfer function defined in Equation (15): Figure 10 shows a microgrid inspired from the IEEE 9bus test feeder that differs slightly from the microgrid of Figure 2 and contains an additional CPL. In fact, this type of load is increasingly present in microgrids and imposes more severe constraints to the system stability. The considered microgrid is composed of two DGs powering three classical inductive loads, modeled by serial R-L circuits, and a CPL load. They are interconnected by RLC power lines to represent a mesh multi-PCC microgrid. The line connecting the PCC i and PCC j have a resistance r ij , an inductance l ij , and a capacitance of c ij . Their values are determined knowing the distance between the considered PCCs. The capacitance c ij is connected to the PCC j . connecting its and , where and ∈ [1, 2, … , ] and ≠ ( Figure 10). Figure 10. IEEE 9bus test feeder. Figure 10. IEEE 9bus test feeder.

Microgrid Structure and Modeling
In order to study the considered microgrid having two DGs, the state equations are given using a common reference frame for expressing the different variables. The used d-q reference frame is defined in a way that its d-axis is oriented toward the first DG voltage vector (E 1 Figure 11). Then, the d-axis is shifted by δ 1 with respect to the α-axis, where δ 1 is the argument of the DG1 voltage vector (E 1 = E 1 e jδ 1 = E 1 α + jE 1 β ). The rotating d-q reference frame turns with an electrical speed of ω 1 = dδ 1 dt with respect to the fixed reference frame α − β. When the different DGs are synchronized and connected to the microgrid, their pulsations vary simultaneously. Then, the pulsations of the DG1 and DG2 voltages (ω 1 = dδ 1 dt and ω 2 = dδ 2 dt ) are the same (i.e., ω 1 = ω 2 = ω com )). Knowing that the phase voltages of DG1 and DG2 (connected to PCC1 and PCC2) have the RMS values of E 1 and E 2 , the α − β components of E 1 and E 2 are obtained using the Concordia transformation T t 32 in (16) and (17): where : Then, Equations (18) and (19) give the relationship between the d-q components of the voltage vectors at PCC 1 and PCC 2 and their α − β components. DG1 : DG2 : Equation (20) presents a network model of N power lines of a mesh microgrid, interconnecting its PCC i and PCC j , where i and j ∈ [1, 2, . . . , N] and i = j (Figure 10).
E i d and E i q are the d-q components of PCC i voltage. The loads considered in this work were serial passive R-L loads and an active CPL ( Figure 10). With the sum of the current connected to PCC i being zero (Kirchhoff current law), Equation (21) relates to all PCCs except the ones connected to the DGs for which the voltages are imposed by their droop controllers.
It should be noted that ∑ j i ji in Equation (17) is the sum of the currents of the power lines connected to PCC i . i CPLid and i CPLiq are the d-q components of the current absorbed by the active CPL i . The non-linear Equation (22) allows the determination of these components knowing the active and reactive powers P CPLi and Q CPLi absorbed by CPL i .
The d-q components of the R-L load current connected to PCC i (i Loadi d and i Loadi q ) can be determined using Equation (23): where R Loadi and L Loadi represent the resistance and inductance of the RL load. The overall model of the microgrid presented in Figure 10 including its control, defined by the previous equations, has 40 state variables which constitute the elements of the state vector x, described in the following (24): Considering the non-linear Equations (7), (11), (12), (19) and (21), the overall microgrid model is non-linear and can be finally put into the following form: .

Validation of the State Model
The results concerning the evolution of the microgrid active and reactive powers, obtained under the same conditions using the Simscape model and the established state model, were used to validate the established state model. The same simulation scenario was considered for obtaining the simulation results using both models. Both DGs were connected to the microgrid at 0 s. The simulation results of both considered models are presented in Figure 12a,b, illustrating the evolution of the active and reactive powers of DG1 and DG2. Both models led to similar results; therefore, the established state model is considered validated. Thus, it can be used for studying the microgrid stability, integrating its control. To further confirm the validity of the developed model, Figure 12c shows a zoomed-in view of the transitory state after the connection of both DGs at 0 s.

Validation of the State Model:
The results concerning the evolution of the microgrid active and reactive powers, obtained under the same conditions using the Simscape model and the established state model, were used to validate the established state model. The same simulation scenario was considered for obtaining the simulation results using both models. Both DGs were connected to the microgrid at 0 s. The simulation results of both considered models are presented in Figure 12a,b, illustrating the evolution of the active and reactive powers of DG1 and DG2. Both models led to similar results; therefore, the established state model is considered validated. Thus, it can be used for studying the microgrid stability, integrating its control. To further confirm the validity of the developed model, Figure 12c

Jacobian Matrix Eigenvalues
In order to study the stability of the microgrid, defined by the continuous nonlinear differential system (24), and to analyze its dynamic behavior, which is affected by several control parameters, an eigenvalue sensibility analysis was performed by means of the calculation of the eigenvalues of the Jacobian matrix at the operating point (x ) = 0.
The system, i.e., the microgrid, will be asymptotically stable around the operating point x if the real parts of all eigenvalues are strictly negative.
The system parameters are given in Tables 3-5. The classical loads were modelled by serial R-L loads. The values of R and L were chosen depending on the operating point of these loads. The CPL was defined by constant values of its absorbed active and reactive powers, which may vary depending on its chosen operating point. The parameters of the DGs corresponded to their rated powers (available powers) and the parameters imposing the dynamic behavior of their control ( , , , permissible variations of pulsation ∆ and voltage ∆ ).
The eigenvalues of the system at the considered operating point, corresponding to the parameters of the loads and DGs given in Tables 4 and 5, are shown in Table 6. It can be observed that four of the Jacobian matrix eigenvalues are the negative real numbers, and that the other 36 eigenvalues are two by two the complex conjugate numbers, having negative real parts.
All of the Jacobian matrix eigenvalues have negative real parts, therefore the microgrid including its DG's control and loads can be considered as a stable system around the considered operating point. The study of the considered mesh microgrid stability around other operating points can also be conducted based on the established state model. Table 3. Parameters of the considered microgrid power lines.

Jacobian Matrix Eigenvalues
In order to study the stability of the microgrid, defined by the continuous nonlinear differential system (24), and to analyze its dynamic behavior, which is affected by several control parameters, an eigenvalue sensibility analysis was performed by means of the calculation of the eigenvalues of the Jacobian matrix The system, i.e., the microgrid, will be asymptotically stable around the operating point x 0 if the real parts of all eigenvalues are strictly negative.
The system parameters are given in Tables 3-5. The classical loads were modelled by serial R-L loads. The values of R and L were chosen depending on the operating point of these loads. The CPL was defined by constant values of its absorbed active and reactive powers, which may vary depending on its chosen operating point. The parameters of the DGs corresponded to their rated powers (available powers) and the parameters imposing the dynamic behavior of their control (ω f , ω c , ξ, permissible variations of pulsation ∆ ω and voltage ∆ E ). Table 3. Parameters of the considered microgrid power lines.

Points of Connection
The eigenvalues of the system at the considered operating point, corresponding to the parameters of the loads and DGs given in Tables 4 and 5, are shown in Table 6. It can be observed that four of the Jacobian matrix eigenvalues are the negative real numbers, and that the other 36 eigenvalues are two by two the complex conjugate numbers, having negative real parts. All of the Jacobian matrix eigenvalues have negative real parts, therefore the microgrid including its DG's control and loads can be considered as a stable system around the considered operating point. The study of the considered mesh microgrid stability around other operating points can also be conducted based on the established state model.
In order to perform the sensitivity analysis and find the origin of different frequency components, depending on the parameters of the established model, the parameters of the different constituents of the microgrid, including its control parameters were modified separately. In the following sections, the evolution of eigenvalue trajectories with variation of the parameters of each microgrid constituent is discussed to verify their impact on the system stability.

Eigenvalue Trajectory and Sensibility Analysis for Different Values of Power Line Parameters
To vary the parameters of the 6 power lines of the microgrid, the values of the lengths of all of these power lines are modified by 10 steps of 10% (for each line, from its initial length to two times this length), the Jacobean matrix eigenvalues are determined, and their evolution shown on Figure 13. The variation of the parameters of the power lines only impacts the eigenvalues λ 1 to λ 12 , which are two by two complex conjugate numbers; the other eigenvalues λ 13 to λ 40 do not vary significantly with the evolution of the power line parameters. The real parts of the impacted eigenvalues (λ 1 to λ 12 ) become closer to zero but remain sufficiently negative to not impact system stability when the line parameters are increased within the considered interval. However, if the lengths of the power lines become excessively long, the real parts of the impacted eigenvalues may become positives and cause system instability. To vary the parameters of the 6 power lines of the microgrid, the values of the lengths of all of these power lines are modified by 10 steps of 10% (for each line, from its initial length to two times this length), the Jacobean matrix eigenvalues are determined, and their evolution shown on Figure 13. The variation of the parameters of the power lines only impacts the eigenvalues to , which are two by two complex conjugate numbers; the other eigenvalues to do not vary significantly with the evolution of the power line parameters. The real parts of the impacted eigenvalues ( to ) become closer to zero but remain sufficiently negative to not impact system stability when the line parameters are increased within the considered interval. However, if the lengths of the power lines become excessively long, the real parts of the impacted eigenvalues may become positives and cause system instability.

Eigenvalue Trajectory and Sensibility Analysis for Different Values of DG Parameters
The same approach as in Section 4.2. is used to study the impact of other parameters that may contribute to the instability of the system. For the considered DGs, the parameters which may impact the system stability are the cut-off frequency for power calculation (Equation (13) and (14)), the value of the filter's cut-off frequency of the second order filter (Equation (15)), as well as the frequency drop ∆ and the voltage drop ∆ in Equations (6)- (8).
First, the cut-off frequency is modified in 10 steps (from its initial value to 10 times this value), the Jacobean matrix eigenvalues are determined, and their evolutions are illustrated in Figure 14. The variation of impacts the real negative eigenvalues , , and as well as the eigenvalues , which are two complex conjugate numbers. The other eigenvalues do not vary significantly with the evolution of . It can be seen that when increases, the eigenvalues 33 , 34 , 36 , 37 and 38 are shifted to the left while 35 moves very slightly toward zero, but its value does not change considerably and stays largely negative (from −4 to −3.6). The system always remains stable when increases, but the system dynamics become increasingly faster. To examine the impact of in Equation (13), representing the delay imposed by the voltage source inverter (VSI) controlling the DG, on the system stability, its value is largely decreased (from its initial value 1000 rad/s to 50 rad/s). The eigenvalues most sensitive to the variation of are to as well as and , which are two by two of the complex conjugate numbers. As shown in Figure 15, the real parts of the impacted eigenvalues move to the right when decreases and tend to be positive when becomes inferior to 50 rad/s, which causes the system instability.

Eigenvalue Trajectory and Sensibility Analysis for Different Values of DG Parameters
The same approach as in Section 4.2. is used to study the impact of other parameters that may contribute to the instability of the system. For the considered DGs, the parameters which may impact the system stability are the cut-off frequency ω f for power calculation (Equations (13) and (14)), the value of the filter's cut-off frequency ω c of the second order filter (Equation (15)), as well as the frequency drop ∆ ω and the voltage drop ∆ E in Equations (6)- (8).
First, the cut-off frequency ω f is modified in 10 steps (from its initial value to 10 times this value), the Jacobean matrix eigenvalues are determined, and their evolutions are illustrated in Figure 14. The variation of ω f impacts the real negative eigenvalues λ 35 , λ 36 , λ 37 and λ 38 as well as the eigenvalues λ 33 , λ 34 which are two complex conjugate numbers. The other eigenvalues do not vary significantly with the evolution of ω f . It can be seen that when ω f increases, the eigenvalues λ 33 , λ 34 , λ 36 , λ 37 and λ 38 are shifted to the left while λ 35 moves very slightly toward zero, but its value does not change considerably and stays largely negative (from −4 to −3.6). The system always remains stable when ω f increases, but the system dynamics become increasingly faster. To validate these results, the Simscape simulations were made to verify that any value of inferior to 50 rad/s can cause system instability. The results of the first simulation, made using the Simscape model with = 50 rad/s, are illustrated on Figure 16a, and those of the second simulation with = 45 rad/s are shown on Figure 16b. These results confirm the system stability for > 50 rad/s, as predicted by the stability analysis based on the established state model.  To examine the impact of ω c in Equation (13), representing the delay imposed by the voltage source inverter (VSI) controlling the DG, on the system stability, its value is largely decreased (from its initial value 1000 rad/s to 50 rad/s). The eigenvalues most sensitive to the variation of ω c are λ 27 to λ 32 as well as λ 39 and λ 40 , which are two by two of the complex conjugate numbers. As shown in Figure 15, the real parts of the impacted eigenvalues move to the right when ω c decreases and tend to be positive when ω c becomes inferior to 50 rad/s, which causes the system instability. To validate these results, the Simscape simulations were made to verify that any value of inferior to 50 rad/s can cause system instability. The results of the first simulation, made using the Simscape model with = 50 rad/s, are illustrated on Figure 16a, and those of the second simulation with = 45 rad/s are shown on Figure 16b. These results confirm the system stability for > 50 rad/s, as predicted by the stability analysis based on the established state model.  To validate these results, the Simscape simulations were made to verify that any value of ω c inferior to 50 rad/s can cause system instability. The results of the first simulation, made using the Simscape model with ω c = 50 rad/s, are illustrated on Figure 16a, and those of the second simulation with ω c = 45 rad/s are shown on Figure 16b. These results confirm the system stability for ω c > 50 rad/s, as predicted by the stability analysis based on the established state model. To evaluate the impact of droop control parameters on the system stability, the frequency drop ∆ is increased in 10 steps (from its initial value to 10 times this value). The eigenvalues mostly impacted by this variation are to as well as and which are two by two complex conjugate numbers. As shown in Figure 17 when ∆ increases the eigenvalues to move toward the left, while and move toward the right side of the real axis and become positive for ∆ ≥ 5 rad/s. By following the same approach as Simscape, simulations were made to verify that any value of ∆ superior or equal to 5 rad/s can cause the system instability. The results of the first simulation, made using the Simscape model with ∆ = 4.5 rad/s, are illustrated In Figure 18a, and those of the second simulation with ∆ = 5 rad/s are shown In Figure 18b. These results confirm the system instability for ∆ ≥ 5 rad/s, as predicted by the stability analysis based on the established state model.
The voltage drop ∆ is the last parameter of the DGs for which its impact on the system stability is studied. This has been done by varying ∆ from its initial value to 10 times this value. The eigenvalues mostly impacted by this variation are and , which are the complex conjugate numbers with a negative real part as well as which is a real negative number. When the voltage drop increases, Figure 19 shows that moves to the left while and move slightly to the right but remain negative. Thus, contrary to ∆ , the voltage drop ∆ has a negligible impact on the system stability. To evaluate the impact of droop control parameters on the system stability, the frequency drop ∆ ω is increased in 10 steps (from its initial value to 10 times this value). The eigenvalues mostly impacted by this variation are λ 23 to λ 26 as well as λ 33 and λ 34 which are two by two complex conjugate numbers. As shown in Figure 17 when ∆ ω increases the eigenvalues λ 23 to λ 26 move toward the left, while λ 33 and λ 34 move toward the right side of the real axis and become positive for ∆ ω ≥ 5 rad/s. By following the same approach as Simscape, simulations were made to verify that any value of ∆ ω superior or equal to 5 rad/s can cause the system instability. The results of the first simulation, made using the Simscape model with ∆ ω = 4.5 rad/s, are illustrated In Figure 18a, and those of the second simulation with ∆ ω = 5 rad/s are shown In Figure 18b. These results confirm the system instability for ∆ ω ≥ 5 rad/s, as predicted by the stability analysis based on the established state model.    The voltage drop ∆ E is the last parameter of the DGs for which its impact on the system stability is studied. This has been done by varying ∆ E from its initial value to 10 times this value. The eigenvalues mostly impacted by this variation are λ 23 and λ 24 , which are the complex conjugate numbers with a negative real part as well as λ 36 which is a real negative number. When the voltage drop increases, Figure 19 shows that λ 36 moves to the left while λ 23 and λ 24 move slightly to the right but remain negative. Thus, contrary to ∆ ω , the voltage drop ∆ E has a negligible impact on the system stability.

Eigenvalue Trajectory and Sensibility Analysis for Different Values of Load Parameters
In order to explain the effect of load variation on system stability, a study of the variation of dominant eigenvalue trajectories under different load conditions is established. Starting with the CPL type, Figure 20a,b show the trajectory of the eigenvalues under the effect of the change in active P c and reactive Q c powers absorbed by the CPL. As the active power or the reactive power absorbed by the CPL increases up to 120 kW (with fixed Q c = 0) and 72 kVAr (with fixed P c = 100 kW), respectively, the eigenvalues λ 5 to λ 8 are the most sensitive to the change of P c or Q c of the CPL. The other eigenvalues do not vary with the evolution of P c or Q c . When the active power (P c ) absorbed by the CPL increases from 100 kW up to 120 kW (with fixed Q c = 0), the eigenvalues λ 7 and λ 8 move to the left far away from zero while λ 5 and λ 6 move toward the right but remain sufficiently negative. In the same manner, when the reactive power (Q c ) absorbed by the CPL increases from zero up to 72 kVAr (with fixed P c = 100 kW), the eigenvalues λ 7 and λ 8 move to the left while λ 5 and λ 6 move toward the right. For the considered intervals of the active power and reactive power, the real part of λ 5 and λ 6 remains sufficiently negative and the system stability is not affected. However, when the active power P c or the reactive power Q c become greater than certain limits, the real parts of the eigenvalues λ 5 and λ 6 become positive which causes system instability. For example, when Q c is fixed to zero, the limit value of active power P c absorbed by the CPL is 124.5 kW, beyond this value, the real part of the eigenvalues λ 5 and λ 6 becomes positive (Figure 21a for P c = 124kW and P c = 126kW) and the system cannot be considered asymptotically stable around the considered operating point. In the same manner, when P c is fixed to 100 kW, the limit value of reactive power Q c absorbed by the CPL is 74.3 kVAr; beyond this value, the real part of the eigenvalues λ 5 and λ 6 becomes positive (Figure 21b for Q c = 73 kVAr and Q c = 74 kVAr). effect of the change in active and reactive powers absorbed by the CPL. As the active power or the reactive power absorbed by the CPL increases up to 120 kW (with fixed = 0) and 72 kVAr (with fixed = 100 kW), respectively, the eigenvalues to are the most sensitive to the change of or of the CPL. The other eigenvalues do not vary with the evolution of or . When the active power ( ) absorbed by the CPL increases from 100 kW up to 120 kW (with fixed = 0), the eigenvalues and move to the left far away from zero while and move toward the right but remain sufficiently negative. In the same manner, when the reactive power ( ) absorbed by the CPL increases from zero up to 72 kVAr (with fixed = 100 kW), the eigenvalues and move to the left while and move toward the right. For the considered intervals of the active power and reactive power, the real part of and remains sufficiently negative and the system stability is not affected. However, when the active power or the reactive power become greater than certain limits, the real parts of the eigenvalues and become positive which causes system instability. For example, when is fixed to zero, the limit value of active power absorbed by the CPL is 124.5 kW, beyond this value, the real part of the eigenvalues and becomes positive (Figure 21a for = 124kW and = 126kW) and the system cannot be considered asymptotically stable around the considered operating point. In the same manner, when is fixed to 100 kW, the limit value of reactive power absorbed by the CPL is 74.3 kVAr; beyond this value, the real part of the eigenvalues and becomes positive ( Figure  21b for = 73 kVAr and = 74 kVAr). To confirm the results concerning the stability limits imposed by the CPL, the Simscape simulations were made for two values of the CPL active power, one less than the active power limit ( = 124 kW) and another one higher than this limit ( = 125kW) and the results are shown in Figure 22a   To confirm the results concerning the stability limits imposed by the CPL, the Simscape simulations were made for two values of the CPL active power, one less than the active power limit (P c = 124 kW) and another one higher than this limit (P c = 125kW) and the results are shown in Figure 22a,b. These results are in perfect accordance with the conclusions drawn from the established state model relating to the impact of the CPL on the stability of the considered system. In order to maintain the microgrid stability for higher values of active power absorbed by the CPL, one solution consists of adding a capacitor to the PCC to which the CPL is connected. Then, the established state model can be used to dimension the necessary capacitor value for the required CPL active power. For the considered microgrid, the CPL is connected to PCC3, and based on the above analysis, the system stability is affected when the CPL active power becomes higher than 124 kW. The latter active power limit can be augmented; for example, to = 620 kW by imposing an adapted value to the capacitor connected to PCC3. This value can be determined using the established state model. When = 0, the limit value of active power absorbed by the CPL becomes 620 kW when the capacitor is modified from 0.4 µF to 2 µF. Figure 23 shows that the real part of the eigenvalues and is negative for = 620 kW and positive for = 625 kW. Then the microgrid stability is ensured up to = 620 kW. This result is in perfect accordance with the Simscape simulation results shown in Figure 24a,b. Hence, the stability analysis based on the established state model can be used as a powerful tool to In order to maintain the microgrid stability for higher values of active power absorbed by the CPL, one solution consists of adding a capacitor to the PCC to which the CPL is connected. Then, the established state model can be used to dimension the necessary capacitor value for the required CPL active power. For the considered microgrid, the CPL is connected to PCC3, and based on the above analysis, the system stability is affected when the CPL active power becomes higher than 124 kW. The latter active power limit can be augmented; for example, to P c = 620 kW by imposing an adapted value to the capacitor connected to PCC3. This value can be determined using the established state model.
When Q c = 0, the limit value of active power P c absorbed by the CPL becomes 620 kW when the capacitor c 3 is modified from 0.4 µF to 2 µF. Figure 23 shows that the real part of the eigenvalues λ 5 and λ 6 is negative for P c = 620 kW and positive for P c = 625 kW. Then the microgrid stability is ensured up to P c = 620 kW. This result is in perfect accordance with the Simscape simulation results shown in Figure 24a,b. Hence, the stability analysis based on the established state model can be used as a powerful tool to adjust mesh microgrids to guarantee their stability when their loads or even their architectures change.
Concerning the serial R-L loads in the microgrid, the impact of their variations on the system stability is studied by analyzing the evolution of the eigenvalues mostly impacted by the change under constant power factor of each RL-load connected to a given PCC. In fact, the power factor of classical loads is usually higher than 0.92 and its variation can be neglected. For the three serial RL-loads of the considered microgrid shown in Figures 10 and 25-27 show the evolution of the eigenvalues mostly impacted by varying (from initial value to 50% of this value) the impedances of RL-load3, RL-load4 and RL-load5 under constant power factor, respectively. The most impacted eigenvalues in this case are λ 1 to λ 12 , which correspond to line current equations and λ 13 to λ 20 as well as λ 25 and λ 26 , which correspond to the PCC voltage equations. It can be noticed that when the value of each RL-load decreases under constant power factor, the eigenvalues mostly impacted move to the right but remain sufficiently negative and do not impact the system stability as long as the load variation is within the indicated power range. adjust mesh microgrids to guarantee their stability when their loads or even their architectures change. Concerning the serial R-L loads in the microgrid, the impact of their variations on the system stability is studied by analyzing the evolution of the eigenvalues mostly impacted by the change under constant power factor of each RL-load connected to a given PCC. In fact, the power factor of classical loads is usually higher than 0.92 and its variation can be neglected. For the three serial RL-loads of the considered microgrid shown in Figure 10, Figures 25-27 show the evolution of the eigenvalues mostly impacted by varying (from initial value to 50% of this value) the impedances of RL-load3, RL-load4 and RL-load5 under constant power factor, respectively. The most impacted eigenvalues in this case are to , which correspond to line current equations and to as well as and , which correspond to the PCC voltage equations. It can be noticed that when the value of each RL-load decreases under constant power factor, the eigenvalues mostly impacted move to the right but remain sufficiently negative and do not impact the system stability as long as the load variation is within the indicated power range.  Concerning the serial R-L loads in the microgrid, the impact of their variations on the system stability is studied by analyzing the evolution of the eigenvalues mostly impacted by the change under constant power factor of each RL-load connected to a given PCC. In fact, the power factor of classical loads is usually higher than 0.92 and its variation can be neglected. For the three serial RL-loads of the considered microgrid shown in Figure 10, Figures 25-27 show the evolution of the eigenvalues mostly impacted by varying (from initial value to 50% of this value) the impedances of RL-load3, RL-load4 and RL-load5 under constant power factor, respectively. The most impacted eigenvalues in this case are to , which correspond to line current equations and to as well as and , which correspond to the PCC voltage equations. It can be noticed that when the value of each RL-load decreases under constant power factor, the eigenvalues mostly impacted move to the right but remain sufficiently negative and do not impact the system stability as long as the load variation is within the indicated power range.

Mesh Microgrid Control Robustness
Based on the established and validated state-model, including the modified droopcontrol that ensures active and reactive power sharing, the considered mesh microgrid control robustness with respect to the modifications caused by the connection of DGs and large load variations is studied. In this simulation, at 0 s, the first DG sets the frequency of the microgrid and the voltages of the PCCs while only the first RL load is connected. Then, the second DG is interconnected to the microgrid at 5 s after being synchronized from 1 s to 5 s. At 8 s, the second RL load is also connected, applying a high positive load step to the microgrid. Finally, at 11 s, the CPL load absorbing an active power of 100 kW is also connected to the microgrid to verify its effect on the active and reactive power sharing.
Due to the modified droop control strategy given by Equations (9) and (10), the active and reactive power sharing is ensured (Figure 28a,b). The convergence of the reactive power of the two DGs under different load conditions is convincingly verified in Figure 28b. In addition, and as foreseen by the stability analysis based on the state-model, the microgrid stability is ensured in different microgrid structure configuration when the CPL active power load does not exceed 124 kW. sharing.
Due to the modified droop control strategy given by Equations (9) and (10), the active and reactive power sharing is ensured (Figure 28a,b). The convergence of the reactive power of the two DGs under different load conditions is convincingly verified in Figure  28b. In addition, and as foreseen by the stability analysis based on the state-model, the microgrid stability is ensured in different microgrid structure configuration when the CPL active power load does not exceed 124 kW.

Discussion
An improved droop control strategy for synchronization and power sharing was applied on a mesh multi-PCC islanded microgrid. With only the information on the measured RMS voltage value at a pilot node, this strategy has proven its efficiency for fast synchronization as well as active and reactive power sharing between all microgrid DGs. This makes the mesh multi-PCC microgrid eligible to the "plug and play" feature.
In order to study the microgrid stability and its droop control robustness, a state model of the complete system has been proposed and validated based on simulation results achieved using both the Simscape model and proposed state-model. Then, the stability analysis of the considered non-linear system at each operating point was performed using the Jacobean matrix of the state model. In practice, the dynamic behavior of different state variables and the origin of different frequency components were found by studying the evolution of the system eigenvalues with respect to the parameters of the power lines, DGs and different loads. It was found that the most influent system parameters on its stability were the DGs control parameters and the CPL active power level.

Discussion
An improved droop control strategy for synchronization and power sharing was applied on a mesh multi-PCC islanded microgrid. With only the information on the measured RMS voltage value at a pilot node, this strategy has proven its efficiency for fast synchronization as well as active and reactive power sharing between all microgrid DGs. This makes the mesh multi-PCC microgrid eligible to the "plug and play" feature.
In order to study the microgrid stability and its droop control robustness, a state model of the complete system has been proposed and validated based on simulation results achieved using both the Simscape model and proposed state-model. Then, the stability analysis of the considered non-linear system at each operating point was performed using the Jacobean matrix of the state model. In practice, the dynamic behavior of different state variables and the origin of different frequency components were found by studying the evolution of the system eigenvalues with respect to the parameters of the power lines, DGs and different loads. It was found that the most influent system parameters on its stability were the DGs control parameters and the CPL active power level.
In addition, it has been shown through an example that the established state model, which allows the study of system stability, can be used as a powerful tool to adjust certain mesh microgrid parameters to guarantee its stability when the absorbed active power level of its CPLs becomes too high. Finally, the validated microgrid state-model, including the modified droop-control that ensures active and reactive power sharing, also achieved mesh microgrid control robustness with respect to the modifications caused by the connection of DGs and large load variations, which means that the new droop control is always effective under different network configurations.
The developed modified droop control strategy can be adapted and used for power sharing in grid-connected mesh microgrids, and the synchronization strategy can also be adapted for a seamless transfer from an islanded microgrid to a grid-connected microgrid. As mentioned in the Introduction, the developed and validated state-space model can be the base for many future studies in mesh microgrids.