A Novel Approach for Searching the Upper / Lower Bounds of Uncertainty Parameters in Microgrids

In this study, a novel method based on μ analysis is presented to search for the upper/lower bounds of uncertainty parameters in microgrids (MGs). It is well known that uncertainty parameters have important effects in a MG, and they may cause instability. Previous studies have mainly focused on identifying the stability of a MG with its uncertainty parameters, but they did not address the problem of the upper/lower bounds of uncertainty parameters, i.e., how far the uncertainty parameters can be extended while the system remains stable in the small-signal sense. Thus, we developed an approach for identifying the bounds of uncertainty in MGs. In the current paper, first, a method is proposed for linear fractional transformation (LFT) configuration to express the uncertainty parameters, which makes the stability of the nominal MG system independent of any extension of the bounds. An algorithm based on this configuration is then designed to find the upper/lower bounds for both single parameter and multiple uncertainty parameters in a MG. Finally, the two cases are discussed, and the accuracy of the proposed method is confirmed using the conventional eigenvalue method.


Introduction
Microgrids (MGs) have attracted much attention throughout the world because of their significant benefits for utility networks in terms of loss reduction and reliability improvement [1].Most of the generators in an MG are distributed generation (DG) units, which interface with the utility grid through an ac/dc inverter [2].However, the inverter does not have the same moment of inertia as a synchronous generator, so the problem of MG stability has attracted considerable attention [3][4][5][6][7][8][9][10][11][12][13][14][15].Previous studies regarding this issue have focused mainly on the stability of (1) systems with fixed parameters; and (2) systems with uncertainty parameters.In the former systems, the values vary in a range around the nominal value [4], and they are used to represent perturbations in the equipment parameters, load changes, and controller parameters [4][5][6][7].The major approaches used to investigate the stability of systems with fixed parameters are the eigenvalue method [8][9][10][11] and the impedance-based method [12][13][14][15][16].In systems with uncertainty parameters, the parameters can be regarded as parameters that vary within a certain range, but the aforementioned approaches are not suitable in this case because the computational burden increases considerably for a sequence of Energies 2018, 11, 1035 2 of 23 parametric changes.In particular, multiple uncertainty parameters are present in the system, and thus, µ analysis must be be employed [5][6][7][17][18][19].
Previous studies have focused on the stability of systems.However, it is also very important to establish the stability margin for a system in terms of its parameters and to calculate its range of stability.The main techniques employed to address this problem for fixed parameters are the eigenvalue method [20,21] and the bifurcation method [22][23][24][25][26].The eigenvalue method needs to determine the dominant eigenvalue as a function of the parameters.In order to reduce the computational burden associated with the eigenvalue method, modified approaches have been proposed, such as the eigenvalue sensitive algorithm [27,28] and the continuation of invariant subspaces method [29].Other techniques for addressing this problem are bifurcation methods, such as the Hopf bifurcation [22,23], saddle-node bifurcation [24,25], and structure-induced bifurcation [26].An advantage of these methods is that they can directly calculate the variable values of parameters that make the eigenvalue cross the stability boundary.
To the best of our knowledge, no previous studies have addressed this issue with respect to uncertainty parameters.In fact, uncertainty parameters have important effects on the stability of MGs, and they are considered fully in stability estimates [16][17][18][19][30][31][32][33].Therefore, it is necessary to establish the stability margin of a system in terms its the uncertainty parameters and to obtain its upper/lower bounds.It should be noted that although the super upper value of the largest structured singular value (SSV) of a system is an accurate index for determining the stability of a system in the presence of uncertainty parameters, this value cannot be treated as the stability margin with respect to the uncertainty parameters.This is because the inverse of the super upper value of the largest SSV only indicates the allowable margin of deviation from the nominal value.The system is determined as unstable regardless of whether the lower bound of the parameter decreases or the upper bound of the parameter increases.Moreover, if multiple uncertainty parameters are considered, the relationship between the inverse of the super upper values of the largest SSV and the stability boundaries of the parameters will be much more complex.
In the present study, we developed an approach to search for the bound of an uncertainty parameter in an MG.In the current paper, first, a linear fractional transformation (LFT) method is introduced to express an uncertainty parameter.The frameworks for the DG unit and MG are then constructed.We also propose an algorithm for identifying the upper/lower bounds of an uncertainty parameter.The proposed method is used to study the parameters for a droop controller in an MG.The selection of these parameters is based on the following consideration: the parameters for a droop controller are determined during the operation of the system, but the parameters have no defined values and they are only are subject to a few constraints [4], and thus, different designers may select different values.Therefore, the droop parameter can be treated as an interval variable, and it is necessary to evaluate the system's margins in terms of different parameters.Finally, in order to provide a comprehensive evaluation of the proposed method, two case studies are performed.In the first case study, we search for the upper/lower bounds of a single uncertainty parameter.In the second case study, we demonstrate the capacity of the proposed approach to identify the stability boundaries for multiple uncertainty parameters.The stability boundaries obtained are confirmed using the eigenvalue method.
The main contributions of this study are summarized as follows: (1) An algorithm is proposed to search for the upper/lower bounds of uncertainty parameters in MGs.In addition, a method is proposed for LFT to configure the uncertainty parameters, which has general applicability to the determination of uncertainty parameters.(2) The method for configuring DG units is described in detail, and it can readily be extended to an MG with multiple sources.Based on theoretical procedures, we demonstrate how to study the uncertainty parameters in other MG structures.(3) If this method is applied correctly, it is possible to accurately identify the stability boundaries for uncertainty parameters.

