Virtual Impedances Optimization to Enhance Microgrid Small-Signal Stability and Reactive Power Sharing

Microgrid instability poses critical issues to the power delivery following a load change or a tripping event. In island operating mode lack of grid intensifies this challenge. This study aims at controlling several converter-based distributed generations (DG) sharing the power in an island microgrid (MG). At first, the microgrid model including virtual impedances and phase-locked loop (PLL) is introduced. Afterwards a novel small-signal stability analysis for island microgrids is proposed. Finally, an optimization algorithm based on particle swarm optimization (PSO) is proposed to design the virtual impedances. The optimization algorithm analyzes all possible operating points and aims at maximizing the microgrid stability index while keeping the reactive power mismatches at minimum level. The fractional objective function facilitates reaching at these objectives simultaneously. The proposed optimization algorithm is implemented in two separate case-studies and the corresponding virtual impedances are drawn in any microgrid. On the other hand, The voltage drops are checked as a condition in the optimization process. The results drawn from two separate case-studies verify that the proposed algorithm effectively maximizes the microgrid stability index and minimizes the reactive power mismatches.


I. INTRODUCTION
The worldwide concern over greenhouse gas emissions and future energy resources has motivated the energy sectors towards renewable energy sources (RES).Microgrids (MGs) are critical building blocks of future power systems which facilitate the integration of RES into traditional power systems.However, MGs pose imperative challenges to power systems, such as the complexity of control and stability, lack of intrinsic inertia, a higher number of generating units and intermittency of prime movers.Similar to control fundamentals of traditional power systems including multiple generators [1], droop control has been deployed in MGs to imitate the power sharing regime of the machine-based power systems [2].The voltage and frequency droop equations are The associate editor coordinating the review of this manuscript and approving it for publication was Min Wang .employed in separate control loops.However, active and reactive power sharing coupling in complex MGs was reported causing non-accurate power sharing [3].
The idea to deploy virtual impedance to correct the non-ideal reactive power sharing; known as reactive mismatch, was presented in [4], [5].Virtual impedances are chosen in specific intervals considering the MG small-signal stability analysis, voltage limits, reactive power sharing and demanded damping status [6].The virtual admittance was also adopted in [7] to be applied to parallel current-controlled DGs (either in grid-connected or island MGs) to compensate harmonics currents and to reduce transmission losses.The adaptive virtual impedance proposed in [8] tends to amend the active and reactive power sharing in a meshed island MG.
The novel droop control proposed in [9] deploys virtual impedances and resolves the coupling effect and decreases the power circulation among DGs and stabilizes the MG.The VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License.For more information, see https://creativecommons.org/licenses/by/4.0/complex virtual impedance was suggested in [10] to stabilize the voltage and to have a correct power sharing in island MG.The robust virtual impedance in [11] was introduced to improve the MG stability, mitigating voltage distortions and enhancing the post-fault behavior of converters.However, MG instability following a load change or tripping of a DG in island MGs is a current challenge [12].MG stability can be evaluated by small-signal stability analysis around an operating point; several parameters such as controller coefficients and phase locked loop (PLL) coefficients play major roles in the MG stability [13], [14].The lower order small-signal model in [15] and [16] were introduced for both isolated or grid-tied MGs which facilitates the small-signal stability analysis.The optimal ranges for MG parameters based on MG small-signal modeling were determined in [17] to enhance the dynamic behavior of converters.The model-based voltage and current controllers were compared to conventional PI controllers regarding MG small-signal stability in [18].
Moreover, the small-signal model of a two-layer MG control structure was presented in [19] which contains two main layers (MG layer and cluster layer) and MG parameters were designed accordingly.The real case microgrid of Dongao island was under study in [20], where MG frequency control method includes three operational zones, stable, precautionary and emergency; The MG stability analysis was applied within second and third zones.MG reconfiguration practices can enhance the stability margin considerably as reported in [21].
The presence of both constant power loads and induction motor loads in MG affects the small-signal stability of the MG as reported in [22].The distributed secondary control in [23] was proposed to present an enhanced dynamic performance based on the MG small-signal stability.Moreover, the advantages of distributed control was scrutinized in [24] using a theoretical basic framework for small-signal stability.
The power sharing effects on small-signal stability of a hybrid MG including low-inertia DGs and diesel generators was analyzed in [25].By cascading lead compensators, a technique was proposed in [26] to enhance MG small-signal stability.The Popov's Absolute Stability Criterion was applied to analyze the MG stability conditions while there are constant power loads installed in the MG [27].The small-signal characteristic equation of MG can be drawn to measure the MG low-frequency stability as in [28] which applies Padé approximation and dynamic phasor model.It was theoretically approved that the common DG connection topology is harmful to small-signal stability of MG [29].
The distributed consensus methodology in [30] was proposed to adaptively share the harmonic load among DGs.The differential algebraic-equation model of MG was used to draw the stability region of a MG in [32].The bifurcation theory analysis was applied in [33] to analyze the parameter stability region of a MG including different types of loads.The MG stability analysis in harmonic conditions was studied in [34] which uses the concept of dynamic phasor (DP) to explain the harmonic components of an ac signal as dc variables.
This study firstly proposes a thorough small-signal dynamical model for the island MGs including PLL and virtual impedances.Then the small-signal stability analysis is performed based on the dynamic model of MG.Afterwards, a novel optimization algorithm to design the virtual impedance for converters in the MG is introduced that maximizes the MG stability index and minimizes the reactive power mismatches.The proposed optimization algorithm analysis the microgrid stability in all possible operating points and the fractional form of the proposed objective function facilitates reaching the objectives simultaneously.This study: • Applies an enhanced dynamical model of MG which consists of PLL dynamics and virtual impedances to enhance the studies in [13], [14] and [18].Both PLL and virtual impedances models are critical elements in MG stability analysis.
• Proposes a novel fractional objective function to design virtual impedances to minimize reactive power mismatches between converters and to enhance the critical eigenvalue of the microgrid.The fractional objective function removes the dual problem of setting the weighting coefficients in a multi-objective objective function as proposed in literature.
• The microgrid small-signal stability is analyzed in all operating points and the stable operation of microgrid in presence of virtual impedances is assured.
The remainder of this paper is arranged as follows.
Section II introduces MG modelling and small-signal analysis.Section III describes the virtual impedance design algorithm.The simulation results and following discussions are provided in Section IV.The conclusions are drawn in Section V. Appendix and References appear subsequently.

