Simultaneous Optimization of Virtual Synchronous Generators Parameters and Virtual Impedances in Islanded Microgrids

An islanded microgrid (MG) including low-inertia converter-based distributed generations (DGs) is subjected to instability. The virtual inertia concept was proposed to alleviate the stability issues by imitating the synchronous generators behavior. This paper spotlights on the optimization of virtual synchronous generator (VSG) parameters and virtual impedances (VI) in islanded MGs using particle swarm optimization (PSO). A small-signal model for MG is developed at first. The permissible ranges of virtual inertia (J) and virtual damping (D) based on MG small-signal stability are scrutinized afterwards. Moreover, VI are considered to lower the reactive power mismatch between converters. Finally, considering the permitted intervals for these parameters, an optimization method and objective function are defined to calculate VSG parameters and VI in the islanded MG. The proposed optimization method enhances the small-signal stability of the MG, decreases the current overshoot and minimizes reactive power mismatches. Simulation results drawn by the “VSG + VI” control include three scenarios. The effectiveness of the proposed “VSG + VI” control method in comparison with “droop” control, “droop + VI”, “non-optimal VSG + VI”, and “VSG ” is verified through simulation studies.


RES
Renewable carbon-free power systems. Two different modes of operation, namely grid-connected and islanded, have been defined for MGs. In grid-connected mode of operation, the gridsupporting role of MG is notable. However, islanded MGs typically supply special loads such as distant rural costumers and cruise ships. The renewable energy sources (RES) are connected to the MG via power electronic converters. However, the lack of inertia in converter-based DGs comparing to conventional machine-based generations can menace the stability of MGs [1], so the virtual synchronous machinery was introduced as a breakthrough [2].
Different control methods to embed virtual inertias in MGs including synchronous generator model based, swing equation based, frequency-power response based, and droop based approaches have been suggested [3]- [4]. A thorough comparison between VSG control method and rate of change of frequency (RoCoF)-based droop control from the mathematical, output characteristics and small-signal stability point of views is found in [5].
The values of the virtual inertias and virtual dampings of an individual VSG must be chosen within stable intervals [6]. Moreover, VSGs can facilitate the damping of low-frequency oscillations created as a result of interactions among synchronous machines [7]. It has been reported in [8] that a stable feedback controller based on robust D-stable region of VSG can enhance the transient response of VSG. However, the reactive power sharing issue in resistive-inductive MGs has not been addressed in this research.
It has been reported that choosing alternating virtual inertias and virtual dampings based on MG frequency deviation and RoCoF can enhance the MG frequency stability [9]. Moreover, the adaptive virtual inertias and virtual dampings can facilitate the power and frequency damping and enhance the dynamic behavior of VSGs [10].
A dual adaptive virtual inertia control strategy has been suggested in [11] to accomplish a balance between active power regulation and frequency regulation, however, the MG reactive power sharing problem has not been resolved in this work. Reference [12] proposes that VSG controllers including governor dynamics can demonstrate power decoupling effect and improve the performance of MG.
An enhanced virtual inertia control strategy has been introduced in [13] to achieve a precise steady state active power regulation and better dynamic response, simultaneously. However, the small-signal stability of MG and reactive power sharing have not been discussed, which is a weakpoint.
Moreover, it has been reported that removing the phase locked loop (PLL) in an advanced VSG controller could positively affect the MG stability [14].
In order to deal with the contradiction between the moment of inertia and frequency response speed, a coordinated adaptive moment of inertia and VI control strategy has been designed in [15] which accelarates the frequency response speed. However, the small-signal stability of MG, optimization of VI and VSG have not been addressed.
The parallel operation of virtual synchronous machines in islanded and grid-connected applications has been analyzed in [16]. However, fair reactive power sharing has not been achieved.
An adaptive data-driven optimal VSG control method has been proposed in [17] by reinforcement learning which limits the frequency deviation and damps the frequency and power oscillations. However, the small signal stability of the MG and reactive power sharing have not been optimized.
A fuzzy-augmented VSG controller has been introduced in [18] which increases the inertia during transients to damp the perturbation. The fuzzy control adds a power correction term to the governer's output power during transients. However, the small-signal stability of the MG and the reactive power sharing have not been considered.
It has been reported that a multi-VSG damping controller using residue index and hybrid PSO algorithm can improve the damping characteristics of multi-VSG grid connected systems [19]. However, the previous research works haven't deal with the reactive power sharing issue. The direct tuning of synchronverter parameters has been performed in [20] to achieve the desired dynamic response.
The other major elements of MGs are VI which have been mostly applied to correct the power sharing and to deliver desired bus voltages. The virtual complex impedance application in [21] aims to enhance the voltage quality and to correct the power sharing. However, the virtual inertia and damping adjustment have not been discussed in the previous work.
The droop equivalent impedance has been introduced in [22] to correct the reactive power sharing using VI. This work has not consider multi-bus complex MG structure and the virtual inertia and damping tuning have not been dealt with.
The adaptive VI in [23] has been proposed to correct the active and reactive power sharing in meshed MG; however, it lacks any MG stability analysis practices. The previous researches such as [24] has provided more accurate modeling for the MG including the PLL dynamic equations. However, it has not considered reactive power sharing, virtual inertia, and virtual damping.
On the other hand, the MG stability is a critical concept considering the small size and the uncertainty of generating units. It includes the primary voltage and frequency control, power sharing, and islanding detection [25]. The MG smallsignal modeling and stability analysis in [26] paved the way for more advanced researches in the field.
The fair reactive power sharing has been attained by adaptive VI in [27] and the method requires communication links and has not analyzed the MG stability. The transient behavior of converters and the small-signal stability of islanded MG have been enhanced using the internal model-based controllers in [28]. The recent research [29] has proposed an optimization algorithm to tune the VI in islanded droopcontrolled MGs to enhance MG stability and reactive powers sharing.
The small-signal stability analysis under distributed control method have been implemented in [30] which requires a communication link between DGs. The cascading of lead compensators have been proposed in [31] to enhance the stability margin of conventional droop control. The small-signal stability of interconnected MG clusters have been analyzed in [32] which has been controlled by a two layer distributed control.
The feasible ranges of VI in a droop-controlled MG have been determined in [33]; however, the PLL, virtual inertia, and virtual damping have not been considered in the model. The secondary distributed voltage control of MGs has been implemented in [34] by input-output feedback linearization. This method enhances the MG stability but requires a complex communication link which is a weakpoint.
Based on the previous literature review, these research gaps are spotlighted.
• The small-signal models of converter-based MG did not include VSGs, VI, and PLL, simultaneously. The MG stability can be violated by the aggregated model of these elements. • The previous literature to draw permissible ranges of VSG parameters have not considered the VI while tuning the VSG parameters and vice versa. The inappropriate tuning of VSG and VI could jeopardize MG stability. • A generalized optimization method to calculate the optimal values of VI and VSG parameters simultaneously could not be found in literature. Considering the previous research gaps, this research is aiming to present the following contributions.
• The permissible ranges of virtual inertia and virtual damping and simultaneous effect of them on MG stability are analyzed considering a thorough MG model including VSGs, loads, lines, PLL and control loops, simultaneously. • The current study proposes a novel PSO-based optimization algorithm and a new objective function to tune the VSG parameters and VI simultaneously to enhance the MG stability while keeping the reactive power mismatches in the minimum level. The proposed objective function excels the multi-objective and the fractional objective functions which suffer from dual problem of assigning weighting factors and instability near the root of denominator, respectively. • The proposed optimization algorithm applies eigenvalues analysis to draw the dominant and non-dominant eigenvalues of MG. The damping ratio of the weakest dominant eigenvalue is used as a measure of microgrid small-signal stability. Moreover, any solution which causes to have an unstable eigenvalue (with positive realpart) is removed by the algorithm. The application of eigenvalues analysis and the damping ratio of dominant eigenvalues (as a measure of small-signal stability) in tuning of VSG and VI in islanded MG has been a novel contribution. • The proposed optimization scheme analyzes all operating points and tunes the VSG parameters and VI, simultaneously. The main positive point is that the small-signal stability of the whole MG (not only the control loops) in all operating points is analyzed and therefore this method can ensure the MG stability. The rest of this paper is organized as follows; The dynamic modeling of islanded MG is presented in Section II. Section III explains the proposed "VSG + VI" optimization  method. The simulation results are presented in Section IV. Finally, Section V terminates this study by presenting the conclusions.