Description of an MG and Its Controller
Figure 1 shows a schematic illustrating an MG and its control structure.The MG comprises two electronically-coupled DG units and a three-wire balanced load.In each DG unit, the inverter is supplied by a DC voltage source, and it is integrated into the MG through an LCL lower pass filter and inductive wire.During normal operation, the MG is connected to the utility through an intelligent bypass switch (IBS), and it supplies power together with the grid.When the grid is absent, the IBS is open and the MG enters the islanded mode, which supports the frequency and voltage of the system.The control strategy for the DG unit in island mode is based on three loops [31], i.e., an outer loop control shares the load by employing the droop strategy, a middle loop follows the reference voltage using a proportional-integral regulator, and an inner loop enhances the system's dynamic response capacity via a proportion controller.

Description of an MG and Its Controller
Figure 1 shows a schematic illustrating an MG and its control structure.The MG comprises two electronically-coupled DG units and a three-wire balanced load.In each DG unit, the inverter is supplied by a DC voltage source, and it is integrated into the MG through an LCL lower pass filter and inductive wire.During normal operation, the MG is connected to the utility through an intelligent bypass switch (IBS), and it supplies power together with the grid.When the grid is absent, the IBS is open and the MG enters the islanded mode, which supports the frequency and voltage of the system.The control strategy for the DG unit in island mode is based on three loops [31], i.e., an outer loop control shares the load by employing the droop strategy, a middle loop follows the reference voltage using a proportional-integral regulator, and an inner loop enhances the system's dynamic response capacity via a proportion controller.Details of the mathematical derivations and the sign convention are explained in Appendix A.  As shown in Figure 2, an analogue block diagram of a DG unit and the control loop are built in the d − q rotating reference frame.In the block diagram, the currents of two inductors, îLdq and îodq , the voltage of the capacitor, vodq , the outputs average power, P and Q, and the virtual variable, Φdq are taken as the state variables.The external input and output variables are îodq and vbdq , respectively.Details of the mathematical derivations and the sign convention are explained in Appendix A. As shown in Figure 2, an analogue block diagram of a DG unit and the control loop are built in the d q − rotating reference frame.In the block diagram, the currents of two inductors, ˆLdq i and ˆodq i , the voltage of the capacitor, ˆodq v , the outputs average power, P and Q , and the virtual variable, ˆdq Φ are taken as the state variables.The external input and output variables are ˆodq i and ˆbdq v , respectively.
Details of the mathematical derivations and the sign convention are explained in Appendix A.
Analogue block diagram of a distributed generation (DG) unit and its controller.

The Impact of Droop Parameters on Microgrid Stability
To consider the impact of the droop parameters on the stability of the MG, an eigenvalue analysis of the MG was conducted based on the small-signal model in Figure 1. Figure 3 shows the locus plots of the dominant eigenvalues as functions of the frequency droop parameter.Clearly, as the frequency droop parameter increases, the eigenvalues move toward the unstable region, which indicates that the stability margin of the MG is reduced.If the eigenvalues are sufficiently close to the positive axis, then even small fluctuations in the droop parameters can have considerable effects on the system's stability.Large droop parameters can make the MG more oscillatory and lead to instability, but they may effectively improve the transient response of the DG unit [2].In general, designers balance the tradeoff between the transient response and the stability when setting droop parameters.Therefore, considering that droop parameters play important roles in a system's stability, and because of the system's actual characteristics, it is necessary to study the stability boundaries for the droop parameters.

The Impact of Droop Parameters on Microgrid Stability
To consider the impact of the droop parameters on the stability of the MG, an eigenvalue analysis of the MG was conducted based on the small-signal model in Figure 1. Figure 3 shows the locus plots of the dominant eigenvalues as functions of the frequency droop parameter.Clearly, as the frequency droop parameter increases, the eigenvalues move toward the unstable region, which indicates that the stability margin of the MG is reduced.If the eigenvalues are sufficiently close to the positive axis, then even small fluctuations in the droop parameters can have considerable effects on the system's stability.Large droop parameters can make the MG more oscillatory and lead to instability, but they may effectively improve the transient response of the DG unit [2].In general, designers balance the tradeoff between the transient response and the stability when setting droop parameters.Therefore, considering that droop parameters play important roles in a system's stability, and because of the system's actual characteristics, it is necessary to study the stability boundaries for the droop parameters.

Methodology
3.1.μ Analysis of the Stability of a System with Uncertainty Parameters μ analysis [32] is a method used for examining whether a structured perturbation system is stable.Figure 4 shows the standard M − Δ configuration of a structured perturbation system, where