II. MICROGRID MODELING AND SMALL-SIGNAL STABILITY ANALYSIS
Distributed generation is typically connected to a MG via a power converter, which is either a voltage source or a current source converter.In island MG applications, voltage source converters often connect DGs to MG and they are called to be in grid-forming mode.It means that they assign the voltage and frequency set points and form an island MG.A MG is stable if subsequent to a disturbance, all state variables reach at steady-state values which satisfy operational constraints [31].While there are several DGs injecting power to the MG, they must be controlled on a common reference frame.The reference frame can either be synchronized to one of available DGs or an external reference based on a common time [36].The latter requires no communication links among DGs.For any DG in a MG, the effective angle is calculated as the angle between the reference frame of the voltage of the DG and the common reference frame as depicted in Fig. 1.It can be concluded that if one of the converters is chosen as the  common reference frame, the effective angle for itself is zero; but the others have deviation angles with respect to it [14].The proposed control block diagram is illustrated in Fig. 2. A PLL device is required to measure the frequency or synchronize a converter to the MG.
As the MG is controlled on a common reference frame, transformation from local to common reference frame and vice versa seems necessary.All over this research work the local dq values are annotated by dq indices while the common frame values are annotated by DQ values.The bus voltages are main common variables (v bD , v bQ ), which are used in lines and loads dynamic equations.But, while writing dynamic equations for any converter, local bus voltages (v bd , v bq ) are required which should be calculated using Equation (1).In case of currents, when the output currents are calculated for any converter local currents (i od , i oq ) should be translated to common reference frame yielding delivered currents (i oD , i oQ ) to the MG.The angle δ com for any converter is the angle between its own local d-axis and the common D-axis in radian degree, as it is seen in Fig. 1.
Hereafter, the dynamic equations of all blocks applied in the proposed control block diagram (Fig. 2) are described.

