Virtual Impedance-Based Droop Control Scheme to Avoid Power Quality and Stability Problems in VSI-Dominated Microgrids

Droop-based techniques have widely been utilized in controlling voltage source inverter based distributed generation units within a microgrid. Advantages of droop based control have been extensively discussed, while its effects on stability and power quality of the inverter based microgrids have not been explicitly studied. This paper addresses this gap by analyzing the stability of islanded microgrid with focusing on power quality problems. From this viewpoint, current paper studies the effect of parameters change on the stability margin of the microgrid using a bifurcation analysis approach. To this end, this work firstly derives a dynamic model of a microgrid including controllers of the inverter, network and loads. Based on this model, the loci of the system eigenvalues have been analyzed, showing that, by changing the droop coefficients, the hopf and period doubling bifurcation can occur in such a system. Furthermore, the droop tuning procedure may also cause the system to experience the flicker phenomenon in advance of the instability. Next, a virtual impedance-based solution is presented to avoid the flicker as well as instability problems. Finally, the efficacy of the proposed method is proved on a practical network.


INDEX TERMS
The growing interest and increasing penetration levels of renewable Distributed Generation (DG) have breathed a new life into the Micro-Grid (MG) concept as an accessible practical solution to effectively manage distribution networks (DNs) [1], [2]. The renewable DGs are usually interfaced with the grid via the voltage source inverter (VSI) being responsible for converting the DC output to a smooth sinusoidal AC power [3], [4]. The VSI also needs a controller to regulate the output power corresponding to the network requirements [5], [6]. This controller empowers the renewable-based MGs to have higher flexibility against new concerns and challenges VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ in voltage and frequency stabilities of the islanded network. In other words, normally, under grid-connected mode, there is no significant concern about the voltage or frequency of the MG, i.e. the up-stream network dictates both frequency and voltages of the MG due to its semi-infinity nature and higher short-circuit capacity of its buses, respectively [7], [8]. However, in isolated MGs (islanded mode), any power inequality will lead to frequency or voltage instabilities [9]. Consequently, VSI controllers should additionally monitor both frequency and voltage to ensure reliable power supply to customers [10], [11].
To this end, the VSI controllers are usually designed to have two main sections, namely ''power sharing'' and ''inverter output'' sub controllers. The former tries to properly share active and reactive power among the sources, while the latter is responsible for regulating output voltages and currents of the inverters [12]. Furthermore, given there are numerous controllers in an MG, they should be coordinated with each other using one of the two approaches, namely, centralized control strategy and local control strategy [13], [14].
In the centralized control strategy, a central controller is responsible for collecting, validating and analysing all the necessary information of the network in order to issue appropriate commands to the objective sources. Despite the high accuracy, this relatively expensive approach lacks reliability due to its dependency on the wide band communication infrastructure as well as the central controller.
The local control strategy usually utilizes common droop-based characteristics for each DG unit to control the frequency/voltage signal without aid from the central controller [15], [16]. In this control strategy, the designers mainly focus on keeping stable power supply to the customers while trying to optimize the simplicity, reliability, cost and expandability of the structure. Although the local strategies may have negligible delays or inaccuracies, they have a lower cost and higher reliability [17], [18].
At this point, it is also noteworthy to say that although VSI controllers increase system flexibility to reach the steady-state requirements of MG, the low-inertia nature of VSIs will decrease system stability against probable disturbances [19], [20]. In other words, the VSI has a number of low frequency and lightly damped modes that may excite the unstable dynamics. This is why a significant portion of the literature, has focused on the stability problems of VSI-based MGs [21]- [24]. For instance, in [23], in order to increase the inertia and the control performance of MG, a cascaded control strategy has been proposed for a voltage source grid-supporting inverter. In this reference, a virtual synchronous generator based control strategy and a sliding mode control based vector control strategy are applied in the droop based control loop and in the inner control loop, respectively. In [24], in order to improve the stability, a Pulse Width Modulation based Sliding Mode Control (PWM-SMC) controller for three-phase grid connected inverter with LCL filter is proposed.
As the lightly damped modes are usually dependent on the slow procedure of the droop mechanism [25], any droop coefficient variation should be consciously and carefully performed. Droop coefficient correction is a common method in solving the inherent power sharing error of droop-based controllers [26], [27]. Although changing the droop coefficients will solve the available steady-state errors in the frequency and voltage of MGs, it may threaten the network stability. On the other hand, as some economic considerations make the MGs operate within a narrow margin of instability, it is necessary to exactly investigate the impacts of the parameter variation on the stability of an islanded inverter-based MG. Small signal analysis is widely used to study the stability of the MGs [21], [22], [25], [27]. However, by changing the parameters, the operating point of the system changes. Therefore, the operating point is deviates from the predefined linearized operating point [28]. The bifurcation theory is applied to predict the parameter stability margin of nonlinear system along with the change of parameters. Using the bifurcation theory, it can be shown that different bifurcations may also lead to the flicker phenomenon in addition to instability.
Various strategies have been proposed in the literature to solve the problems caused by utilization of the conventional droop technique in the islanded microgrids. Virtual impedance is an efficient solution for reduction of power coupling, improvement of load sharing among DG units, adjustment of DG output impedance, improvement of microgrid voltage regulation, improvement of microgrid stability [29]- [31]. The virtual impedance-based control methods focus on reduction of the feeder impedance of sources. In this systems, the voltage reference of the voltage-sourced inverter (VSI) is reduced proportionally to its output current [32], the VSI control system simulates the behavior of an actual impedance.
In [17], [30], control strategies based on adaptive virtual impedance control are proposed. In [30], an adaptive virtual impedance control scheme is proposed for overcome the unbalanced and harmonic power sharing in islanded microgrids. In [17], an adaptive complex virtual impedance method is proposed for ideal reactive power sharing in the islanded microgrids. In [31], an improved droop control method is applied for power decoupling and accurate reactive power sharing. The method is based on virtual power source and composite virtual negative impedance for low-voltage microgrid.
In [33], the resistive virtual impedance is proposed to correct the reactive power sharing caused by the unbalanced loads. In this method, the parameters of the virtual impedance are optimized by the genetic algorithm. An adaptive virtual impedance-based method has been presented in [34] to compensate the voltage drop on the DG feeders. Using the communication link, the reactive power of each DG is sent to the control center for regulation of virtual impedance. Another virtual impedance-based techniques work based on the analysis of small signal stability. In [35], [36], the optimal value of virtual impedance is determined to improve the stability indexes and the load sharing.
In this paper, first, the bifurcation theory is utilized to investigate the probable impacts of droop coefficient variations on the healthy operation of the islanded MG. To this end, a comprehensive dynamic model of inverter-based microgrid is derived. Next, using the virtual impedance theory, a control methodology is presented to guarantee both stability and power quality of the network against common droop correction processes.
The main contributions of this paper are the following: 1) Using continuous variation of power controller's parameters as well as eigenvalue monitoring, different bifurcations are investigated. 2) Both stability and power quality aspects are simultaneously monitored.
3) The impact of common droop correction methods on the stability and power quality of the network is considered. 4) The sensitivity of flicker to each control parameter is explored. 5) The virtual impedance theory is successfully utilized in solving stability and power quality problems.
The rest of the paper is organized as follows: Section 2 proposes a comprehensive model of inverter-based MG. In Section 3, eigenvalues of the model are analyzed to find out how controllers will affect them. Section 4 studies how droop coefficient correction will lead to bifurcation and consequently flicker. Finally, section 5 presents the virtual impedance-based control method to avoid any stability and power quality problems.