µ Analysis of the Stability of a System with Uncertainty Parameters
µ analysis [32] is a method used for examining whether a structured perturbation system is stable.Figure 4 shows the standard M − ∆ configuration of a structured perturbation system, where M(s) denotes the transfer function of a nominal system, and ∆ represents the "parametric uncertainty" of the system, and both are stable and bounded.It is assumed that the block structure, ∆, has been normalized, i.e., ∆ ∞ < 1.According to the definition of the µ analysis, if and only if the inequality relationship, µ ∆ (M(s)) < 1, is satisfied, then the structured perturbation system is stable with respect to ∆, where µ ∆ (M(s)) := sup ω∈ µ ∆ (M(jω)) represents the super upper value of the largest SSV of µ ∆ (M(jω)) in the frequency domain.In contrast, if µ ∆ (M(s)) > 1, then the structured perturbation system is unstable with respect to ∆.Therefore, the super upper value of the largest SSV can be treated as an index for identifying the stability of a perturbation system.
It should be noted that two limitations of the µ analysis need to be addressed: the system needs to be represented as the M − ∆ configuration before the µ analysis, and a stable M(s) is the premise for the µ analysis.These limitations constrain the configuration of M, and they are fully considered in the proposed algorithm.In order to satisfy these limitations, we introduced a method to build the LFT configuration for the uncertainty parameters.
Energies 2018, 11, x 5 of 24 considered in the proposed algorithm.In order to satisfy these limitations, we introduced a method to build the LFT configuration for the uncertainty parameters.
 is expressed as [31]: where As explained previously, the parameters m and n are treated as interval variables in this study.In general, a parameter, c , that is bounded in the region [ ] , a b can be represented as a perturbation of its nominal value, 0 c [32]: where 0 ( )/2 c a b . In Equation ( 2), the perturbation of a parameter is described in two parts, where c δ denotes the mathematical perturbation and c λ is the named size coefficient representing the size of the perturbation.Thus, the frequency droop parameter, m , and the voltage droop parameter, n , can be expressed as Equations ( 3) and ( 4), respectively:

Proposed LFT Configuration for Uncertainty Parameters
Figure 2 shows a block diagram of a DG unit.In general, an MG employs the droop strategy to share the load in the islanded mode.Let P and Q be the terminal active power and reactive power of a DG unit, respectively, and v * o and ω are the output reference voltage and frequency of the droop controller.The transfer function of droop controller, G D , from input is expressed as [31]: where As explained previously, the parameters m and n are treated as interval variables in this study.In general, a parameter, c, that is bounded in the region [a, b] can be represented as a perturbation of its nominal value, c 0 [32]: where In Equation ( 2), the perturbation of a parameter is described in two parts, where δ c denotes the mathematical perturbation and λ c is the named size coefficient representing the size of the perturbation.Thus, the frequency droop parameter, m, and the voltage droop parameter, n, can be expressed as Equations ( 3) and ( 4), respectively: where m 0 and n 0 are the nominal values of parameters m and n, respectively, and λ m and λ n are the named size coefficients which adjust the size of the interval for each parameter.After substituting Equations ( 3) and ( 4) into Equation ( 1), the transfer function, G D , is as follows: By employing the LFT method, the perturbations in the transfer function, G D , are extracted and arranged in a diagonal matrix, ∆.Correspondingly, the transfer function, G D , is expanded and the fictitious vectors are introduced as the increased input and output vectors.The configuration for G D is depicted in Figure 5. Energies 2018, 11, x 6 of 24 - where It should be noted that the matrix D G can be decomposed in different ways, and the matrices where It should be noted that the matrix G D can be decomposed in different ways, and the matrices Q 11 , Q 12 , Q 21 and Q 22 have other solutions.In this study, all of the size coefficients were arranged in block Q 12 where the output was the fictitious vector, e.This arrangement ensured that the stability of the nominal system was independent of the coefficients, λ m and λ n , as explained in Section 3.3.
It should be noted that although this model is similar to the small-signal models presented Mohamed and El-Saadany [3], the difference between the two models is obvious because the uncertainty is integrated into the droop parameters of the DG unit in the present study, which changes the configuration of the DG unit.The configuration for the DG unit with the uncertainty parameter is shown in Figure 6, where G DGi , represented by Equations ( 7)-( 9), denotes the transfer function matrix of the nominal model of the DG unit.The perturbation block, ∆ d , represents the variations in the droop parameters.2_ 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 It should be noted that although this model is similar to the small-signal models presented Mohamed and El-Saadany [3], the difference between the two models is obvious because the uncertainty is integrated into the droop parameters of the DG unit in the present study, which changes the configuration of the DG unit.The configuration for the DG unit with the uncertainty parameter is shown in Figure 6, where DGi G , represented by Equations ( 7)-( 9  To obtain the small-signal model for the MG, the DG unit described by Equations ( 7)-( 10) is combined with those for the load and network.This combination is shown in Figure 7, and the mathematical derivations are presented in Appendix B. The complete state-space expressions for the MG are as follows: In Equations ( 11)-(13), the system state matrices, A con , B con , C con and ∆, are as follows:  11)-( 13) is shown in Figure 8.
It should be noted that coefficients λ m and λ n are excluded from matrix A con due to the proposed arrangement for the uncertainty parameter.In fact, after LFT divides the perturbation function, G D , into the nominal function, Q, and uncertainty blocks, ∆, the fictitious output vector, e, of the nominal function does not participate in the calculation of the state equations of the DG system Furthermore, Figure 7 shows that the fictitious vector, e DG , which comprises all of the fictitious vectors, e, in each DG unit still does not participate in the calculation of the state equations for the nominal system, G MG .Therefore, the stability of the nominal system, G MG , is independent of coefficients λ m and λ n under this configuration.The proposed configuration accurately represents the system with parameter variations, and it is suitable for the subsequent µ analysis.