A. POWER CALCULATOR
The instantaneous active and reactive powers (p, q) are calculated based on output dq components of voltage and currents of the converter.The instantaneous active and reactive powers pass a low-pass filter (LPF) afterwards.The cut-off frequency for the filter is ω c .Two dynamic equations are yielded afterwards, which will be used in the MG small-signal stability analysis.The power calculator and LPF inputs and output are depicted in Fig. 2.
The virtual impedance is considered in voltage loops which causes a voltage drop as seen in Equation ( 7).Where R v , X v and V vir are virtual resistance, virtual reactance and virtual voltage drop, respectively.The d-axis voltage component is fixed at zero according to the control strategy.
The positive values of virtual impedances will be determined in an optimization algorithm in the next section.
The small perturbation around any operating-point will affect the voltage drops on virtual impedances as seen in Equation ( 8).
C. DROOP CONTROL The droop control in Fig. 2 includes two traditional droop control equations as follows.Where ω * and V * oq are frequency set-point and q-axis voltage set-point respectively.The d-axis voltage set-point is zero in the proposed control method.m and n are active power droop coefficient and reactive power droop coefficient respectively.V vir is the virtual voltage drop on virtual impedance calculated by Equation (7).While the active and reactive power set-points are zero, the voltages and frequencies deviate from nominal values (ω n , V oq n ).However, if a certain active and reactive powers (P 0 , Q 0 ) are assigned for the converter, these values should be considered in droop equations.

D. VOLTAGE CONTROLLER
As the voltage controller Equations ( 13), ( 14) were exactly drawn from research in [13], the corresponding equations are mentioned here respectively for the sake of simplicity.However, any converter has its specific ω * during the time.The proportional (P) and integral (I) coefficients of the proportional-integral (PI) controller are k pv and k iv respectively.Two auxiliary variables ϕ d and ϕ q are defined and used to simplify the dynamical modeling of system.ω PLL is the MG frequency measured by PLL (rad/sec).The ω * and v * oq are frequency and voltage set-points of MG according to droop equations.Two first order derivations are annotated by φd and φq which will be used in small-signal stability analysis of the MG.The outputs of voltage controller block in Fig. 2 are dq current commands (I * ldq ).

E. CURRENT CONTROLLER
The current controller generates voltage commands or references for the converter in dq reference frame as explained by Equations ( 15), (16).It is notable that k pc and k ic are corresponding PI controller coefficients in the current controller.ω n and L f are nominal frequency of MG and the filter inductance respectively.Two auxiliary variables γ d and γ q are defined to simplify the dynamical modeling of the MG.The outputs of the current controller are two dq voltage commands (v * id , v * iq ).It is assumed that the generated voltage components of converter are identical to these voltage commands.

F. PLL MODEL
A PLL model proposed in [13] in dq reference frame is utilized which forces the d-axis voltage towards zero.It means that when the PLL tracks the phase angle, in steady state operating condition the d-axis component of voltage is zero and q-axis component of voltage is aligned with the Q axis.
The following PLL model uses a low-pass filter with the cut-off frequency of ω c,PLL and a proportional-integral controller with coefficients k p,PLL and k i,PLL .The parameter (f n ) in Equation ( 22) is the nominal frequency of the system (Hz).Equations ( 19), ( 21) and ( 23) are used in small-signal stability analysis of MG.

G. LOAD MODEL IN THE PROPOSED MODEL
An impedance load model is considered in this study.Applying traditional Kirchhoff's circuit laws and DQ frame concept, one can extract following equations on common reference frame: where R load and L load are the resistance and the inductance of series resistive-inductive (RL) load and ω PLL is the frequency of the MG (rad/sec) measured by PLL.The load currents are functions of themselves and the bus voltages.A common practice is to use virtual resistor method to define bus voltages according to state variables.The advantage of this method is that when reaches to extracting state space matrices ''B'' is zero and on the other hand the load dynamic effect on bus voltages is taken into account.Using virtual resistor idea, the bus voltages which are inputs to the system are formulated as functions of state variables as follows: where r N is the virtual resistor value which is assumed to be greater than all load resistances in the MG.The ± sign suggest that when the line current is flowing out of the bus, the minus sign is used in Equations ( 26), ( 27) and if the line current is flowing into a bus, the plus sign is used in Equations (26), (27).The bus voltages in common DQ reference frame (v bDQ ) are defined as a function of output currents (i oDQ ), load currents (i load DQ ) and line current components (i line DQ ).