II. VSI-BASED MG MODELLING
This section presents a d-q reference frame dynamic model of the islanded MG that is shown in Fig. 1a. The islanded MG mainly consists of three parts, i.e., inverter-based DG units, lines and loads.
In a general MG system, the DG units contain a DC voltage source, an inverter bridge, an LC filter, a coupling impedance and control loops (including the power, voltage, and current controllers), as shown in Fig. 1b. The inverter bridge uses a switching mechanism to convert the DC voltage source to a semi-sinusoidal AC output power. Accordingly, the designers usually consider an LC filter in parallel to high switching frequency to reduce the number of harmonics. Moreover, the control loops adjust the output voltage and frequency of the inverter bridge. Consequently, it is assumed that the whole VSI is similar to an ideal balanced three-phase source from the network's viewpoint. In addition, other assumptions that are considered in the model are listed as below: 1) Nonlinear dynamic of each DG unit is presented in its own d-q reference frame. 2) The reference frame of each DG unit is rotates at the angular frequency of ω.
3) To obtain the state-space model of MG, the reference frame of a pre-selected DG unit is assumed as the common reference frame with a frequency of ω com . According to these assumptions, the reference frame angle of each DG with respect to the common reference frame is calculated as shown in (1): As shown in Fig. 2, the power controller is mainly based on frequency and voltage droop curves. To achieve appropriate power quality and an adequate time delay between the power control loop and voltage and current control loops responses, it is necessary to use average signals of instantaneous active and reactive power [37], [38]. Accordingly, as is illustrated in Fig. 2, the fundamental components of p and q signals, that are passed through low-pass filters with the cut-off frequency of ω c , are calculated through the following equations.
wherev odq andî odq are the output voltage and current in the d-q reference frame, respectively. The differential equations corresponding to (2) and (3) can be written asṖ As mentioned earlier, the DG shown in Fig. 1b uses two droop-based controllers for frequency and voltage in power controller that can be obtained as where, v od_ref and v oq_ref are the reference values for the voltage controller and ω is the angular frequency for the inverter bridge. According to the above equations, the output voltage magnitude reference is aligned to the d-axis of the inverter reference frame while the q-axis reference is set to zero.

