Microgrid small-signal stability analysis considering dynamic load model

Microgrids suitable operation depends on the system’s proper design, based on accurate component models. A complete model considering entire system dynamics is a vital necessity for a successful design. Using dynamic load in microgrid small-signal model results in a model that shows transient and steady-state dynamics, since designing a low-inertia system like microgrid need extra accuracy. In this paper, an inverter-based microgrid’s small-signal model in islanded operation mode containing a dynamic load model has been achieved. The Exponential Recovery Load (ERL) model is presented here to study the dynamic behavior of load. Microgrid stability analysis has been carried out for both static and dynamic load models. To attain this goal, the entire component equations are obtained and linearized around an operating point, then the state-space model is formed and eigenvalues are derived. System stability condition has been investigated by plotting eigenvalues and root locus method. In the ﬁnal step, a simulation on Matlab/Simulink is performed to show the microgrids operation in the presence of a static and dynamic load.


INTRODUCTION
Today, with energy consumption increase and related environmental contaminations, the necessity of renewable clean energy sources can be felt. Using renewable distributed energy sources close to load centers not just decrease transmission and distribution costs but also reduce power system dependency to large power plants. A microgrid is a group of distributed sources and loads connected through a local network, that can operate either connected or disconnected to the power grid, respectively, known as grid-connected and islanded operation modes. Along with all the positive effects of the microgrid, there is a lot of design complexity that comes with its utilization. Any proper design requires a complete model that describes a full picture of the microgrid. Therefore, the primary goal is to obtain a precise model of the microgrid.
Most of the researches in this field have been focused on microgrid operating characteristics. Distributed Generations (DGs) and power storages are usually connected to the local network through Inverter-Interfaced DGs (IIDGs) due to the This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2021 The Authors. IET Renewable Power Generation published by John Wiley & Sons Ltd on behalf of The Institution of Engineering and Technology flexibility of operation and control. But using power electronic devices instead of conventional methods, make systems inertia very low [1]. Lack of inertia results in a fast response to any change in microgrids condition [2]. Studies have shown IIDGs should be employed most carefully due to their low-inertia and lack of sufficient synchronizing torque, to stay in sync with other power supplies and rest of power system in islanded and connected operation modes [3]. In [4], a study has been carried out to analyze the impact of inertia on power and frequency. This study proposes a dual-adaptivity inertia control to minimize power and frequency variation.
Stability analysis formed a major section of microgrids studies which divide into small and large signal stability analysis. In small-signal analysis, non-linear system linearized around an operating point, then the stability of the derived model was investigated using linear methods [5].
In [6] and [7], effective control methods on some assumed disturbances have been studied. Each article tries to eliminate disturbances effect by proposing control methods. The most popular and adaptable control method in IIDGs is droop control which originates from the principle of conventional synchronous generators by decreasing frequency/voltage in the response of any active/reactive load power increase [8], [9]. By employing droop control, an accurate power balance between power sources will be established [1]. Consequently, many researches have been carried out based on droop control optimization [2], [10]. Reference [11] proposed a two-level hierarchical controller consisting of a primary decentralized local droop controller in each DG and a centralized secondary controller to compensate load's voltage and frequency variation caused by load disturbances. One of the reasons that cause system instability is the load fluctuation [12]. Therefore, an uninterrupted connection must be established between generation and demand. Reference [13] studied the accurate load model's effect on voltage stability analysis in different power systems. Most researches mainly focused on distributed DC sources combined with inverter interface to supply load demand [2], [14]. It has been proven that load composition can influence MG's stability especially in islanded mode [15]. Therefore, an extra care in load modeling with maximum coverage should be utilized.
Recently a review in [16] has been done to classify load modeling studies. Dynamic models describe active/reactive power as functions of voltage and time [17], on the other hand, static load models are just voltage-dependent or even constant power. The importance of the load model on voltage stability has been carried out in [13], [17]. Almost all studies on microgrid stability such as [1], [2], [10] used static load models and mostly RL models. In [19] active load model has been introduced, active load consists of a resistor as load consumer, current and voltage controller to manage load's power demand, and a controlled rectifier. In 1993 David. J. Hill studied voltage stability and associated load modeling [20]. The proposed load model was designed to capture the usual non-linear behavior of load. Based on laboratory measurements, a typical load response to a step in the voltage is of the exponential form. One year later, Hill and Karlsson research parameter estimation for some typical loads at different times of the day and year [21].
Gathering information about load behavior is a preliminary step in load parameter estimation. Measurement-based and component-based are two approaches to data collection. In [22] measurement-based approach used data collection of a static ZIP model and a dynamic exponential recovery model. Afterward, parameter estimation was performed through an improved particle swarm optimization (PSO) algorithm. Besides model accuracy, one of the most important issues about load models is minimizing computational costs in the parameter estimation process [23]; dynamic load models have this advantage with the Exponential Recovery Load (ERL) model due to its fewer parameters.
Composite load models became popular in recent studies, in [24] a parameter estimation for a composite load consist of ZIP and Induction Motor (IM) models has been carried out. In a similar study, [25] present a modeling approach to estimate the composite load model proposed by the Western Electricity Coordinating Council (WECC) for the first time. In [26] transient response of disturbances in different load models is derived, then simulation and experimental results are com-pared. In [27], IM load model in microgrid stability analysis is presented. Although the dynamics of the load have been considered, lack of generality in the load model seems to be a lag. In a similar study in [28], a state-space model has been derived by linearization of the differential equation of three major sub-modules in the DQ reference frame. In this study, loads are a combination of IM and active load models and showed that output voltage and frequency of microgrid in IM load bus have major oscillations which conduct system to instability.
Utilization of microgrids comes with uncertainties generated power fluctuations in some energy sources such as wind and solar [29]. It's important to study any possible disturbance scenarios; in [30] an inverter-based microgrid under harmonic condition is studied. Using Dynamic Phasor (DP) theory, the fundamental and harmonic component of output ac waveform via dc variables is investigated. The main contribution of this paper is developing a simple droop controller with harmonic decomposer and virtual impedance blocks to compensate harmonics. Another disturbance scenario considered in [31] is a microgrid under unbalanced condition; the DP model obtained then eigenvalues and participation factors were derived. Based on its results, main controller parameters for desirable system performance are determined and performed on Matlab/Simulink simulations.
A dynamic IM, Constant Power Load (CPL) and resistive load models are considered in [32] to study microgrid stability by taking system uncertainties into account. This study showed that the presence of the IM load model affects system dynamics which particularly at low-power demands conducts system to an unstable region based on the eigenvalue spectrum. In the recent work in [33], parameter stability region based on bifurcation theory studied an ac microgrid consisting of a combination of IM and ZIP load model. Based on its results, droop controller parameters stable region is determined. As can be seen, most studies used various load combinations. Even though the IM load model is considered in some references, ERL model is neglected in studies. Therefore, conducting a study considering the ERL model, that can show the dynamics of thermostatic load and also has fewer variables than the IM model as the only corresponding dynamic model, seems to be necessary.
In this paper, we intend to cover the following targets: 1. State-space modeling of dynamic ERL model is derived, and its effect on the dynamic-stability of an islanded microgrid is investigated. The ERL model is more realistic and reasonable for variable and complex nature of the dynamic loads such as electrical heating devices which has not been conducted previously. 2. Based on eigenvalues sensitivity analysis, the effect of system parameters (droop coefficients, load model parameters, network parameters etc.) on the stability of the whole system is derived in the form of participation factors (PFs), and 3. Controller parameters adjustment according to PFs values is presented in the presence of dynamic load model. It is worth noting that the design of microgrids with static load models while the nature of loads is dynamic may lead to microgrid instability. 4. The sensitivity analysis of the system stability and its parameters to changes in dynamic load parameters have been evaluated.
Remainder of this paper is organized as follows; in Section 2, by presenting all sub-modules state equations including Voltage Source Inverter (VSIs), network and static and dynamic load models, the small signal model of an islanded microgrid is obtained. Afterwards, the eigenvalue analysis of the model and the PF calculating method are introduced. By observing the behavior of the system in the presence of Resistive-Inductive Load (RLL) and ERL models, control objectives are set to be applied in Section 4 using the extracted PFs. Finally, Section 5 contained a discussion about case study results following by a guide to future works.