H. LINE MODEL IN THE PROPOSED MODEL
A typical resistive-inductive line is considered between two adjacent buses.It is noteworthy that the common reference frame is located on all MG buses, so corresponding voltages are used while writing dynamic equations for lines.For a line connected between bus i and bus j the dynamic equations can be expressed as: where r line and L line are the resistance and inductance of transmission line in the MG, respectively.The measured frequency of MG is ω PLL in Equations ( 28), (29).Moreover, the direction of line current is assumed arbitrarily, so while calculating bus voltages for two adjacent buses as in Equation (25) and Equation ( 26) the sign of last term for two buses will be positive and negative respectively.

I. MICROGRID SMALL-SIGNAL STABILITY ANALYSIS
Considering previous dynamic equations for the voltage controller, current controller, PLL, load, LPF and interfacing lines, following state variables are considered for a typical two node MG.Afterwards, small-signal equations for all these state variables are written accordingly.The small-signal equations for the MG are nonlinear and they will be linearized around an operating point.These small-signal equations are not brought here because of the page limit but the A matrix which is the final representation of all small-signal equations is written in the Appendix.
The ''A'' matrix for a 2-Bus MG is a square 36 × 36 matrix which will be explained in the Appendix.For a n-Bus MG the small-signal stability analysis and the A can be drawn similarly.As the small-signal equations are non-linear, the linearization around an operating point is necessary.The sign denotes a small change of any state variable which is linearized around the operating point.The linearization process was described in previous researches (e.g. in [10], [13]) as well, but it is not repeated here because of page-limit.
The following operating point is used as a starting point for load-flow analysis, but for the small-signal stability analysis which is done off-line, all possible operating points will be examined which is well-explained in the next section.

III. VIRTUAL IMPEDANCE DESIGN
The proposed design algorithm to design virtual impedances is demonstrated in Fig. 3.It is run off-line and the designed virtual impedances are then installed in converter controllers.
The voltage set-points are assigned at First.Then, a positive interval is defined for the virtual resistance and inductance.It was proposed in [8] to choose the virtual resistance equal to 0.2 of virtual reactance.However, equal intervals are considered to enable optimizing the objective function.The maximum values are calculated based on maximum impdenace mismatches.Particle swarm optimization (PSO) is a powerful stochastic tool for solving optimization problem, the description of its application in power system issues is explained in [35].The higher convergence capability, an enhanced global search, the relative simplicity of parameters tuning and the higher precision of solutions motivated the authors to apply PSO algorithm.The initial population consist of virtual resistances and inductances.The PSO algorithm in any iteration calls the small-signal analysis and the MG load-flow to draw the MG stability index, injected currents and consequently the objective function.
The operational flowchart for the proposed virtual impedance design is depicted in Fig. 3 which can be summarized as following stages: 1) Initialization of the optimization variables which are L v1 , L v2 , . . ., L vn , and R v1 , R v2 , . . ., R vn .It should be noted that the initial values are also chosen inside the stable and permitted interval of any variable.
2) The number of PSO population, and the iterations and the other parameters of the algorithm are assigned and the PSO algorithms is run.
3) The load-flow analysis is run in the desired time interval applying MATLAB which can include load changes in any bus at any moment.Several operating points are drawn from the power flow analysis depending on the events occurred during the time interval.The number of operating points depend on the solver time step.4) For all operating points considering the virtual impedances, the eigenvalues stability analysis is performed to examine the microgrid small-signal stability.
A microgrid stability index which is the absolute value  of real part of worst eigenvalue (the nearest eigenvalue to the vertical axis of the complex plain) among all non-zero eigenvalues of the microgrid.If all eigenvalues have negative real parts and there is only one zero eigenvalue, the MG is asymptotically stable and the stability condition is satisfied.Otherwise the algorithm returns to stage 2. It should be noted that if there are m eigenvalues for the microgrid and the load-flow yield n operating points, m × n stability indexes are drawn.5) If the stability condition is satisfied then the λ and the set of optimization solutions are saved and the algorithm advances to step 6. 6) The simulation is run in a certain time interval, if the time exceeds t f then the load flow is terminated.If the time is still lower than t f the simulation continues to analyze the other operating points.7) The convergence criterion could either be the number of iterations or the magnitude of the objective function.The former is considered in this study.However, the nominator of the objective function is the summation of the reactive power mismatches for all converters multiplied by their reactive power droop coefficients and the denominator is the minimum λ among all operating points and eigenvalues.8) The bus voltages in microgrid must remain inside the acceptable interval; so applying the optimal virtual impedances, if the voltage limits are not served, the algorithm returns to stage 1 and lowers the upper bounds of virtual impedances and virtual inductances in order to reduce the voltage losses.One another possibility is to cover the voltage drops by changing the voltage set-points of converters around 1± 0.05 p.u. to deliver the satisfactory voltage to the loads at buses.9) If the voltage limits are served, the the algorithm generates optimal values for virtual inductances and resistances.The applied PSO algorithm has several parameters which are set as seen in Table 1.After applying the algorithm presented in Fig. 3 the optimal inductances and resistances are calculated and they are applied in the converter control system.