B. VOLTAGE CONTROLLER
The voltage controller enables MGs to track the input reference voltage. As shown in Fig. 3a, a standard PI controller, G Vdq (s) = k p_vdq +(k i_vdq /s), compares the measured voltage at the output of DG (v odq ) with the reference voltage given by the power controller (v odq_ref ) to regulate the final output voltage.
Compensating the cross coupling terms, known as ω b C fvod and ω b C fvoq , this controller decouples voltage control loops with respect to d and q axes. Moreover, in order to decrease the impact of load fluctuation on the output voltage, this controller also utilizes a feedforward term with a constant gain of K [27]. Accordingly, the output current of the controller can be obtained as where ω b denotes the nominal angular frequency and C f is the output filter capacitance.

C. CURRENT CONTROLLER
A current controller is used to improve the power quality, internal stability and capacity of MGs. The block diagram of a sample current controller is shown in Fig. 3b. This controller uses a similar PI controller, defined as G Idq (s) = k p_idq + (k i_idq /s), to regulate the output current. To this end, the PI controller compares the measured current of inductor (î ldq ) with the reference current given by the voltage controller (i odq_ref ). Similarly, the current controller also decouples current control loops with respect to d and q axes by compensation of the cross-coupling terms, known as ω b L fîld and ω b L fîlq . Finally, the output of the current controller is used as the input voltage signal of the inverter. According to the block diagram, the output voltage equation can be obtained as: where L f is the inductance of the output filter.

D. OUTPUT LC FILTER AND COUPLING IMPEDANCE
As shown in Fig. 1, LC filters are usually utilized in DGs to improve power quality and the voltage profile. This filter is usually connected to the network via a coupling impedance, assuming v idq = v dq_ref . The voltage and current equations of the output LC filter as well as the current equation of the coupling impedance can be calculated as follows: L c dî odq dt = −r cîodq − jωL cîodq +v odq − v bdq (12) where r c and L c are the resistance and inductance of the coupling impedance, respectively.

E. STATE-SPACE EQUATIONS FOR LINES AND LODES
The RL model is used for representing lines and loads. Accordingly, the current passing through the line, connecting bus i to bus j, and the local current of each bus is calculated as follows:

III. EIGENVALUE ANALYSIS
Equations that are presented in Section 2 illustrate the nonlinear dynamics of the inverter-based DGs, lines and loads of an MG in the d-q reference frame, i.e. these equations have the ability to show the system dynamic model in the state-space equation form. Accordingly, the state space equations of the sources, the loads and the lines can be summarized as below: where ables, vector of inputs and known disturbances, respectively. Furthermore, x 4i , x 5i , x 6i and x 7i are the variables of the ith DG, which are used in the integral terms of (8) and (9). Analysing the small signal stability of an islanded inverterbased MG, it is necessary to linearize the state-space equations of (15) around an equilibrium operating point as shown in (16).
where, A sys represents the state matrix of the system. Accordingly, it is possible to form the respective system state matrix for the MG with any number of inverter-based sources, loads and lines. In this paper, we use the sample network which includes four inverter-based DGs to demonstrate the impact of droop coefficients on the stability margin of MGs. In this network, nominal voltage and frequency are assumed 380 V and 50 Hz, respectively. Fig. 4 shows the single line diagram of the network under study.
As shown in Fig. 4, four inverter-based DG units are connected to this network in which the first two (DG 1 and 2) are 45 kVA, while the others (DG 3 and 4) are 34 kVA. In this network, three lines connect the buses together and two local loads are supplied from bus 1 and 3. It is assumed that the loads are linear and constant impedance. Detailed information about the DG units and their controllers is presented in Table 1. Fig. 5 illustrates the locus of dominant eigenvalues of the network under study when the frequency droop coefficient of the first DG units (m p1 ) has changed from 9.4 × 10 −5 to 30.8320 × 10 −5 with a 0.1 step size. As shown in Fig. 5a,  the dominant eigenvalue λ 1,2 is equal to −6.5766 ± j54.6932 when m p1 = 9.4 × 10 −5 . Increasing the m p1 , this mode first gets away from the jω axis, next, for m p1 larger than 17.860 ×10 −5 , it will gradually get close to the imaginary axis. Finally, when the m p1 becomes larger than 22.560 × 10 −5 , starting from a new point (λ 1,2 = −3.1686±j58.5162), this mode will cross the jω axis. Fig. 5b shows the eigenvalue locus to illustrate how the dominant eigenvalues, λ 3,4 , react as m p1 varies. As shown, when m p1 is equal to 9.4 × 10 −5 , the dominant eigenvalues meet −8.6330 ± j46.2153. Increasing m p1 , while it is still smaller than 21.62×10 −5 , will make the mode move towards the jω axis. For larger m p1 s, the eigenvalues move to a farther point as −6.4858 ± j55.1246. Reaching this point, m p1 will not have any more impact on this mode. Fig. 6 shows the locus of the system dominant eigenvalues (λ 1,2 and λ 3,4 ), when the voltage droop coefficient (n q1 ) is changed from 1.3 × 10 −3 to 4.16 × 10 −3 with a 0.1 step size. As shown in this figure, when n q1 increases, the damping and frequency of λ 1,2 remain unchanged in −6.5766 ± j54.6932. This is while, λ 3,4 gradually move away from the jω axis. Thus, selecting a larger n q1 will increase system stability. However, the operator should be aware of voltage regulation; this fact may lead to small, weak or even unacceptable output voltages in MG.