SMALL-SIGNAL MODELING OF ISLANDED MICROGRID
An islanded microgrid generally divide into three major submodules as shown in Figure 1: VSIs, network and loads. Each inverter is modeled on its individual reference frame. Rotational frequency sets by the inverter power controller. To achieve a complete model on a common reference frame, one of the inverter frames is taken as common, as shown in Figure 1. In this paper, without losing generality, the first inverter frame is selected as the common reference frame. Figure 2 shows the reference frame transformation from each individual frame to a common frame.
Axis (d − q) i , (d − q) j are the individual reference frames of the i th and j th DG unit, respectively, and (D-Q) represents a common reference frame. i is the angle between the i th and common reference frame.
References transfer equations from (d − q) i frame to (D-Q) are as in (1) and (2), [ f dq ] represent any state equation which transferring from any individual (d-q) to common (D-Q) reference frame and [T i ] is the transfer matrix [34]. . (2)

State-space model of VSIs
As shown in Figure 3, each VSI includes a three-leg inverter, output LC filter, coupling inductance, and controller section. The power source is connected to the LC filter through a threephase inverter controlled by the control section. Finally output current goes through coupling inductance. It is to be noted that every equation in this section and consequently state-equations which form the state-space model are for the i th inverter but for keeping briefness, subscription ′′ i ′′ has been ignored unless specification of the inverter is necessary (as in (11) and (12)). In the following, each sub-system is discussed and state-equations are derived.