IV. SIMULATION RESULTS
A 2-Bus MG introduced in [13] and a 3-Bus MG introduced in [14] are implemented in MATLAB as two separate case studies to validate the effectiveness of the proposed optimization algorithm.The results for any case study are analyzed and discussed in detail.

A. CASE STUDY 1: 2-BUS MICROGRID
A 2-Bus test microgrid shown in Fig. 4 which was under study in [13] is implemented as the first case study.Table 2 demonstrates full specifications of the controllers and the 2-Bus MG under study [13].The proposed virtual impedance designing algorithm was applied to obtain optimal virtual resistances and inductances in this MG.The corresponding virtual impedances are presented in Table 2. Afterwards, At Time= 2 seconds, R pert1 + jX pert1 are switched in to the MG and is placed in parallel with R Load1 + jX Load1 at bus 1.The local load of bus 2 is R Load2 +jX Load2 during the scenario.Hereafter, different characteristics of the MG following this load change scenario are scrutinized.

1) MICROGRID STABILITY INDEX
The MG stability index is defined as the absolute real part of the worst eigenvalue, which is the most subjected to instability.Installing an virtual impedance can mitigate this  mode and maximize the absolute real part of corresponding stability index.Fig. 5 compares the major eigenvalues of the studied MG in two separate scenarios.The MG stability index in former study [13], denoted as λ; and the corresponding MG stability index applying the proposed control framework; denoted as λ v , are seen in Fig. 5.The MG stability index rises from 1.4 to 3.9 deploying the designed virtual impedance.In other words, applying the optimal virtual impedance in the proposed control framework makes the MG more stable.
The MG eigenvalues are listed in Table .3.The MG eigenvalues are listed in second column and damping ratios are written in third column.As it can be seen in the third column, all the modes are damped and stable.However, some modes  are well-damped (e.g.mode 1,2) and have the greater damping ratios (ζ = 100%).The damping ratio corresponding to the MG stability index is 16.155%.As the MG has one eigenvalue located on zero point and the other eigenvalues have negative real part, the MG is asymptotically stable [1].

2) MICROGRID FREQUENCY ANALYSIS
The MG frequency is measured by installed PLLs on the buses.The applied PLLs use PI controllers in control loops and the parameters of the PLLs are mentioned in Table .1.Following a load change at t= 2 sec, the frequency decreases base on active power-frequency droop characteristic.Fig. 6 presents the frequency of the MG following the load change scenarios, using the novel control strategy and applying the previous method in [13], respectively.The rate of change of frequency (ROCOF) and the minimum point (Nadir) are almost identical applying the former strategy and the novel control method, although, the Nadir value for the novel strategy is a little less than the former method.
3) CONVERTERS OUTPUT POWERS Fig. 7 provides output active and reactive powers injected by converters in two scenarios; adopting the novel control strategy and deploying the control method in [13], respectively.Active power generated by converters applying the proposed control framework (P out ) in Fig. 7 settles faster than the corresponding active power in [13].Moreover, the overshoot is less in the former than in the later.The steady-state active powers are identical applying either former or proposed control strategy.The reactive powers served by the converters adopting the proposed control framework (Q out ) for converters 1 and 2 are almost identical before the load change.After the load change event at Time=2 sec both converters share the reactive power equally (Q 1 = Q 2 = 100 Var).As both converters are identical, this reactive power sharing is perfect.However, while applying the previous control method in [13] the reactive power shares after the load change scenario were not ideal, the injected reactive powers are far different (150 VAr and 50 VAr respectively).