IV. MG STABILITY AND POWER QUALITY ANALYSIS USING BIFURCATION THEORY
It should be noted that the system is stable and works within its own nominal operational range. However, for a microgrid VOLUME 9, 2021 in operation, determining the critical value of the parameters can be used for predicting the stability margin due to its current state, which is useful for the online scheduling of control parameters.
In order to enhance the accuracy of power sharing, increasing the slope of the voltage and frequency droop characteristics is performed according to the permissible voltage and frequency range. Typically, the permissible range is obtained through the small signal analysis. However, the correct selection of these coefficients within the allowed range is carried out through trial and error due to the inherent compromise between voltage regulation and reactive power sharing as well as frequency regulation and active power sharing. Therefore, there is the possibility of bifurcation occurrence as a consequence of changing the obtained parameters within the permissible range.
In general, a bifurcation is a change of the topological type of the system when the parameters pass through a bifurcation (critical) value. The continuation method is also one of the aspects of bifurcation theory being applied to the stability analysis of a dynamic system with respect to a specific parameter variation [39]. In this method, for a predefined operating (equilibrium) point of a dynamic system, either continuous time or discrete time, periodic or non-periodic, some parameters of the systems are continuously being changed to monitor the following variables in the new operating point:

1) Eigenvalues (for non-periodic continuous time systems) 2) Floquet Coefficients (for periodic continuous time systems and their Poincare mappings) 3) Multipliers (for discrete time systems)
Depending on the location of the above variables in the complex plane, different bifurcations (Saddle-node, transcortical, Pitchfork, Period-doubling, Hopf, Neimark-Sacker) may occur leading the system towards the two following inappropriate phenomena: Instability: As shown in [11], [40], Saddle-node, Hopf and Period-Doubling bifurcations will gradually lead the system towards chaos or voltage instability.
Flicker: Sometimes, a bifurcation may not lead the system towards instability; however, the customers will experience a bother able flicker.
In the continuation method, the first value of the parameter, in which one or two eigenvalues reach the jω axis is called as bifurcation parameter. Fig. 7a shows the steady state waveform of the voltage when m p1 is equal to 9.4 × 10 −5 and consequently the dominant eigenvalues (λ 1,2 ) move to a new point as (−6.5766 ± j54.6932). As shown in this figure, the sinusoidal waveform has a 372.6 V magnitude and a 313.68 rad/s angular frequency. The direct bold line in this figure is the upper voltage envelope.
Given in this case study, only the droop characteristics are responsible for maintaining voltage and frequency within their respective allowable ranges. As illustrated, the voltage and frequency are set near their nominal values. Fig. 7b shows the system phase portrait in a neighborhood of the equilibrium point. As shown in this figure, this equilibrium point is linearly stable and the damping is sufficiently large. Therefore, the phase portrait will quickly move towards the equilibrium point, which can be named as stable focus. Furthermore, Fig. 7c illustrates the frequency spectrum of the voltage waveform for the corresponding m p1 , in which the 50 Hz frequency has a 100 percent amplitude. To obtain this spectrum, the FFT method with a 15 cycle period is selected according to the IEC standard [41]. It is noteworthy that the other nonzero amplitudes are due to the numerical errors in the FFT operation, known as spectral leakage.
As mentioned in the previous section, m p1 specifies the closeness of λ 1,2 to the jω axis. For instance, when m p1 = 24.44 × 10 −5 , these eigenvalues are equal to −2.3936 ± j60.1084. For this m p1 value, the phase portrait, frequency spectrum and voltage waveform are shown in Fig. 8. As shown in Fig. 8a, when m p1 is equal to this value, the bold voltage envelope fluctuations are gradually damped within 0.6 seconds. According to what is shown in Fig. 8b, although the system is still linearly stable, the phase portrait moves toward the equilibrium point slower than in Fig. 7b.
As mentioned earlier, the rapid fluctuations in the rms voltage of a power supply lead to the flicker phenomena.
Generally, the following equation shows a sinusoidal waveform with an angular frequency of ω 0 which is modulated by another waveform with the frequency ω m [42], [43].
Expanding equation (17), the following relation is obtained: This equation represents a modulated signal containing three frequencies namely ω 0 , ω 0 + ω m and ω 0 − ω m . Fig. 8c shows the voltage waveform spectrum for a 15 cycle sampling frequency and within a 0.1 to 0.4 s period. As shown, the voltage waveform consists of a fundamental frequency (50 Hz) with a 370 V amplitude, a subharmonic of 40 Hz with 4.2 % amplitude and an inter-harmonic of 60 Hz with a 3.7 % amplitude. According to different standards, e.g. IEEE Std. 1159-2009, IEC 61000-4-15, an annoying flicker has a frequency range of [0. 1 30] Hz and a voltage range of [0. 1 10] percent. Given the system frequency is 50 Hz, the two 40 and 60 Hz harmonics make a 10 Hz flicker in the studied network.