Controller section
The controller section consists of three subsystems: 1) power sharing controller, 2) voltage controller, and 3) current controller, as shown in Figure 3.

Power sharing controller:
The main idea of the power controller is to imitate the governor of the conventional synchronous generator [1]. The power controller sets the fundamental voltage magnitude and frequency of inverter output voltage, according to droop characteristics of active power-frequency and reactive power-voltage magnitude.
The droop controller responds to any load increase by decreasing frequency and voltage magnitude, as follows (3), (4) are the droop characteristic equations where 0 , V 0 are nominal angular frequency and voltage values corresponding to P 0 , Q 0 , and m p , n q are droop characteristics slopes which are  (5) and (6): where max and min , and V od min and V od min denote, respectively, the maximum and minimum values of angular frequency, and d-axis output voltage on droop characteristic. The instantaneous active and reactive powersp andq are calculated from measured output voltage and currents as (7) and (8).p As shown in Figure 3, calculated powers passed through lowpass filter, resulting in filtered P and Q powers: where c is the cut-off frequency of low-pass filter. As frequency in the inverter is obtained from (3), i will be: And therefore:̇i wherėi is the rotational frequency of i th inverter seen from the reference frame rotating at com . Voltage magnitude obtained from (3) aligned to d-axis (v * od ) and q-axis value is set to zero. Then v * od and are the outputs and , P and Q are the statevariables of power controller block.
Voltage controller: Voltage controller mainly consists of a standard Proportional Integral (PI) controller as shown in Figure 3. Inputs of this block are feedback signals v o , i o from filter's output and power controller's output signal. PI controller compares feedback with a reference signal and create a feed-forward gain to compensate for output current disturbances. d and q are the state-variables of the voltage controller, thus state-equations can be described as: where i * ld and i * lq are block's feed-forward output signals to the current controller, K pv and K iv are proportional and integral gains of the voltage controller, C f is the filter capacitance and G is the feed-forward control gain.
where i ld ,q and i * ld ,q are, respectively, inverters output current and voltage controllers output signals, and v * id ,q are the current controller's feed-forward signal to the inverter. K pc and K ic are proportional and integral gains of the PI controller, respectively, and L f is the filter inductance value.

Inverter
Inverter's function assumed to be ideal, so this block's input signal will be present as its output:

LC filter and coupling inductance
LC filter is implemented to remove high-frequency harmonics caused by the switching operation of the inverter. Stateequations of this block presented here: where r f and L f are resistance and inductance of LC filter, v id ,q are inverter's output voltage, r c and L c are coupling inductance parameters, and v bd ,q are bus voltages.

The complete model of VSIs
A complete model of each VSI can be obtained by combining all sub-system models. State-equations are obtained and by linearizing them, a state-space model is derived: where A inv , B inv and B com are presented in Appendix. As can be seen, there are totally 13 states, three inputs, and three outputs in each inverter model. Inputs consist of bus voltages (v bD,Q ) and common reference rotational frequency ( com ), and outputs are individual inverter's rotational frequency ( i ) and output currents (i oD,Q ) (Chosen reference inverter has two inputs and three outputs).
To describe all parallel function of all VSIs state-space model of each individual inverters is combined as below: and complete state-variables are: The resulting combinations A INV , B INV and C INVc are given in Appendix.

State-space model of the network
The network is considered as an overhead line, thereupon modeled as a series RL as shown in Figure 1. State-equations of the line between buses i j and i j on global reference frame are [1], [18]: By linearizing (34) and (35): where R line and L line are, respectively, lines resistance and inductance, and v bD and v bQ are bus voltages, and i lineD and i lineQ are line currents, in D − Q frame. It is important to be noted that if the network was of underground, line capacitance was considerable. Consequently state equations of the network were different and had a different impact on system stability.

State-space model of load
Load models used in other stability analysis are mainly a simple series resistance and inductance (i.e. RL); this static constant impedance model adds no dynamic to the small-signal model, same as the network model. In the designing process, neglecting any dynamic drag entire system to instability and without suitable control, even voltage collapse. Here we implement two different load models in the same microgrid, to see how a less accurate load model can result in microgrid instability.

Static constant impedance load (CIL) model
As mentioned before, this model has no independent dynamic and two considered dynamic are under the influence of network currents dynamics. State variables of constant impedance are load currents (i loadD,Q ): Linearized state equations are:

Dynamic ERL model
Having an accurate picture congruous with actual operating load data is a vital prerequisite for the designing process. Dynamic models give the designer this accurate picture during the identification process which performs to identify load parameters. The ERL model is presented as a set of non-linear relations between active/reactive power and voltage variations: where v 0 , P 0 , and Q 0 are pre-disturbance voltage and power values, P r (t ) and Q r (t ) are power recovery of the load, P L (t ) and Q L (t ) are total load response to voltage variation, p and q are active and reactive recovery time constant, s and s are loadvoltage steady-state dependency, and t and t are load-voltage transient dependency. Figure 4 shows exponential recovery load model response to a step-change in load voltage. Loads behavior characterized by its parameters, time constant p and q shows the time it takes for load to recover 63% of final load power from disturbance occurrence moment. Steady-state loadvoltage dependency s and s represents power that the load can finally recover, zero value in this parameter means total load recovery to pre-disturbance. Transient load-voltage dependency t and t show load response to voltage variation in the change moment, load model won't show any reaction in case of voltage change if this parameters considered equal to zero.  [20] State-equations of ERL model linearized in (46) and (47).
where v bDQ0 are pre-disturbance bus voltage values.

Complete state-space model of the microgrid
As can be seen in (31), (36), (37) and both load model statespace models, bus voltages are considered as input for load and DGs. Therefore, to ensure bus voltages are well defined, a virtual resistor (r N ) is assumed to be between each bus node and ground. This resistor is chosen large enough to make sure to have minimum influence on the stability of the system. Bus voltages can be achieved in terms of currents of DGs, network, and load as below: For each bus, these equations calculate by considering input or output currents from/to network and load current in load buses. In the static load model, i loadD and i loadQ are loads statevariables, on the other hand, in the dynamic load model currents can be obtained from (50) and (51): Above equations linearized as: Δi loadQ = a 5 ΔP L + a 6 ΔQ L + a 7 Δv bD + a 8 Δv bQ .
And coefficients a 1 , … , a 8 are: By linearizing (48) and (49), and using (52) and (53), stateequations of dynamic load model are expressed based on other state-variables. Finally, the complete microgrid small-signal state-space model and system state-matrix can be obtained by combining individual subsystem models. (60) ΔX INV and ΔX line are discussed, respectively, in Sections 2.1 and 2.2. In case of load, for the static model, load currents are chosen as state-variables ΔX load , and in the dynamic model, as discussed, state-variables consist of ΔP L and ΔQ L . By putting each state-equations in (60), A MG for both load models obtained.