II. THE DYNAMIC MODELING OF ISLANDED MG
The dynamic modeling of MG includes the modeling of all elements on a common dq reference frame. The converterbased DGs in Fig. 1 have controllers as depicted in Fig. 2. The dynamic models of different blocks deployed in Fig. 2 are explained hereafter. An islanded MG is modeled as a nonlinear system in state-space. Then the system equations are linearized around an stable operating point and the small-signal stability is analyzed by eigenvalues analysis [26]- [29] and considering the Lyapunov's first method [35], the permissible ranges for parameters are drawn so as the system remains asymptotically stable. The abc/dq transformation is explained in [29] and is not repeated here because of page limit.

A. Power Calculator and Low-Pass Filter
The instantaneous active and reactive powers injected by a converter (p, q) pass through a low-pass filter (LPF) with the cut-off frequency of ω c and yield active and reactive powers, respectively. The dynamical equations describing this subsystem are as follows:

B. Virtual Impedances
The reactive power loop of a VSG block assigns dq voltage references of the converter based on voltage droop control. The virtual resistances and virtual reactances will be calculated by an optimization method in Section III. The dq voltage references are dictated to the voltage controller as follows:

C. VSG Controller
The VSG block in the converter control diagram is supposed to assign the frequency reference in dq reference frame based on following dynamic equations [38].
Any converter has its frequency reference ω * based on its active power, J, and D values; where ω n and P 0 are nominal frequency and active power of the converter, respectively.