4) CONVERTERS TERMINAL CURRENT COMPONENTS
The terminal currents which are denoted as i ld and i lq are dq components of I l−abc in Fig. 1.Obviously, these currents are pre-filtering values.Terminal currents of converters employing the novel control structure and adopting the previous one in [13] are depicted in Fig. 8. Comparing the d-axis components verifies that applying the novel control strategy, the d-axis currents of two converters are equal after load changing event, while these components were non-similar in the previous method.This fair current sharing avoids the overloading of one converter and under loading of the other.Afterwards, examining the q-axis current components admits that the current fluctuations and settling times while applying the novel control strategy are lower than the previous case adopting the other control strategy.However, the steady state q-axis currents in both scenarios are identical.

5) CONVERTERS OUTPUT CURRENT COMPONENTS
As it can be seen in Fig. 9, the proposed control strategy turns out to be more effective to share current components.d-axis current components of two converters in Fig. 9 are almost equal either before the load change or after it (i.e.0.8 A).However, two converters in the previous study [13] inject different d-axis currents even though two converters are similar (i.e.1.2 & 0.4 A).In case of q-axis output currents, the steady state values are identical, but the current overshoot and settling time while applying the proposed novel control strategy are shorter (Fig. 9, bottom graph)

6) OUTPUT VOLTAGE COMPONENTS OF CONVERTERS
The output voltages components for two different scenarios; applying the novel control and employing the strategy in [13] are demonstrated in Fig. 10.According to this figure, the d-axis components of voltages in both methods are forced to remain zero, a negligible fluctuation is seen at t= 2 sec which is mitigated after half a second.Fig. 10 verifies that installing virtual impedances cause voltage drops on these impedances; Hence, a voltage drop of nearly 1.5 V is seen on the q-axis output voltages while applying the novel strategy.It is worthy to mention that, to deliver an identical output voltage the voltage set-point is set on 89 V.This is also mentioned as one of the optimization steps in Fig. 3.

B. CASE-STUDY 2: 3-BUS MICROGRID
A well known 3-Bus low voltage microgrid (220 V RMS) depicted in Fig. 11 introduced in [14] is chosen as the second case-study because it is the most complicated case with resistive inductive transmission lines among DG units and the lines have diferent parameters.Several researches have been  validated by this MG from 2007 until now.The parameters of the 3-Bus MG are listed in Table .4.First of all, applying the proposed optimization algorithm, the virtual impedances are drawn for three converters in the 3-Bus MG.Table .4provides the virtual impedance values for 3-Bus MG which have been drawn applying the proposed optimization algorithm.The proposed values insure that the MG remain stable and the whole modes are damped completely.Moreover, the reactive power mismatches among converters are minimized.The important point is that while applying the algorithm, the 29 kW load was installed in bus 1 of 3-Bus MG and the other buses do not have local loads.At Time = 2 sec the 29 kW load (R pert2 = 5 , X pert2 = 1000 ) is switched in at bus 1 and different characteristics of the MG are analyzed hereafter.

1) MICROGRID STABILITY INDEX
The absolute minimum real part of non-zero eigenvalues has been considered as the MG stability index (λ) in this study.The proposed optimization algorithm enlarges this value as much as possible.Fig. 12 demonstrates the major eigenvalues of the 3-Bus MG while applying the proposed algorithm and install the virtual impedances of Table.4 in the MG and also the case without this virtual impedances which is the control method in [13].The full eigenvalues of the 3-Bus MG while applying the proposed optimal virtual impedances are demonstrated in Table.5.It is seen that the minimum damping is related to mode 25,26 which is 25.528 % and FIGURE 12. Major eigenvalues of the 3-Bus MG in two scenarios, applying the proposed control method (black squares) and the method in [13] (red diamonds).most of the other damping modes are damped perfectly and have the damping ratio of 100 %.It is seen that the MG stability index before installing optimal virtual impedances is λ = 0.311 and the MG stability index after applying the proposed virtual impedances is λ v = 2.01.The proposed virtual impedances could apparently enhance the MG stability index.