Eigenvalues analysis
The eigenvalue analysis of a system is a common tool in control theory to determine the stability of the power systems. Eigenvalues of the system can be derived from the linearized state matrix and show the different modes. Former studies concluded that low-frequency eigenvalues near the imaginary axis are the most dominant and crucial for system stability, while the eigenvalues closeness to the real axis will affect the dynamic performance and damping properties of the system [16]. Participation factors (PFs) scale the efficiency of any system parameter on any system mode. It is important to know how any desirable changes of the system such as stabilizing unstable modes or responsibility improvement, can be achieved.
The PF of the i th state and j th eigenvalue is derived by left eigenvectors (w i j ) and right eigenvectors (e ji ) of the state matrix. p i j represent the participation factor of i th state on j th mode of the system and calculate from (61): For easy comparison between values, participation factors transform into a scale of 0 to 1. Thus the sum of each row is equal to 1 as it is the sum of all effects of that state on any mode. The higher participation factor of the i th state represents more influence on j th mode, so for control of each particular mode, the state with the highest participation factor will be used.
The flowchart of the proposed analysis method is illustrated in ten steps as presented in Figure 5. Note that, steps 9 and 10 which present the adjustment procedure of the controller parameters for stabilizing the system, in the case of microgrids' instability, will be further discussed in Section 4.

CASE STUDY AND RESULTS
A test case system with 220 V (per phase RMS), 50 Hz, similar to the one in [1] has been simulated to verify results and then replace the CIL model, used in paper, with new dynamic ERL model. As can be seen in Figure 6, test system consists of three similar inverter-based DGs connected to each other and load through the network. Here, network consists of two lines modeled with a simple RL model. Load in bus 1 has 5.8 kW initially and will be faced with a step change in active and reactive pow-   Table 1. As mentioned, m p and n q parameters are related to sharing powers between all three DGs and because of parameter similarity, power shares equally. The system has been simulated in SIMULINK/MATLAB to obtain an equilibrium point as the initial condition for eigenvalue analysis and in the same time check system performance under certain configurations.

CIL model
As mentioned, the CIL model will not show load dynamics. State-variables of this model are load's current, which is injected to the network and loads by DGs. Therefore, load-related coefficients in state matrix are generally dependent on DGs, just as network models. In order to investigate the stability of the test system in the presence of a 5.8 kW active load, the eigenvalue plot is shown in Figure 6. The eigenvalues are clustered based on the location and similar response with respect to changes in controller parameters. Figure 7 shows that with no disturbances in load, the dominant eigenvalues of the system have relatively good stability and damping properties. The slope of the droop characteristics, especially m p , is very low, indicating very slight frequency varia- tions in response to the active power variation. It is important to be noted since droop coefficients can determine the stability of the system in case of load disturbances. In the following, the disturbance scenario in load is applied. The purpose of these loads step changes, firstly, is to check the system stability under load disturbance and, secondly, to study the system ability to track the load increase and observe the individual behavior of the sources in the course of production matching consumption. In any case, the criterion of the designer's assessment of the success in control system design is the stability of the frequency and voltage, and, further, the absence of overshoot in voltages and frequencies, which will greatly affect the quality of the power. It should be noted here that the final values and voltage or frequency drop from the base values during load increase are not considered since these are specified by the droop characteristic.
As shown in Figure 8 load step change performed in t = 0.7 s and DG1 instantly react to supply new demand, active power shares equally between each DG and reactive power is mostly supplied by DG1. It can be seen that closer DG reacts faster to load change, in this case, DG1. It is worth noting that in all sim-  ulation results first ripple is related to initial loading. Afterward outputs will be converged to a fixed quantity. Figure 9 shows the voltage and frequency values of each DG. The voltage values were reduced according to the droop characteristic of the controller section of DGs. As mentioned before about powers, the accountability of a source close to the load is much faster, that the reaction speed of second and third sources is also inversely related to distance. Due to a bigger change in the reactive power of DG1, its voltage will also decrease more in the moment of the load step. Similarly, DG 2 and 3 also have changes that are proportional to reactive power change, although in the shortest time again all voltages return to an equal value dictated by the Q-V characteristics of the droop controller.