Stable Upper/Lower Bound Analysis for a Single Uncertainty Parameter
After the MG with uncertainty parameters has been described using Equations ( 11)-( 13), the value ( ) can reflect the stability of the MG with respect to a given parameter within a certain range.For convenience, we consider parameter c instead of parameter m or n in the following explanation.First, we discuss the case with only one uncertainty parameter in an MG.According to the μ definitions, if the value of

Stable Upper/Lower Bound Analysis for a Single Uncertainty Parameter
After the MG with uncertainty parameters has been described using Equations ( 11)-( 13), the value ( ) can reflect the stability of the MG with respect to a given parameter within a certain range.For convenience, we consider parameter c instead of parameter m or n in the following explanation.First, we discuss the case with only one uncertainty parameter in an MG.According to the μ definitions, if the value of

Stable Upper/Lower Bound Analysis for a Single Uncertainty Parameter
After the MG with uncertainty parameters has been described using Equations ( 11)-( 13), the value µ ∆ (G MG ) can reflect the stability of the MG with respect to a given parameter within a certain range.
For convenience, we consider parameter c instead of parameter m or n in the following explanation.First, we discuss the case with only one uncertainty parameter in an MG.According to the µ definitions, if the value of µ ∆ (G MG (c 0 , λ c )) is less than 1, the MG is stable and parameter c takes any value within the range of [c 0 − λ c , c 0 + λ c ].In contrast, when µ ∆ (G MG (c 0 , λ c )) is greater than 1, this indicates that this range is greater than the stability boundary for parameter c.Based on this fact, we propose the following method for identifying the upper/lower bounds of an uncertainty parameter.First, the value c 0 is initiated as the operating value for the MG and the size coefficient, λ c_0 , is set.The size coefficient, λ c_0 , is then increased in steps of t λ to extend the interval of parameter c.The iterative process is completed up to a value of λ c_n , which makes µ ∆ (G MG (c 0 , λ c_i )) greater than 1.As shown in Figure 9a, two possible options can trigger µ ∆ (G MG (c 0 , λ c_n )) to exceed 1: the extension of the interval to reach the lower bound of parameter c if the initial value, c 0 , is near the lower bound; or if the value of c 0 is close to the upper bound, the upper bound of parameter c is reached first.  .This selection will make the error no more than step value t λ .For the other bound of uncertainty parameter c , we first consider the case where the previous trigger is the upper bound of parameter c .To identify the lower bound of c , the size coefficient is fixed at  To determine the most appropriate possible option, the size coefficient is fixed at λ c_n , and the nominal value, c 0 , is reduced by a step, t λ .If the value of µ ∆ (G MG (c 0 − t λ , λ c_n )) changes to less than 1, the trigger is the upper bound of the parameter, and it must fall within the interval of [c 0 + λ c_n−1 , c 0 + λ c_n ]; otherwise, if the value of µ ∆ (G MG (c 0 − t λ , λ c_n )) remains greater than 1, the trigger is the lower bound, where the value is limited within [c 0 − λ c_n , c 0 − λ c_n−1 ].Theoretically, the size coefficient that makes µ ∆ (G MG (c 0 , λ c_i )) closest to 1 denotes the maximum size of the perturbation of parameter c.However, it is difficult to satisfy this condition because it incurs considerable computational burden.Thus, it is replaced by the maximum size coefficient to make µ ∆ (G MG (c 0 , λ c_i )) no more than 1.For instance, λ c_n−1 is selected as the maximum size of the perturbation, and the upper bound of c is taken as c 0 + λ c_n−1 .This selection will make the error no more than step value t λ .
For the other bound of uncertainty parameter c, we first consider the case where the previous trigger is the upper bound of parameter c.To identify the lower bound of c, the size coefficient is fixed at λ c_n according to the steps shown in Figure 9a, and the nominal value, c 0 , is reduced in steps of t c .As shown in Figure 9b, as c 0 decreases, the interval [c j − λ c_n , c j + λ c_n ] slides to the left so the system enters a stable state again.Up to the point where the lower bound of the uncertainty parameter is included in the interval, the value of µ ∆ (G MG (c m , λ c_n )) will change to exceed 1; the lower bound of the parameter is subsequently within the interval of [c m − λ c_n , c m + λ c_n ].To further identify the accurate point for the lower bound, c m is then fixed and λ c_n decreases in steps of t λ to narrow the interval.According to the theoretical µ analysis, if the lower bound is outside the range of [c m − λ c_k , c m + λ c_k ], then the value of µ ∆ (G MG (c m , λ c_k )) becomes smaller than 1.Hence, the lower bound of the parameter can be taken as c m − λ c_k .Similarly, if the previous trigger is the lower bound of the parameter, the upper bound of parameter c can also be obtained in this manner.This process is shown in Figure 9c.
The key to this method is satisfying the requirement of the µ analysis that the nominal system, G MG (c j , λ c_i ), should be stable.This requirement can be classified with two cases, as follows.In the first case, we consider how changing the size coefficient, λ c_i , affects the stability of G MG (c j , λ c_i ) when the nominal parameter, c j , is fixed.Clearly, the stability of nominal system G MG (c j , λ c_i ) is determined by matrix A con of system state space expressions.However, in the proposed method, coefficients λ m and λ n are not introduced into matrix A con , and thus the internal stability of the nominal system, G MG (c j , λ c_i ), remains the same as the steady state regardless of how the size coefficients change.
In the second case, we consider that the size coefficient is fixed at λ c_n , and we need to ensure the stability of nominal system G MG (c j+1 , λ c_n ) when the nominal parameter varies from c j to c j+1 .According to the µ analysis, if the size coefficients have no impact on the stability of the nominal system, then the conditions in which the perturbation system is stable with respect to c in the range of c j − λ c_n , c j + λ c_n can ensure that the nominal system G MG (c j , λ c_n ) is stable with respect to parameter c j in the same range.Thus, this condition is satisfied provided that step t c does not exceed step λ c_n .In the proposed method, the value of t c is set as equal to λ c_n .