D. Current Controller
Two PI-controllers translate current commands to voltage commands in dq reference frame based on (5) and (6). These voltage commands will then be translated to the abc reference frame by a dq/abc transformation. It is noteworthy to say that γ d and γ q are auxiliary state variables that are used to simplify the state space representation of the system. The decoupling terms which are related to mutual effect of d and q axes are seen which have ω n .L f coefficients.

E. Voltage Controller
The current commands are assigned by voltage controller in Fig. 2. The PLL unit sends the measured frequency to the voltage controller. The q-axis voltage is accordingly controlled. The reason behind definition of γ d and γ q in current controller holds for φ d and φ q in voltage controller.φ

F. PLL Model
The PLL model used in this study was proposed in [24] is demonstrated in Fig. 3 and its dynamic equations are written as follows:

III. THE PROPOSED "VSG + VI" OPTIMIZATION METHOD
The reactive power sharing, the bus voltages, and the stability of MG are affected seriously by the values of VI and VSG parameters. The application of VSG leads to inertial behavior of inverters and this inertial performance can be controlled by J and D. On the other hand, the application of VI could correct the impedance mismatches and enhance the reactive power sharing. However, both VSG and VI will affect the MG eigenvalues and small-signal stability. consequently, simultaneous application and optimization of VSG and VI can cause a better inertial performance, enhanced small-signal stability and fair reactive power sharing. Fig. 4 demonstrates the role of "VSG + VI" algorithm in control of an islanded MG including p VSGs, q lines, and r loads. This novel control method entitled as "VSG + VI" is aiming at calculating J, D, and VI within their corresponding stable intervals for all VSG units, simultaneously. The PSO algorithm is one of the powerful heuristic optimization algorithms which is applied to solve different problems in power systems with a high number of variables to be optimized. The lower number of tuning parameters, the favorable global convergence, and simple tuning of PSO parameters are our motivations to apply this algorithm [36]- [37]. The individual "VSG + VI" method and the "MG small-signal stability analysis" which is a necessary procedure for the "VSG + VI" method are explained in separate subsections as follows.

A. The Individual "VSG + VI" Method Description
Each VSG of a MG has its J, D, virtual resistance (R v ), and virtual inductance (L v ) parameters. An objective function is defined to calculate these parameters to minimize reactive power mismatches and the damping ratio of critical damping point (a dominant eigenvalue with the lowest damping ratio) is maximized. It is suggested that the dominant eigenvalue with minimum damping is the critical eigenvalue which can be a representative of the stability level of the MG. The damping ratio of an eigenvalue which is a complex number (Re + j.Im) is defined by (13) in terms of the real and imaginary parts of it. The more the real part of an eigenvalue, the more the damping ratio of this eigenvalue [35].
The MG at any operating point has n dominant eigenvalues and s non-dominant eigenvalues as seen in (14). The minimum damping ratio among dominant eigenvalues at operating point i is named as ζ i min as seen in (15).
The more the damping ratio of the weakest dominant eigenvalue of MG, the less the percentage overshoot of system response.
where the dominant eigenvalues are those which have small real-parts and the damping of the MG is prominently determined by these eigenvalues. Since there is not a strict definition for dominant eigenvalues, the eigenvalues with real part in a certain interval (−300 < Real(ζ d ) < 0) are considered as dominant eigenvalues [22], [23]. Considering all m operating points of the MG, the total minimum damping of the MG is defined by (16).
The summation of the per unit reactive power mismatches between p converters is calculated by (17).
The proposed objective function is defined in (18) which aims at minimizing the reactive power mismatches and enhancing the small-signal stability of the MG, simultaneously. The voltage drops on buses should be limited to a standard value which is stated by a condition as follows.
Two per-unit terms are multiplied in the objective function. This objective function does not require weighting coefficients and facilitates reaching at two objectives, simultaneously. The flowchart drawn in Fig. 5 describes different steps of individual "VSG + VI" method as follows.  Table I. The optimization variables are considered as particles in PSO, as represented by (19).
The locations of PSO particles are updated based on (20) where the velocity of this movement is calculated by (21) [36].
where y are the location of i th particle in iteration k + 1, the location of i th particle in iteration k, and the velocity of i th particle in iteration k + 1, respectively.
i is the value of objective function for i th particle in k th iteration which is calculated by (18) and P (k) best is the minimum value of objective function experienced by all population in k th iteration. The constant values ω, r 1 , r 2 , c 1 , and c 2 are found in Table I and explained in Nomenclature. The locations which can be occupied by particles are constrained by minimum (y i,min ) and maximum (y i,max ) permissible values. This condition can be summarized by (22). It is assumed that i = 1. • Step 5: This step is shown as a condition block in the flowchart. If all eigenvalues at operating point x i are stable, the ζ i min and Q(p.u.) are reported to the data storage (step 6). Otherwise, if one or several eigenvalues have positive real parts flag = 1 is set and the algorithm goes to step 7. Moreover, since that PSO solution is not realizable, a large value (from the order of 10000) is considered for the objective function so as PSO automatically ignores this solution. • Step 6: The ζ i min for all operating points are calculated by (15) and saved in this data storage and the ζ min is chosen as the minimum value among them as seen in (16). The corresponding Q(p.u.) is also saved in this data storage. • Step 7: The objective function is calculated considering the data available in step 6 and also the flag received from step 5. If the flag is 1 the PSO removes this unstable solution by considering a large value for objective function. • Step 8: If all the operating points have been analyzed, the algorithm proceeds to step 9, otherwise the next operating point is analyzed starting from step 4. • Step 9: Checking the convergence criterion which can be the value of objective function or the number of iterations (it). If the algorithm is not converged, the next iteration will begin (it = it + 1). • Step 10: The optimal values for J, D and VI are exported by the algorithm and these values are applied to the VSG controllers.