ERL model
In the following section, dynamic ERL model replaces the CIL model in bus 1. ERL model used is introduced in Section 2.3.2. Table 2 shows the dynamic load parameters estimated in [19].
Parameters are chosen such that the model shows similar behavior in steady-state and this article is just focused on that, neglecting the transient behavior of load can bring what consequences. As shown in Figure 10, the system, which was capable of response sever loading condition, is clearly unstable consider- ing the new load model that has an independent state (contrary CIL model). In the following, using the state-space model, the eigenvalue spectrum of the system is plotted in Figure 11. The eigenvalues are clustered based on the relative location and their behavior in response to changes in controller parameters. Referring to Figure 7 and compared to Figure 11, it can be seen that by replacing the dynamic model, some of the eigenvalues of the system are placed on the unstable region, which confirms the results of the non-linear simulation. As mentioned earlier, the objectives of this study are: firstly to stabilize and secondly to improve the system operating conditions, thus 1 as unstable modes, 2 as steady-state error modes and 3 as dominant modes need more attention.
The reason for this instability is the addition of load dynamics to the VSI's dynamics. In the modeling of static load, there are also two-axis current terms for load, but the dynamics of these states depend on the dynamics of input current from the network and ultimately inverter's dynamics. Using static load model and neglecting load dynamics, is a setback in achieving the exact insight of the system.

CONTROLLER ADJUSTMENT
Participation factors of the studied system are obtained and by applying necessary changes on the state variables with higher participation factors, takes the system into stable performance. In Table 3, the coordinates of the unstable modes and corresponding state variables with most participation are shown.
Regarding participation factor values in Table 3, it can be seen that changes in the PI controller parameters along with droop controller coefficients have the most effect on the unstable modes of the system. Thus, by applying 1 to 100 times decreasing changes in the parameters of K pv , K pc , K iv and K ic , their effect on the eigenvalues for system stabilization are investigated. By applying a 1-100 time decrease on the K pv parameter, unstable modes move toward the imaginary axis. The effect of these changes on the dominant stable modes is to move toward the imaginary axis and slightly far from the real axis which brings these modes more vulnerable to instability and more damping ratio. The desirable result is to firstly place modes on the negative side of the real axis and secondly, move far from the real axis to gain more damping.
The behavior of the modes by decreasing of K iv is: (i) unstable modes are fixed and no movement was seen, (ii) modes on origin point move in different directions some to unstable and some toward a stable region, and (iii) dominant modes of the system are not showing any significant movement.
Changing K pc results in stability of unstable modes, but on the other hand, it caused moving most stable modes toward the right side of the real axis. Thus first a balance between unstable modes stabilization and stable modes approaching the imaginary axis should be performed. Then using other parameters, a reasonable margin between modes and imaginary axis should be created.
Decreasing K ic has different effects on the location of the modes, as shown in Figure 12(d). Unstable modes respond in two ways to these changes, some of them show a decrease in imaginary part and place on the real axis and some move towards more instability. Dominant stable modes converge to the imaginary axis and originate from the axis. Aside from the problem of unstable modes, these changes resulted in less damping on dominant modes.
It was seen that a decrease of 1-100 times of each parameter will change the location of the modes; the purpose of changing system parameters is first to stabilize the system and second, to achieve the appropriate damping values, especially on the dominant modes. It was seen that system stability can be achieved by reducing the K pc and K pv parameters; on the other hand, by increasing the K iv and K ic parameters, the unstable modes will be more stable and dominant modes will gain more damping. Figure 13 is obtained by simultaneous changes on K pc , K pv , K iv and K ic as described in Table 4. Note that selected parameters are determined based on the changes in the eigenvalues location, as shown in Figure 12. Since the the presented work has been focused on the dynamic load effects on the microgrid stability, no objective function has been optimized to evaluate the value of stability improvement. Therefore, in the step of controller   parameters calculation, only appropriate performance improvement has been considered. Figure 13 shows that system modes are in the stable region and dominant modes have acceptable damping properties. Unstable modes converged to the origin of the axis which results in a steady-state error as a side effect of system stabilization. It can be solved by optimizing a cost function based on eigenvalue placement, but it requires a reduced-order model of the entire microgrid. Now it is time to check the performance of the system under new control coefficients by applying the same scenario applied to the static load. Results of 16.8 kW and 12 kVar step change in load are shown in Figures 14 and 15. As shown in Figure 14, all three DGs responding to load change, DG1 which in the static load test had overshot to 10 kW here is not so far from final values. As mentioned before, the model has a steady-state error but this error is much smaller than the static one, because of controller adjustment. As shown, voltage and frequency converged in an acceptable time and almost no