Analysis of Stability Boundaries for Two Uncertainty Parameters
In this problem, we note that for an uncertainty parameter, c1, there is a stability boundary, c2, for the uncertainty parameter under the MG, as shown in Figure 10.The stability boundaries of the two uncertainty parameters comprise many of these areas.Therefore, obtaining the stability boundary, c2, for a given uncertain parameter, c1, is crucial for solving this problem.Therefore, we assume that parameter c1 can be expressed as where By substituting c1 0 and λ c1 from Equation ( 14) into Equations ( 11)-( 13), the problem of searching for stable boundaries for c1 and c2 can be simplified as identifying the stability boundary of a single parameter, c2.The details of the algorithm are given in Figure 11, and the search process is depicted in Figure 12, where the vertical axis denotes the predetermined range of c1, and the horizontal axis represents the range of c2.boundary, 2 c , for a given uncertain parameter, 1 c , is crucial for solving this problem.Therefore, we assume that parameter 1 c can be expressed as where  λ from Equation ( 14) into Equations ( 11)-( 13), the problem of searching for stable boundaries for 1 c and 2 c can be simplified as identifying the stability boundary of a single parameter, 2 c .The details of the algorithm are given in Figure 11, and the search process is depicted in Figure 12, where the vertical axis denotes the predetermined range of