B. MG Small-Signal Stability Analysis
The small-signal stability has been examined by Lyapunov's first method in this study [35]. The small-signal stability of a non-linear system is examined by the eigenvalues locations.
1) When the eigenvalues have negative real parts, the original system is asymptotically stable. 2) When at least one of the eigenvalues has a positive real part, the original system is unstable. 3) When the eigenvalues have real parts equal to zero, it is not possible on the basis of the first approximation to say anything in general. Since the islanded MG have one zero eigenvalue [24]- [26], it is asymptotically stable if the first and second conditions are met.
The "VSG + VI" optimization method requires a smallsignal stability analysis to draw the eigenvalues and corresponding damping ratios. A MG includes several VSGs, lines, and loads. In order to draw the linearized state-space equations of the whole MG, the state equations of all VSGs, loads, and lines be written and considering all these state-equations, the A MG matrix is drawn straightforwardly.
A typical MG including arbitrary number of VSGs is considered as seen in Fig. 6. It is assumed that a transmission line is connecting buses k and i and another line connects buses i and j. The state variables of VSGs i and j are chosen as seen in (23) and (24), respectively.
The reference frame of VSG 1 is chosen as the common reference frame. The state of common VSG is explained by (25).
The state variables of the line between buses i and j, and the line between buses k and i and the load at bus i are chosen as seen in (26), (27), (28), respectively. It should be noted that these variables locate on common DQ reference frame.
The typical bus i is considered in an arbitrary islanded MG including p VSGs, q lines, and r loads. The linearized state space equations are drawn around an stable operating point for this VSG and the loads and lines connected to it as it is seen in (29) and (30).
where A 1 -A 14 matrices are explained in the Appendix. The state equations describing a transmission line between buses i and j are written as (31).
The state equations for a static load installed at bus i are seen in (32).
The foregoing formulations explain the state equations for all VSGs, lines, and loads. The final step is to consider the MG as a united system such as (33) and draw the state matrix.
x MG = x VSG,1 , . . . , x VSG,p , x line,1 , . . . , x line,q , A final and straightforward step is to add zero rows and columns to any of A 1 -A 14 matrices in (29)-(32) to express these equations again in terms of MG state variables matrix (33). For example, in a 2 bus MG including two VSGs, two loads and one line, there are 38 state variables. Consequently the matrix x MG is 38 × 1. If the state equations for VSG 1 are going to be written by (29) and (30), to write the equations in terms of MG variables, all A i matrices are rewritten in a suitable way by adding zero rows and columns so as they will finally have 38 columns. For instance, the matrix A o 1 while applying (29) The equation (30) is similarly rewritten using where A o 3 -A o 6 are drawn by the previous zero rows and columns insertion method. The matrix A oo 6 is drawn similar to A o 6 from A 6 and its difference with A o 6 is that their non-zero columns and rows locate in different rows and columns. The same argument holds for A oo 9 , A oo 10 , and A oo 14 . The general state-space equations for VSG i , line ij , and load i are seen in (37), (38), and (39), respectively.
When all state equations for all VSGs, lines, and loads are drawn and written in general form in terms of x MG , then the state matrix of MG can be drawn by (40). The state space representation of MG including p VSGs, q lines and r loads is written as (41).
A MATLAB command (eig) can be used to draw the eigenvalues of the MG as follows.