Time constant analysis
Introducing the ERL model, it was mentioned that the time constant changes control the convergence rate of the power response to the voltage step. Figure 16 shows the effect of 1, 10, 100, 200-time changes for the constant load power on the active power response. The time constant values in [19] chose to have the fastest response, or in the other words the worst possi-

FIGURE 17
Eigenvalue spectrum in response to load time constant variation ble condition for the system. Here, with time constant increase, the system has more time to converge to its final power value. Figure 17 shows the effect of the 200-time variation of the time constant on the dominant modes of the system. In the linearized equation of active power recovery (P r ), time constant parameter can be seen in the denominator of all the coefficients of the corresponding states; accordingly, with increasing time constant, the coefficients became smaller so loads related modes became less dependent to other state variables. The eigenvalue spectrum shows that with increasing active power time constant, a slight change in load-related modes can be seen. On the other hand, by this change, load-related modes will not be responsible for other parameter variations. Figure 18 shows the frequency and voltage response to the load step change. As a result of time constant increase, load change causes a drop in voltage and frequency based on droop characteristic. The voltage drop will effect load power according to ERL model equations, but in the case of small time constants, instead of step-change load power will respond slower which makes enough time for the system to match generation to the load demand.

CONCLUSION AND DISCUSSION
In this paper, a small-signal stability analysis of microgrid, based on the derived model in alternative study in [1] with static RLL model is investigated. Results from the disturbance scenario of the understudy system show that the microgrid is well responsive to load variation and is even favorable for severe loading based on both eigenvalues spectrum and Matlab/Simulink simulations. But when the load model is replaced by standard dynamic ERL model, system stability compromises even for loads with load dependencies at a voltage close to the constant impedance model ( s = 2). ERL model represents a wide range of loads such as heating loads which control the temperature with thermostatic function. The parameters of ERL model used in this study are determined by [19]. To be more cautious, in this paper the worst-case scenario due to the smallest time constant ( p , q ) is assumed. Consequently system operation would be unstable. Then the eigenvalues spectrum of the system with a dynamic load model reveal that some of the eigenvalues are in the unstable region. Afterward, a participation analysis based on right and left eigenvectors is performed. Based on its results, most participated states were found to be related to controller sections. Each controller parameter varied to find a stable setting based on the eigenvalues spectrum. Finally a set of controller parameters determined by Matlab/Simulink simulation, investigate the system operation with new controller setting. Increasing the load demand causes voltage and frequency drop which is based on droop characteristics, therefore they had no changes. On the other hand, low fluctuations are due to the proper operation of the controllers. The simulations show the efficiency of the proposed method for microgrid stabilization in the presence of dynamic loads. Dynamic stability of microgrid studies in the future can focus on suggestion addressed in the following: (i) investigating the effect of underground cable capacitors on the model and it's stability, (ii) adding a tuner part for controller section which tune controller coefficients after any alteration in system condition, and (iii) considering a specific objective function for improving the stability of the system in the step of adjusting the controller parameters.