Study Cases
Next, we present examples to demonstrate how the proposed method can be used to obtain the upper/lower bounds of uncertainty parameters in the MG shown in Figure 1.The nominal values of the system parameters are given in Table 1.Two cases were considered in this study: (1) searching for the stability boundary for a single uncertainty parameter; and (2) searching for two uncertainty parameters in the MG.The eigenvalues of the uncertain MG under the parameters in the boundary obtained were calculated to verify the result.In this case, we assumed that the system in Figure 1 was subject to a single uncertainty parameter, i.e., parameter m or n, and the proposed method was employed to identify the stable boundary for each parameter of interest.It should be noted that the constraints on the parameters were considered, i.e., the values of m DG1 and m DG2 were limited by the rated active power ratio, ε 1 , of DG1 to DG2.Similarly, the values of n DG1 and n DG2 were limited by the reactive rated power ratio, ε 2 , of DG1 to DG2.Thus, the size coefficients of the two DG units were set to λ m_DG1 = ε 1 λ m_DG2 and λ n_DG1 = ε 2 λ n_DG2 .First, we considered the parameter m to be the uncertainty in the MG.To illustrate the search process using the proposed method, the value of µ ∆ (G MG ) as a function of the size coefficient λ m_DG1 is plotted in Figure 13a, where m 0 = 22 × 10 −4 , λ m_0 = 0.06, and t λ = 0.01.As the size coefficient, λ m_DG1 , increased, the value of µ ∆ (G MG ) increased gradually and became close to 1, thereby indicating that the robust margin of the system was reduced.When the size factor, λ m_DG1 , had a value of 0.17, µ ∆ (G MG ) was able to reach 0.9629.As λ m_DG1 continued to increase by one step, the value of µ ∆ (G MG ) was 1.018, which showed that the system was unstable.Therefore, the maximum value of λ m_DG1 at which the system remained stable was 0.17.By substituting this value into Equation (3), we derived the upper bound of the frequency droop parameter m DG1 as 25.74 × 10 −4 .Thus, the system is stable in the small-signal sense provided that the parameter m DG1 is smaller than 25.74 × 10 −4 .
To verify the result, the eigenvalues were calculated for the MG system with the parameter, m DG1 , under this stability range, and they are plotted in Figure 13b.The initial conditions of the system were the same as those described above.As shown, as parameter m DG1 increased from 0 to 25.8 × 10 −4 , some eigenvalues moved toward the unstable region, and only a few eigenvalues were far from the positive real axis.Clearly, all of the dominant eigenvalues fell into the left half plane.However, as the value of m DG1 increased to 26.2 × 10 −4 , a pair of eigenvalues crossed the imaginary axis.Thus, the stability boundary of m DG1 can be taken as 25.8 × 10 −4 , which is close to the computed boundary given in Figure 13a, although it is not equal.
axis.Thus, the stability boundary of    To further illustrate function µ ∆ (G MG (s)), the frequency responses of µ ∆ (G MG (jω)) were plotted by considering the droop parameters, m DG1 , of 25.08 and 25.96.For the first value, µ ∆ (G MG (jω)) was smaller than 1 for the entire frequency domain.In contrast, when parameter m DG1 was set to 25.96, Figure 14 clearly shows that the maximum value of µ ∆ (G MG (jω)) was 1.188 at a frequency of approximately 88 rad/s.Thus, the system was unstable under this value.Next, we considered parameter n as the uncertainty in the MG system.Figure 15a shows the trajectory of µ ∆ (G MG ) as a function of size coefficient λ n_DG1 .The initial values for n 0 and λ n_0 were 1 × 10 −3 and 0.1, respectively, and the step size, t λ , was 0.05.As shown, using the maximum size factor, λ n_DG1 , the value of µ ∆ (G MG ) was smaller than 1 at 0.75.By substituting the value of λ n_DG1 into Equation (4), we derived the upper bound of the voltage droop parameter, n DG1 , as 1.75 × 10 −3 .
In addition, we calculated the eigenvalues of the MG system with the parameter, n DG1 , within this range to verify the result.Figure 15b shows the traces of the dominant eigenvalues when the voltage droop parameter, n DG1 , varied from 0 to 1.8 × 10 −3 .When the droop parameter, n DG1 , was not greater than 1.7 × 10 −3 , all of the dominant eigenvalues fell into the left half plane.An eigenvalue crossed the imaginary axis when the droop parameter, n DG1 , was 1.8 × 10 −3 , which is close to the stability range obtained using the proposed method.

Case 2: Analysis of the Stability Range for Two Uncertainty Parameters
In this case, the system in Figure 1 was subject to simultaneous variations in two parameters, m and n .The proposed method was employed to identify the stable boundaries for these parameters.
As explained previously, we first addressed the problem of how to obtain the stability boundary of uncertainty parameter 2 c when 1 c varied within a certain range.For example, the uncertainty parameter, DG m , had an interval of  To verify the region obtained, we considered as many points ( )

Case 2: Analysis of the Stability Range for Two Uncertainty Parameters
In this case, the system in Figure 1 was subject to simultaneous variations in two parameters, m and n.The proposed method was employed to identify the stable boundaries for these parameters.As explained previously, we first addressed the problem of how to obtain the stability boundary of uncertainty parameter c2 when c1 varied within a certain range.For example, the uncertainty parameter, m DG , had an interval of 19.8 × 10 −4 , 24.2 × 10 −4 ; then, using the proposed method, the stability boundary for parameter n DG was determined as 0, 1.3 × 10 −3 .Thus, the system is stable if (m DG , n DG ) is in the closed area shown in Figure 16.