IV. SIMULATION RESULTS
The simulation studies are performed in three separate scenarios. Scenarios 1 and 2 are performed in a 2-bus MG [24] and scenario 3 is applied to a 3-bus MG in [26].

A. Scenario 1: Comparison of Proposed "VSG + VI" Control With "Droop" Control and "Droop + VI" in 2-Bus MG
A load perturbation is enforced at t = 1.8 s and lasts until t = 5 s. The series impedance R pert +jω.L pert is switched-in at bus 1 in parallel with R load1 + jω.L load1 as it is seen in Fig. 1. The load at bus 2 is R load2 + jω.L load2 . The load perturbation is 440 W + 48 VAR at bus 1. Table I and Table II show detailed parameters of the MG and the values obtained for J, D, and VI by the proposed "VSG + VI" method. Firstly, a VSG control strategy which was introduced in previous section is applied to the converters. The effect of changing J individually, D individually, and simultaneous J and D are analyzed. The permissible ranges of J and D are drawn accordingly. Afterwards, the optimal "VSG + VI" method introduced in Fig. 4 is used to draw the optimal virtual inertias, virtual dampings, and VI for two converters in the MG. In scenario 1, the results for the optimal "VSG + VI" control are compared to the results drawn from the "droop" control [24] and optimal "droop + VI".

1) Virtual Inertias and Virtual Dampings:
The MG can work stable or unstable depending on the values of virtual inertias and virtual dampings. The virtual inertias of inverters are changing from J = 0.1 to J = 10 and the virtual dampings are fixed at D = 100 the major eigenvalues of the MG are drawn in Fig. 7. The more the virtual inertias are increased, the more the MG approaches the instability boundary.
On the other hand, while keeping a fixed virtual inertia (J = 1), the virtual dampings are changing from D = 200 to D = 10 to analyze its effect on the MG's major eigenvalues. As it is seen in Fig. 8, lowering the virtual dampings decreased the real value of major eigenvalues and propelled the MG towards the instability. For D = 10 the MG is in the boundary of instability. However, the simultaneous change of     10. The comparison of MG frequency in load change scenario applying the proposed optimal "VSG + VI" control, droop control [24], and "droop + VI". virtual inertia and virtual damping will push the critical eigenvalue either toward instability (over the black plain) or stable region (beneath the black plain) of operation as it is depicted in Fig. 9.
2) MG Frequency: The MG frequency is demonstrated in Fig. 10 for three different cases controlled by the proposed optimal "VSG + VI" control method, droop control method [24], and "droop + VI" method. As it is seen, the RoCoF for these three cases are 0.351 Hz/s, 0.61 Hz/s, and 0.61 Hz/s, respectively. The lower RoCoF is preferred in order to avoid unnecessary relay trippings. So, the proposed optimal "VSG + VI" method excels the other control methods in this regard. The minimum frequency point (Nadir), while applying "VSG + VI" control method is higher than the corresponding value while using droop control in [24] and "droop + VI". All in all, the "VSG + VI" method demonstrates a better frequency control performance in comparison with the droop control [24] and "droop + VI" method. The enhanced frequency control by the proposed optimal "VSG + VI" method can support the frequency stability of MG effectively.

3) MG Voltage Components:
The output dq voltages of converters are seen in Fig. 11 for three different cases controlled by the proposed optimal "VSG + VI" control method, the "droop" control, and "droop + VI" method. The d-axis voltages, which are seen in the upper part of Fig. 11, are supposed to remain at zero and the voltage oscillation at t = 1.8 s for "VSG + VI" control method is less than the corresponding value using "droop" control method and "droop + VI" method. The bottom part of Fig. 11 shows the q-axis voltages. The voltage control loop in "Droop + VI" control method includes the VI, therefore the voltage drops for converters using this method are greater than voltage drops of "droop" control as it can be seen in (3). The maximum q-axis voltage drop for "VSG + VI" method, "droop" control, and "Droop + VI" control are 1.02 V, 1.16 V, and 1.98 V, respectively. However, the "VSG + VI" control method excels another two methods and the steady state voltage drop is negligible.