2) MICROGRID FREQUENCY ANALYSIS
The frequency of converters in 3-Bus MG are demonstrated in Fig. 13 for two different cases; one which applies optimal virtual impedances and the other case is exactly using the control method in [13].The load change of 29 kW occurs at bus 1 at Time = 2 sec.Fig. 13 shows that the point of minimum frequency (Nadir) while applying the proposed optimal virtual impedances is 49.84 Hz and on the other hand the point of minimum frequency while applying the control method in [13] is 49.71 Hz which demonstrates the enhanced performance of the proposed control and optimal virtual impedances in MG frequency control.The rate of change of frequency (ROCOF) is also lower while deploying the proposed optimal virtual impedances.In both cases the frequency drop at bus 1 is higher than the other buses because the load change occurs at bus 1.The steady state frequency of MG while applying the proposed method and the method in [13] are 49.87 Hz and 49.85 Hz, respectively which are roughly identical.

3) CONVERTERS OUTPUT POWERS
The powers injected by converters to the MG in two separate cases; applying optimal virtual impedances and using the method in [13] are demonstrated in Fig. 14  top of Fig. 14) in load change scenario has a lower overshoot than the case which uses the method in [13] (red dashed line in the top of Fig. 14).In the steady-state, any converter injects 9.74 kW while applying the proposed optimal virtual impedances while any converter injects 9.74 kW in case the method in [13] is applied.The differences in the steady-state injected powers are negligible.However, the reactive powers injected by three converters while applying the proposed optimal virtual impedances are roughly 285 VAr which is consumed in transmission lines.In case the method in [13] is used the reactive powers for converters 1 to 3 are 5.132 kVAr, −1.292 kVAr, −3.135 kVAr, respectively.The proposed optimal virtual impedances could successfully remove the reactive power exchanges among converters and facilitates the fair reactive power sharing.

4) CONVERTERS OUTPUT CURRENT COMPONENTS
The output current components of three converters are depicted in Fig. 15 while the proposed optimal virtual impedances are installed (i od new) and in case the method of [13] is deployed (i od in [13]).The q-axis current components of converters while deploying the optimal virtual impedances demonstrate a lower overshoot (maximum overshoot = 38.9A) than the case with control method in [13] (maximum overshoot = 50.88A) which is a great advantage of the proposed control method.The d-axis current components while installing the optimal virtual impedances are depicted in the top part of Fig. 15.It is visible that the q-axis output currents while applying the optimal virtual impedances reach at identical values (0.763 A); However the converters 1 to 3 take different q-axis currents 3.349 A, −2.553A and −8.12A while applying the method in [13].The proposed optimal virtual impedances successfully remove the q-component current exchanges among converters.

5) THE OUTPUT VOLTAGE COMPONENTS OF CONVERTERS
The output voltage components of three converters in 3-Bus MG are depicted in Fig. 16.The d-axis voltage component (v od ) either using the proposed method or the method in [13] have a zero steady state value.The q-axis voltage components while installing the optimal virtual impedances start from 400 V and after the 29 kW load change at Time = 2 sec reach 367.9 V, 379 V, 387.4 V for converters 1 to 3, respectively.The voltages after the load change have roughly less than 3 % voltage drop which is acceptable.The point is that the initial voltage was 1.05 p.u. to compensate the voltage drops on virtual impedances.The q-axis voltage components while applying the method in [13] are shown as dashed lines in the bottom part of Fig. 16.The v oq1 , v oq2 , v oq3 while deploying the method in [13] start from 380 V and reach at 374.3 V, 382.3 V,