Case 2: Analysis of the Stability Range for Two Uncertainty Parameters
In this case, the system in Figure 1 was subject to simultaneous variations in two parameters, m and n .The proposed method was employed to identify the stable boundaries for these parameters.
As explained previously, we first addressed the problem of how to obtain the stability boundary of uncertainty parameter 2 c when 1 c varied within a certain range.For example, the uncertainty parameter, DG m , had an interval of  To verify the region obtained, we considered as many points (m DG , n DG ) as possible in this region and investigated the stability of the MG system at each point.Without any loss of generality, m was set to 19.8 × 10 −4 , 22 × 10 −4 , or 24.2 × 10 −4 , and n was taken as being within 0 ≤ n ≤ 1.3 × 10 −3 .The eigenvalues of an MG system considering such values are depicted in Figure 17, which demonstrates that the system remained stable when m and n were within this range.Outside this range, such as when the parameter n increased slightly, an eigenvalue crossed the imaginary axis immediately.Similarly, when parameter m increased slightly, a pair of eigenvalues moved to cross the imaginary axis.
immediately.Similarly, when parameter m increased slightly, a pair of eigenvalues moved to cross the imaginary axis.λ instead of m and n , respectively, to plot the stability range.Figure 18 shows the stability bounds for m λ and n λ , where each point on this black line represents the upper bound of parameter n λ corresponding to each value of m λ .For example, when the size coefficient, m λ , had a value of about 0.08, the maximum size coefficient, n λ , was 0.4.After converting these values to m and n using Equations ( 3) and ( 4), respectively, we found that when the uncertainty of parameter m was smaller than   We then addressed the problem of the stability ranges for m and n.As explained previously, the stability ranges comprised multiple areas, as shown in Figure 10.To clearly illustrate the relationship between the stability boundaries of m and n, we employed the size coefficients λ m and λ n instead of m and n, respectively, to plot the stability range.Figure 18 shows the stability bounds for λ m and λ n , where each point on this black line represents the upper bound of parameter λ n corresponding to each value of λ m .For example, when the size coefficient, λ m , had a value of about 0.08, the maximum size coefficient, λ n , was 0.4.After converting these values to m and n using Equations ( 3) and ( 4), respectively, we found that when the uncertainty of parameter m was smaller than 23.76 × 10 −4 , the system was stable provided that the uncertainty of parameter n was smaller than 1.4 × 10 −3 .λ instead of m and n , respectively, to plot the stability range.Figure 18 shows the stability bounds for m λ and n λ , where each point on this black line represents the upper bound of parameter n λ corresponding to each value of m λ .For example, when the size coefficient, m λ , had a value of about 0.08, the maximum size coefficient, n λ , was 0.4.After converting these values to m and n using Equations ( 3) and ( 4), respectively, we found that when the uncertainty of parameter m was smaller than

Simulation Study
To further demonstrate the validity of the stability boundaries for the droop parameters, a model of the MG considered in this study was simulated with PSCAD/EMTDC software (Version:4.2.1, Manitoba HVDC Research Centre, Winnipeg, MB, Canada).A detailed illustration of the inverter circuit, pulse-width modulation (PWM) schemes, and power control is shown in Figure 1.The parameters and initial conditions for the system are listed in Table 1.
Two size coefficients for the droop parameters were selected from the values shown in Figure 18, λ m = 0.1 and λ n = 0.2, which were in the stability range, and λ m = 0.13 and λ n = 0.3, which were outside the range.The values of the size coefficients were converted into m and n using Equations ( 3) and ( 4), respectively.Figure 19 shows the dynamics of the output voltages and currents for the two sets of droop parameters.The system started from a steady state where the output currents and voltages for both DG1 and DG2 were close to constant.At t = 1.02 s, the first set of data was assigned to m and n.After t = 1.1 s, the other set of data was assigned.Figure 19 shows that the system remained stable for the first set of droop parameters.However, for data outside these points, the instantaneous currents and voltages oscillated for more than 10 cycles, and eventually became unstable.Therefore, the simulation results demonstrate the accuracy of the stability boundaries of the droop parameters predicted by µ analysis.
Energies 2018, 11, x 18 of 24 To further demonstrate the validity of the stability boundaries for the droop parameters, a model of the MG considered in this study was simulated with PSCAD/EMTDC software (Version:4.2.1, Manitoba HVDC Research Centre, Winnipeg, MB, Canada).A detailed illustration of the inverter circuit, pulse-width modulation (PWM) schemes, and power control is shown in Figure 1.The parameters and initial conditions for the system are listed in Table 1.
Two size coefficients for the droop parameters were selected from the values shown in Figure 18 3) and (4), respectively.
Figure 19 shows the dynamics of the output voltages and currents for the two sets of droop parameters.The system started from a steady state where the output currents and voltages for both DG1 and DG2 were close to constant.At 1.02 t = s, the first set of data was assigned to m and n .After 1.1 t = s, the other set of data was assigned.Figure 19 shows that the system remained stable for the first set of droop parameters.However, for data outside these points, the instantaneous currents and voltages oscillated for more than 10 cycles, and eventually became unstable.Therefore, the simulation results demonstrate the accuracy of the stability boundaries of the droop parameters predicted by μ analysis.The system state matrices, A net , B 1net , and B 2net , are as follows: where îloadDQ is the current of the load.
To connect each subsystem model for the MG, a virtual resistor, R vir , is assumed between each node and the ground.Therefore, the voltage of each node is given by: where 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 Based on a common reference frame, the complete state-space small-signal model of the MG can be obtained by combining the state-space models of the DG with the RL load and network.

Figure 2 .
Figure 2. Analogue block diagram of a distributed generation (DG) unit and its controller.

Figure 3 .
Figure 3. Trace of dominant eigenvalue with increasing the frequency droop gain of both DG1 and DG2 units.

Figure 3 .
Figure 3. Trace of dominant eigenvalue with increasing the frequency droop gain of both DG1 and DG2 units.

Figure 4 .
Figure 4. Standard M − Δ configuration of a system with perturbations.

3. 2 .
Figure 2 shows a block diagram of a DG unit.In general, an MG employs the droop strategy to share the load in the islanded mode.Let P and Q be the terminal active power and reactive power of a DG unit, respectively, and * o v and ω are the output reference voltage and frequency of the droop controller.The transfer function of droop controller, D G , from input ˆT w P Q   =   to output