4) Major Eigenvalues of the MG:
The eigenvalues of the MG applying the "droop" control [24], the proposed "VSG + VI" control, and "droop + VI" control are demonstrated in Fig. 12. The dominant eigenvalues have been enhanced by "VSG + VI" optimization method in Fig. 12 in comparison to "droop" control and "droop + VI" method. For example in "droop + VI" method the dominant eigenvalue λ d = −2.43 ± j15.01 (ζ d = 15.98%) in Fig. 12b is displaced to λ d = −0.42 (ζ d = 100%) using the proposed "VSG + VI" control; however, the corresponding eigenvalue for "droop" method is λ d = −1.45 ± j5.39 (ζ d = 25.97%) which is weak and less damped. The other dominant eigenvalues also are enhanced and more damped by "VSG + VI" optimization method in comparison with the other methods. Afterwards, by applying (13), the damping ratios of all eigenvalues are calculated and expressed as p.u. values.  The bar diagram in Fig. 13 demonstrates the damping ratios (ζ ) for MG eigenvalues in three different scenarios; Firstly, applying a "droop" control in [24] (black bar diagram) and then applying the proposed "VSG + VI" control method (green bar diagram), and finally using the "droop + VI" control method (red bar diagram). Analyzing Fig. 13 unveils that the minimum damping ratio of dominant eigenvalues while applying "droop" control, optimal "VSG + VI" control, and "droop + VI" control are 26%, 28.7%, and 18, 6%, respectively. Consequently, the "VSG + VI" control method excels the other methods from the MG small-signal stability point of view.

5) Eigenvalues Participation Analysis:
A participation analysis can reveal which states more affect a certain mode. In this regards, a participation matrix is defined for all states which was well-explained in [35]. Considering the right and left eigenvectors as φ and ψ, the participation matrix is defined as: where the i th participation vector is defined as: where p ni represents the effect of n th state variable on i th mode. Since the stability can be examined by major eigenvalues, the low frequency modes and their major participants are demonstrated in Table III. This table shows the modes with lowest damping ratios (ζ ), major participants which affect these modes, and their corresponding participation factors. For the modes 13 and 14 which are the worst eigenvalues, the major participants are v od1 , v od2 , v oq1 , v oq2 , v od,f 1 , and v od,f 2 .

6) Output Current Components of Converters:
The current components of two converters in the common dq reference frame are depicted in Fig. 14. The current components while applying traditional "droop" control [24], the proposed "VSG + VI" control, and "droop + VI" control method are plotted simultaneously to facilitate the comparison.
The d-components (i od1 and i od2 ) while using "droop" control have different steady-state values (1.166 and 0.421 A) but the corresponding i od1 and i od2 values when deploying the proposed "VSG + VI" control method and "droop + VI" control method are 0.736, 0.853, 0.71, and 0.85 A, respectively. It means that the proposed "VSG + VI" control excels the other methods to share the d-component of the current among converters. It should be noted the performance of "Droop control + VI" method has been better than "droop" control [24] but it includes more current fluctuations.
On the other hand, the q-components, injected by two converters using "droop" control, "VSG + VI" control, and "droop + VI" control eventually reach at identical values (i oq1 = i oq2 = 5.00 A), but the inertial behavior of "VSG + VI" method causes the q-current to have a bigger settling time (green curves in Fig. 14). Moreover, the "droop + VI" has the shortest settling time, however it includes adverse current fluctuations. Fig. 15 demonstrates the active and reactive powers injected by converters 1 and 2 to the MG. The active powers injected by converters 1 and 2 reach at identical values (637.4 W) because two converters have similar specifications while applying "droop" control, "VSG + VI" control, and "droop + VI" control, respectively. The active power while using "droop + VI" has a short settling time which excels the "droop" control [24]. However, when applying "VSG + VI" control, the active power indicates inertial behavior and the settling time is longer than "Droop + VI" control method because of virtual inertia and damping. The "Droop + VI" control method has no virtual inertia and damping in its control loops and a faster response is expected from it. The damping of "VSG + VI" control method can be enhanced by tuning of control parameters such as PI controller coefficients.

7) Output Active and Reactive Powers of Converters:
The lower part of Fig. 15 shows the reactive powers of converters 1 and 2 while applying conventional "droop" control [24], "VSG + VI" control, and "droop + VI" control, respectively. As it is seen, the reactive powers injected by converters 1 and 2 while using "droop" control are totally different (53.6, 148 VAR), but applying "VSG + VI" control method leads to roughly identical reactive power injections by converters 1 and 2 (96.30, 108.93 VAR). The steady-state reactive powers of converters while applying "droop + VI" are 90.27 and 108.21 VAR. Consequently, the "VSG + VI" control method excels the conventional "droop" control [24] and "droop + VI" in reactive power sharing.