A. HOPF BIFURCATION
To clearly highlight the flicker problem in the network under study, m p1 is increased to the critical value as m * p1 = 29.945016 × 10 −5 according to what is shown in Fig. 9a. In such a condition, a pair of eigenvalues is placed on the jω axis (λ 1,2 = ±j64.47796) while the others are still in the left half of the plane. As a result, for the above coefficient, the system is exposed to Hopf bifurcation.
For every m p1 smaller than m * p1 , the equilibrium point of the system is linearly stable making the phase portraits maintain stable focus in the equilibrium neighborhood.
When m p1 is equal to m * p1 , the equilibrium point becomes a nonlinear stable point. Thus, the convergence rate is no longer exponential. For m p1 larger than m * p1 , the equilibrium point vanishes to be replaced by a unique and stable limit cycle. In this case, the system is still stable while there are oscillations. This bifurcation, named a supercritical Hopf bifurcation, occurs when all the neighboring trajectories converge toward the limit cycle. Fig. 9b shows the phase portrait of the system when the Hopf bifurcation happens. Fig. 9c also shows the voltage waveform frequency spectrum for the corresponding m * p1 in which the fundamental frequency is 50 Hz with a 368.1 V amplitude. The subharmonic frequency is 40 Hz with a 8.1 % amplitude and the VOLUME 9, 2021 inter-harmonic frequency is 60 Hz with a 7.9 % amplitude. These results illustrate an annoying 10 Hz flicker on the voltage waveform of a sample MG. Hence, the MGs may easily experience Hopf bifurcation and consequently power quality problems due to a probable parameter tuning process.

B. PERIOD-DOUBLING BIFURCATION
As mentioned in the previous subsection, increasing droop coefficient (m p1 ) lead to a stable limit cycle. Accordingly, this subsection increases this coefficient further to investigate the impact on the system. If m p1 exceeds m * p1 , a flip (or period-doubling) bifurcation will take place at m p1 = 30.7192 × 10 −5 .
In this condition, voltage waveform oscillations increase such that the fundamental decreases to 364.3 V (See Fig. 10a). Fig. 10b also shows the system phase portrait.
According to the increasing oscillations in the voltage waveform, it is obvious that the flicker is significantly increased for the corresponding droop coefficient.
The voltage frequency spectrum of Fig. 10c is proof of the power quality problems of the MG; the 40 Hz sub-harmonic and 60 Hz inter-harmonic have increased to 11.5%, while two new sub and inter harmonics of 30 Hz and 70 Hz are added to the spectrum with a 2.2% amplitude.
Should the m p1 increase more and more, the flicker problem will worsen, i.e. increasing m p1 to 31.0952 × 10 −5 build a new solution (period-two cycle) which is roughly twice the fundamental period (see Fig. 11b). At m p1 = 31.4712 × 10 −5 , period 2 bifurcates to a period 4 solution having a four-time period (see Fig. 12b). As a detailed description, Fig. 11a shows the voltage waveform and the respective envelope when m p1 is equal to 31.0952 × 10 −5 . In this case, the fundamental is equal to 360.7 V. According to Fig. 11c, the amplitude of the 10 Hz flicker is noticeably increased to the value of 14.5%. Power quality problems are more obvious in Fig. 12a where the voltage fundamental decreases to 356.9 V. As shown in Fig. 12c, the 10 Hz flicker has a magnitude of 18%, while there is also a 20 Hz flicker with a 4% magnitude.