Figure 4 .
Figure 4. Standard M − ∆ configuration of a system with perturbations.

Figure 5 .
Figure 5. Representation of droop parameters as the linear fractional transformation (LFT) configuration.

Figure 5 .
Figure 5. Representation of droop parameters as the linear fractional transformation (LFT) configuration.

T
are the state vectors of the DG unit.The input vectors are d DGi = d m d n T and u DGi = vbdi vbqi T , and the output vectors are e DGi = e m e n T and y DGi = îodi îoqi T .ωcom is the angular frequency of the common reference frame, and θi represents the angle between the individual reference frame and the common reference frame.The matrices for A DGi , B 1_DGi , B 2_DGi , B 3_DGi , C 1_DGi , and C 2_DGi are as follows: ), denotes the transfer function matrix of the nominal model of the DG unit.The perturbation block, d Δ , represents the variations in the droop parameters.

Energies 2018, 11 , x 9 of 24 Figure 7 .
Figure 7. Block diagram of the complete state-space model of a MG.

Figure 8 .
Figure 8. M − Δ configuration of a MG with uncertainty parameters.

Figure 7 . 24 Figure 7 .
Figure 7. Block diagram of the complete state-space model of a MG.

Figure 8 .
Figure 8. M − Δ configuration of a MG with uncertainty parameters.

Figure 8 .
Figure 8. M − ∆ configuration of a MG with uncertainty parameters.
trigger is the upper bound of the parameter, and it must fall within the interval of the maximum size of the perturbation of parameter c .However, it is difficult to satisfy this condition because it incurs considerable computational burden.Thus, it is replaced by the maximum size coefficient to make the maximum size of the perturbation, and the upper bound of c is taken as 0

Figure 9 .
Figure 9. Process employed to identify the upper/lower bounds of a single uncertainty parameter.(a) finding a boundary for an uncertainty parameter and determining its type; (b) searching for the lower bound; (c) searching for the upper bound.

Δ
steps shown in Figure9a, and the nominal value, 0 c , is reduced in steps of c t .As shown in Figure9b, as 0 c decreases, the interval enters a stable state again.Up to the point where the lower bound of the uncertainty parameter is included in the interval, the value of will change to exceed 1; the lower bound of the parameter is subsequently within the interval of

Figure 9 .
Figure 9. Process employed to identify the upper/lower bounds of a single uncertainty parameter.(a) finding a boundary for an uncertainty parameter and determining its type; (b) searching for the lower bound; (c) searching for the upper bound.

1 c
, and the horizontal axis represents the range of 2 c .

Figure 11 .
Figure 11.Algorithm for identifying the bounds of multiple uncertainty parameters.

Figure 11 .
Figure 11.Algorithm for identifying the bounds of multiple uncertainty parameters.
close to the computed boundary given in Figure13a, although it is not equal.

Figure 14 .
Figure 14.Frequency responses of
close to the computed boundary given in Figure13a, although it is not equal.

Figure 14 .
Figure 14.Frequency responses of

Figure 16 .
Figure 16.The stable region for droop parameters in the MG.

Figure 16 .Figure 16 .
Figure 16.The stable region for droop parameters in the MG.To verify the region obtained, we considered as many points ( ), DG DGm n as possible in this region and investigated the stability of the MG system at each point.Without any loss of generality, m was set to

Figure 17 .
Figure 17.Dominant eigenvalues of a MG when

4 23 .
76 10 − × , the system was stable provided that the uncertainty of parameter n

Figure 18 .
Figure 18.The stable region for size coefficients m λ and n λ .

Energies 2018 ,
11, x 17 of 24immediately.Similarly, when parameter m increased slightly, a pair of eigenvalues moved to cross the imaginary axis.

Figure 17 ..
Figure 17.Dominant eigenvalues of a MG when

4 23 .
76 10 − × , the system was stable provided that the uncertainty of parameter n

Figure 18 .
Figure 18.The stable region for size coefficients m λ and n λ .

Figure 18 .
Figure 18.The stable region for size coefficients λ m and λ n .
outside the range.The values of the size coefficients were converted into m and n using Equations (

Figure 19 .
Figure 19.Response of the DG unit with two sets of droop parameters: (a) d-axis components of the output voltage, v o ; (b) q-axis components of the output voltage, v o ; (c) d-axis components of the output current, i o ; and (d) q-axis components of the output current, i o .

B
2net = I line1Q −I line1D I line2Q −I line2D T , where R line1 and L line1 are the resistor and inductor of the line, respectively.In this study, only a general RL load is considered.The state equation for the RL load is as follows:ˆ.i loadDQ = A load [ îloadDQ ] + B 1load [ vbDQ ] + B 2load [ ωcom ].
A net , B 1_net , B 2_net , M DG , M NET , M LOAD , T S1 , T S2 , T −1 V1 , T −1V1 , T C2 , and T C1 are given in the Appendix B. The M − ∆ configuration for Equations ( where
4.1.Case 1: Analysis of the Upper/Lower Stability Bounds for a Single Uncertainty Parameter