B. Scenario 2: Comparison of "VSG + VI" Control With "Non-Optimal VSG + VI" and "Non-Optimal VSG" in 2-Bus MG
The simulation studies are performed on a 2-bus MG introduced in [24]. A load perturbation is enforced at t = 1.8 s and lasts until t = 10 s. The series impedance R pert + jω.L pert is switched-in at bus 1 in parallel with R load1 + jω.L load1 as it is seen in Fig. 1. The load at bus 2 is R load2 + jω.L load2 .
In this scenario the effect of optimal parameters (J, D, Lv 1 , Lv 2 , Rv 1 , Rv 2 ) on dynamic response of VSG units in the 2-bus test MG is analyzed in comparison with two non-optimal VSGs. The load change event is considered identical to scenario 1 and different dynamic characteristics of the MG while   having optimal VSG (case 1) are compared with two different cases; non-optimal VSG + VI (case 2) and non-optimal VSG (case 3) while the parameters are seen in Table IV. Different dynamic characteristics are analyzed hereafter.

1) MG Voltage Components:
The voltage components of bus 1 of the MG are depicted in Fig. 16. As it is seen, v od in cases 1-3 show similar small deviations and its steady-state value is zero which is enforced by control methods in all three cases. The q-axis voltages are shown in bottom part of Fig. 16. As it is seen, the steady-state voltage drops in cases 1-3 are 0.99 V, 5.96 V, and 0.07 V, respectively. It means adding VI makes voltage drops 1.15% and 6.9% in cases 1 and 2, respectively. Moreover, the voltage drop in case 1 is less than 5% and it is permissible. The voltage drop in case 3 is negligible, but the other dynamic characteristics of VSGs are not favorable in this case.
2) Eigenvalues of the MG: The minimum damping ratios of 2-bus MG in different operating points have been shown in Fig. 17. Minimum damping ratios in cases 1, 2, and 3 are 15.38%, 15.35%, and 15.38%, respectively. However, the steady state damping ratio of case 1 (i = 10000) is more than corresponding value in case 2 and case 3. Consequently, the proposed optimal VSG exceeds non-optimal VSGs. The dominant eigenvalues of 2-bus MG in three operating points (A, B, C) are shown in Fig. 18. The MG in all three cases has one eigenvalue at (0, 0) and the other eigenvalues have  negative real parts and therefore the MG is asymptotically stable [35]. As it is seen, the dominant eigenvalues in cases 1, 2, and 3 are −6.65±22.18, −6.61±22.25, and −6.66±22.18. The damping ratios are 28.73%, 28.47%, and 28.75%, respectively. Consequently, the optimal VSG (case 1) and case 3 have been roughly similar in enhancing the MG small-signal stability at the operating points A and C.

3) Output Current Components of VSGs:
The current components of VSGs in the MG are shown in Fig. 19.  of view and case 1 has a current mismatch comparable to case 2. The lower current mismatch decreases the pressure on switches which is of our interest.
The q axis current components are shown in bottom graph of Fig. 19. The overshoot of current in case 3 is bigger than the other cases and settling time of case 1 is shorter than the other cases. However, the bigger virtual inertial of case 1 makes its current response inertial. The VSGs in case 2 inject the lowest q axis current to the MG which is because of the highest voltage drops in this case. The active power response of VSGs in case 2 is not favorable and reactive power mismatches of case 3 is not preferred. Therefore the case 1 which is the optimal case has the best performance.

C. Scenario 3: Comparison of "VSG + VI" Control With "VSG" Control in 3-Bus MG
In this scenario, the 3-bus test MG in Fig. 21 is implemented in MATLAB. The MG and VSGs specifications are listed in Table V and Table VI, respectively. In this scenario the load R load3 + jX load3 is switched in at t = 1 s, and the load R load1 + jX load1 is switched in at t = 3 s. The MG characteristics including eigenvalues, voltage components, current components, and active and reactive power are scrutinized herein.
1) Major Eigenvalues of the MG: The eigenvalues of the 3-bus MG are depicted in Fig. 22 at operating point O. The whole eigenvalues in Fig. 22a locate in left half plain (except   SCENARIO 3 one in (0, 0)). Therefore, the MG is asymptotically stable while it is controlled by VSG or "VSG + VI" method.
Focusing on Fig. 22b discloses that the weakest dominant eigenvalue while using "VSG + VI" method has a better damping ratio and overall stability. For instance, the dominant eigenvalue λ d min,1 while using "VSG" control is −4.83 ± 8.17, however, the dominant eigenvalue λ d min,2 while applying "VSG + VI" is −7.25 ± 11.69. Consequently, the damping ratio is enhanced from 50.87% to 52.70% using the "VSG + VI" control method. The non-dominant stable eigenvalues have minor effects on the MG small-signal stability and can be ignored.
The bar diagram in Fig. 23 demonstrates the damping ratios (ζ ) for all MG eigenvalues. The minimum damping ratio while applying "VSG" control and "VSG + VI" control are 50.87% and 52.70%, respectively. Therefore, the proposed "VSG + VI" control enhance the damping ratio of the weakest dominant eigenvalue.
2) MG Voltage Components: The output dq voltage components of VSGs are seen in Fig. 24. The d-axis voltages are meant to become zero in both control methods and Fig. 24a verifies that either using "VSG" or "VSG + VI" the d-axis   voltages are zero, as expected. The q-axis voltages while using "VSG" or "VSG + VI" method should satisfy the voltage drop constraint (18). As it is seen in Fig. 24b, the steady-state voltage drops while using "VSG" and "VSG + VI" methods are 1.36 V and 3.20 V, respectively. In other words, the q-axis output voltages while applying "VSG + VI" method and "VSG" method are 0.97 p.u. and 0.99 p.u.. However, the voltage drop while applying "VSG + VI" method is more than the voltage drop using "VSG" method, but still the percentage of voltage drops in both methods are acceptable. Theses voltage drops can be corrected in secondary voltage control by changing the voltage set-points.