V. VIRTUAL IMPEDANCE BASED CONTROLSTRATEGY
Using droop based control strategies in low-voltage MGs leads to the three following challenges [36], [44]: 1) In case of low X/R ratios, the active and reactive powers are mutually coupled.
2) The reactive power is not properly shared among DGs due to different line voltage drops.
3) The DGs dynamic modes are most likely to be poorly damped. Given these challenges are mainly due to the line impedances connecting the DGs to the MG, implementing the virtual impedance based techniques can effectively solve the issues. To the best of our knowledge, up to now, researchers have mainly focused on the virtual impedance techniques to properly share reactive power among DGs [45], [46] or definitely regulate the voltage profile [47], [48]. This is while, it seems these kinds of techniques have also the ability to solve the aforementioned dynamic and power quality problems.

A. IMPLEMENTING VIRTUAL IMPEDANCE TECHNIQUE
The investigated virtual impedance concept for an inverterbased DG unit is shown in Fig. 13. The virtual impedance, known as Z v = R v +jX v [22], is implemented in the following equation of droop control strategy: where v * odq_ref denotes the reference value forv odq in the voltage controller. As mentioned in Section 2, according to (7), v odq_ref is the output voltage of droop controller which is equal to v * odq_ref , in the absence of virtual impedance. It is also noteworthy that according to (19), implementing the virtual impedance will add new terms to the scalardifferential equation (8). Thus, the virtual-impedance based system state matrix will be different in comparison with (17).

B. VIRTUAL IMPEDANCE CONTROL IMPACT ON BIFURCATION
In this subsection, the Z v = 0.1 + j0.6 is utilized as a virtual impedance in the network under study. As shown in Fig. 14, the virtual impedance has significantly enhanced the system dominant eigenvalues in the complex plane. In other VOLUME 9, 2021 words, the virtual impedance modes are increasingly more damped than simple-droop modes of the system. Therefore, the system stability margin can be significantly improved by an appropriate virtual impedance technique. Fig. 15 illustrates how the dominant modes move within the complex plane when the X v parameter of the system changes from 0.05 to 0.8 ohm while other parameters are fixed with the values of Table 1. As shown in Fig. 15, when the X v increases, the dominant modes become more damped due to the increasing virtual distance between the corresponding DG and the others. Moreover, even for the smallest X v (X v = 0.05), the dominant modes are more stable than the modes related to the simple droop control. Thus, it is possible to say that increasing X v is an appropriate strategy for increasing system stability and delaying the occurrence of bifurcation. Fig. 16a shows the dominant eigenvalues locus for the studied MG, when the m p1 is changed in a predefined range while the X v is fixed at 0.05 ohm. As shown, in such a   condition, the hopf bifurcation takes place at m p1 = 5.2875× 10 −4 which is much larger than the corresponding m p1 without virtual impedance (m * p1 = 29.945016 × 10 −5 ). Similarly, Fig. 16b shows the dominant eigenvalues locus when the X v is equal to 0.6 and m p1 changes from 9.4 × 10 −5 to 137.24 × 10 −5 (a non-practical point only considered for simulation analysis). As shown in this figure, X v has a significant impact on the effect of m p1 in the variation of the dominant modes. For this value of X v , the dominant modes change less than in Fig. 16a when m p1 increases. Accordingly, knowing the variation range of m p1 , it is possible to select the best value for X v to avoid any instability or power quality problems.

VI. CONCLUSION
The MGs usually tune the droop coefficients to properly share power among DG units. This coefficient correction process may lead the VSI-dominated MGs towards a problematic operating point being threatened by instability or power quality problems. Using the bifurcation theory, we have shown the system may be exposed to flicker in advance of instability. Accordingly, it is necessary to tune droop control schemes more consciously to avoid any problem. To this end, first, we analyzed the impact of coefficient variation on stability and power quality using the bifurcation theory. Additionally, we developed a virtual impedance based technique to make the dominant eigenvalues more stable during coefficient correction process. Finally, our proposed technique was successfully applied to a practical network. The simulation studies, first, prove the importance of this less-considered, probable and dangerous problem in the MGs, and next, illustrate the efficacy of our proposed virtual impedance control strategy as an affordable solution.