V. CONCLUSION
A small signal model for the island MG including PLL and virtual impedances is developed in this study.Afterwards, a PSO-based optimization method is introduced which determines the virtual impedances for the converters in the MG based on the MG small-signal stability analysis and reactive power exchanges in a MG.It analyzes the small-signal stability in all different operating points to insure the microgrid stability while applying the virtual impedances.The microgrid stability index is calculated for all operating points yielded from load-flow analysis and the minimum stability index is the critical one which will be enhanced using the designed virtual impedances.The designed virtual impedances are installed in all converters in the MG.While applying the proposed control methodology to converters in the MG, a load changing scenario is enforced to the MG.The MG voltage and frequencies are kept within standard limits, the MG stability index is enhanced and the reactive power mismatches are minimized.Moreover, the transient behavior of current components and active powers are enhanced comparing to another recent research.The results in two separate case-studies approve that an optimal virtual impedance enhances the MG small-signal stability in the presence of PLL, which has been a serious challenge in several previous researches.Moreover, the reactive power mismatches problem which has been a challenge in several recent researches is resolved.In case study 1, the MG stability index (λ) increased from 1.4 to 3.9 when the proposed method is applied rather than another previous control method which is an obvious enhancement.In case study 2, the MG stability index (λ) increased from 0.311 to 2.01 when applying the proposed control technique rather than the other recent control method which again validates the effectiveness of the proposed control method.The reactive powers injected by converters 1 and 2 in case-study 1 while applying the proposed control technique are 92.1 VAr, 92.1 VAr, and on the other hand these powers are 148.2VAr, 53.4 VAr when another recent control method is applied which demonstrates a great enhancement in fair reactive power sharing.In case study 2, the reactive powers injected by converters 1, 2, 3 while applying the proposed control technique are 285 VAr, 285 VAr, 285 VAr and on the other hand, these reactive powers while another recent control method is applied are 5.132 kVAr, −1.292 kVAr, −3.135 kVAr which also validated that the proposed method can successfully share the reactive power among converters.One potential proposal for future research is to integrate this method into an online tuning method for virtual impedances optimization in microgrids.

FIGURE 2 .
FIGURE 2. The proposed control block diagram for a converter in an island MG.

FIGURE 3 .
FIGURE 3. The optimization flowchart for designing virtual impedance of converters in island MG.

FIGURE 5 .
FIGURE 5. Major eigenvalues of the 2-Bus MG in two scenarios, applying the proposed control method (black squares) and the method in[13] (red diamonds).

FIGURE 6 .
FIGURE 6.The comparison of 2-Bus MG frequencies applying the proposed control strategy (full lines) and the method in[13] (dashed lines).

FIGURE 7 .
FIGURE 7. Output powers of converters in the 2-Bus MG, applying the proposed control strategy (full lines) and the method in[13] (dashed lines).

FIGURE 8 .
FIGURE 8. Converters terminal currents in 2-Bus MG applying the proposed control strategy (full lines) and the method in [13] (dashed lines).

FIGURE 9 .
FIGURE 9. Converters output currents in 2-Bus MG applying the proposed control strategy (full lines) and the method in[13] (dashed lines).

FIGURE 10 .
FIGURE 10.Converters terminal voltages in 2-Bus MG applying the proposed control strategy (full lines) and the method in [13] (dashed lines).

TABLE 4 .
Parameters of converters and controllers in the 3-Bus microgrid.

FIGURE 13 .
FIGURE 13.The comparison of 3-Bus MG frequencies applying the proposed control strategy (full lines) and the method in [13] (dashed lines).

FIGURE 14 .
FIGURE 14.Output powers of converters in the 3-Bus MG, applying the proposed control strategy (full lines) and the method in[13] (dashed lines).

FIGURE 15 .
FIGURE 15.Converters output currents in 3-Bus MG applying the proposed control strategy (full lines) and the method in [13] (dashed lines).

FIGURE 16 .
FIGURE 16.Converters terminal voltages in 3-Bus MG applying the proposed control strategy (full lines) and the method in [13] (dashed lines).

385. 1 V
after load change at Time = 2 sec, respectively.Both control methods keep the voltage in a permitted interval after the load-change scenario.The point is that the cost of having zero reactive powers exchanges and enhancing the MG small-signal is to have voltage drops which are mitigated by changing the voltage set-points.

TABLE 1 .
The optimization algorithm specifications.
. The P 1 new which is the active power injected by converter 1 (red line in the