3) Output Current Components of VSGs:
The current components of three VSGs in the common dq reference frame are depicted in Fig. 25. The d-axis currents are shown in Fig. 25a. The maximum d-axis currents overshoot while using "VSG" control and "VSG + VI" control are 2 A and 0.84 A, respectively. On the other hand, the d-axis current mismatches among VSGs are 2.46 A and 0.36 A while using "VSG" and "VSG + VI", respectively. The lower current overshoot and current mismatches in d-axis are prominent advantages for proposed "VSG + VI" method.
The q-axis current components are demonstrated in Fig. 25b in load change scenario. The steady-state q-axis currents while using "VSG" or "VSG + VI" method have identical values. However, the settling time while applying "VSG + VI" method is clearly shorter than the case with "VSG" control. For example, the VSG response is settled at t = 5 s and t = 4 s while applying "VSG" control and "VSG + VI" control, respectively. Therefore the "VSG + VI" response is faster than "VSG" response. Moreover, the maximum q-axis current overshoot while applying "VSG" and "VSG + VI" are 8.32 A and 7.74 A, respectively.
On the whole, the faster current response with lower overshoot and smaller current mismatches are the substantial advantages of "VSG + VI" method in comparison with "VSG" control. Fig. 26 depicts the active and reactive powers injected by VSGs to the 3-bus MG. The maximum active power injected by VSG 1 to the MG in Fig. 26a while applying "VSG" and "VSG + VI" control are 4665 W and 4230 W, respectively. Therefore lower active power overshoot is an advantage for "VSG + VI" method. Moreover, the settling time of active power while applying "VSG + VI" method is shorter in Fig. 26a. For instance, P 1 is settled down at t = 5 s and t = 4 s while applying "VSG" and "VSG + VI", respectively. Consequently, active power response is prominently faster while using "VSG + VI" control method.

4) Output Active and Reactive Powers of VSGs:
The reactive powers injected to the 3-bus MG are seen in Fig. 26b. The maximum reactive power injected to the MG while applying "VSG" and "VSG + VI" control methods are 1059.92 and 548.75 VAR, respectively. On the other hand, the maximum reactive power mismatches while applying "VSG" and "VSG + VI" methods are 1203 and 207.25 VAR, respectively. Therefore, smaller reactive power mismatch is an outstanding pros for "VSG + VI" method.
All in all, the "VSG + VI" control method excels "VSG" control by having a shorter settling time, smaller power overshoot, and a negligible reactive power mismatch.

V. CONCLUSION
The present paper develops a small-signal model for a VSG-based MG including PLL unit and VI in the control structure. The permissible intervals for virtual inertias and virtual dampings are analyzed and drawn from the MG small-signal stability point of view. Then a novel PSO-based optimization algorithm is introduced which aims at enhancing small-signal stability of MG and minimizing the reactive power mismatches, simultaneously. The introduced objective function removes the dual problem of adjusting the weighting coefficients in a multi-objective function. The optimization algorithm analyzes the MG stability in all operating points inside the operating time interval as an off-line procedure. The optimization algorithm calculates the optimal values of virtual inertias, virtual dampings and VI so as the damping ratio of the weakest dominant eigenvalue of the MG is enhanced as much as possible. The settling-time of the current and power responses are shortened using the proposed "VSG + VI" control method. Moreover, this method decreases the percentage overshoot of current and active power and succeeds the other control methods in enhancing the MG small-signal stability and reactive power sharing among VSGs.

APPENDIX MG SMALL-SIGNAL STABILITY MATRICES
The following matrices are used in general small-signal stability analysis in Section III-B. It should be noted that the I od,i , I oq,i , I ld,i , I lq,i , V od,i , and V oq,i are the dq components of output current, terminal current, and output voltage of VSG i in the corresponding stable operating point.