A Novel, Stable, and Economic Power Sharing Scheme for an Autonomous Microgrid in the Energy Internet

Abstract

With a higher penetration of distributed generation in the power system, the application of microgrids is expected to increase dramatically in the future. This paper proposes a novel method to design optimal droop coefficients of dispatchable distributed energy resources for a microgrid in the Energy Internet considering the volatility of renewable energy generation, such as wind and photovoltaics. The uncertainties of renewable energy generation are modeled by a limited number of scenarios with high probabilities. In order to achieve stable and economical operation of a microgrid that is also suitable for plug-and-play distributed renewable energy and distributed energy storage devices, a multi-objective optimization model of droop coefficients compromising between operational cost and the integral of time-weighted absolute error criterion is developed. The optimization is solved by using a differential evolution algorithm. Case studies demonstrate that the economy and transient behavior of microgrids in the Energy Internet can both be improved significantly using the proposed method.

1. Introduction

Due to the environmental and geopolitical costs of fossil fuels, a federal program was created to stimulate essential innovation in energy technologies, such as wind and solar power [1,2,3,4]. In China, the cumulative installed capacity of wind power will reach 100 GW by 2015 and will easily meet the target of 200 GW by 2020 [5,6]. However, renewable energy sources have a fluctuating and unpredictable nature. Due to these characteristic uncertainties, the ever-increasing level of renewable energy sources in the power system makes system operation a challenging task and can introduce higher operational costs and stricter network constraints [7,8,9,10,11].

Recently, the concept of an “Energy Internet,” has been proposed by the Industrial Internet Consortium (IIC) and it is thought that microgrids in the Energy Internet can provide promising solutions [12]. The advantages of microgrids in the Energy Internet include supplying more reliable power to loads, integrating more renewable energy into power system, increasing energy utilization efficiency, and offering plug-and-play capabilities of distributed energy resources (DERs) and distributed energy storage devices.

To realize proportional power-sharing among distributed generations (DGs), some theoretical optimization methods and control schemes for Energy Internet microgrids have been proposed. Liu et al. [13] presented a model of multi-objective economic dispatch of the microgrid system, including Vehicle-to-Grid using the judgment matrix method for converting the multi-objective functions to a single objective function. In [14], Hawkes et al. developed a linear programming-based annual cost minimization model for unit commitment of generators and storage in a microgrid. A new distributed real-time electricity allocation scheme for the smart grid/microgrids energy management system (EMS) was proposed in [15]; its objectives are to reduce the electricity bills of residential customers, to increase the overall social benefit of the smart grid/microgrids community, and to increase the energy efficiency and reliability of the microgrids to rely on locally-generated electricity. Mohamed et al. [16] created a cost optimization scheme for a microgrid taking into account the costs of fuel, operation, maintenance, start-up, and emission. However, the implementation of the optimum economic solutions derived by [13,14,15,16] requires a centralized control system and additional communications infrastructure, which would result in a rise of investment cost and also limit the expansion of microgrids. The power sharing schemes addressed in these papers do not include any stability-related analysis and, thus, there is no guarantee the stability of microgrids could be maintained in the stated economic supporting operation conditions. In fact, the economic justifications may sometimes collide with the stability requirements.

Several control methods have been proposed to improve the stability of an autonomous microgrid. Nguyen et al. [17] proposed a new flywheel energy storage system based on a doubly-fed induction machine and a battery for use in microgrids that can control its power to compensate for fluctuations in wind power in grid-connected mode, in addition to controlling the frequency and voltage of the microgrid in islanded mode. Su et al. [3] proposed an adaptive robust total sliding-mode control system for the microgrid stability controller. It is proven that the adaptive, robust total sliding-mode control system is insensitive to parametric uncertainties and external disturbances. An improved dynamic model for the energy storage system in the context of grid-connected microgrids is proposed in [18]. Based on the model, a Proportional-Integral-Resonant DC link voltage controller is introduced to maintain the DC link voltage through the charging-discharging control of ultracapacitors and is capable of working properly under all operating conditions. However, these papers focus on advanced control methods to improve the stability of the microgrid and most of these rely on the centralized and communicated power-sharing schemes.

Compared to the centralized and communicated power sharing schemes mentioned above, a droop control is more desirable due to its reliability and robustness [19,20,21]. Under this circumstance, dispatchable distributed energy resources (DDERs) are controlled based on the droop characteristics to share the unbalanced power between outputs of non-dispatchable distributed energy resources (NDERs) and loads using local measured information. Inspired by the aforementioned consideration, this paper presents a multi-objective optimization model of droop coefficients compromising between the operational cost and the integral of time-weighted absolute error (ITAE) criterion. The detailed ITAE criterion is defined in [22] and has the advantage of producing smaller overshoots and oscillations. The model improves the dynamic behavior of the microgrid and is suitable for plug-and-play of distributed renewable energy and distributed energy storage devices. The proposed power sharing scheme can preserve not only the economic dispatch among DDERs, but also the stability of microgrids at all possible operating conditions. Due to the nonlinearity and complexity of the proposed optimization problem, the differential evolution (DE) algorithm is applied to solve the problem. The effectiveness of the proposed method has been demonstrated via numerical simulations on MATLAB R2013a/Simulink (The MathWorks, Natick, Massachusetts, MA, USA).

The original contribution of this paper is to provide an insight as to how to manage a decentralized droop control-based microgrid, from both an economic- and a stability-based viewpoint, via simply scheduling the droop coefficients. The remainder of this paper is organized as follows. In Section 2, description of the microgrid in the Energy Internet is presented. In Section 3, the multi-objective optimization model of droop coefficients considering the volatility of wind and solar power generation is proposed. Simulation results are provided to demonstrate the proposed method in Section 4. Conclusions are drawn in Section 5.

2. Microgrid in the Energy Internet Description

Figure 1 shows the structure of the microgrid in the Energy Internet used in this study. It is a cluster of distributed resource units and loads, serviced by a distribution system, and can operate in (1) grid-connected mode; (2) islanded (autonomous) mode; (3) ride-through between the two modes. Under normal circumstances, the microgrid shown in Figure 1 operates in autonomous mode. However, when it is subjected to large disturbances, such as loss of generating units or major lines, the system must operate in grid-connected mode. As a benefit of the droop control strategy, the plug-and-play of distributed renewable energy, distributed energy storage devices, and added loads can be easily realized. Therefore, there are two parts of the studied structure: the exited part and the extended part. In the exited part, three DDERs, a wind turbine generator (WTG), a photovoltaic plant (PV), an energy storage system (EES) and the commercial loads are installed at the main bus. For achieving a more refined load management, the commercial loads are divided into two groups: interrupted loads (IL) and non-interrupted loads (NL). The inverter controllers of WTG and PV employ the maximum power point tracking (MPPT) concept to output the captured wind and solar power to the system. All of the inverter controllers of DDERs employ the droop control concept to automatically undertake the unbalance power when the wind and solar power and the commercial loads change. The inverter controller of the EES employs V/f droop control to support the voltage and frequency of the system in the transient process through adjusting the output power by testing the change of terminal voltage and frequency. The extended part can include one or more DGs or other added commercial loads. The important note is that the microgrid must operate in grid-connected mode when a new unit has been interfaced into the system. Appendix A shows the detailed configuration of the microgrid.

Figure 1. Structure of a microgrid in the Energy Internet.

Figure 1. Structure of a microgrid in the Energy Internet.

3. The Economic and Stable Power Sharing Scheme

3.1. Modeling of Wind and Solar Power Uncertainties

3.1.1. Wind Power Forecast

Wind is a fluctuating power source. Therefore, it is necessary to predict the power output of wind generators for the next few hours to days ahead as accurately as possible. As wind power cannot be predicted with 100% accuracy, it is possible that the actual wind power could be different from the forecasted one. In this paper, we assume that the wind power forecast is available and the wind power forecast error is simulated based on the Auto-Regressive and Moving Average Model (1, 1) (ARMA) model, which was originally proposed in Wind Power Integration in Liberalised Electricity Markets (WILMAR) project [20,21,22,23,24].

It is assumed in the ARMA (1, 1) model that the wind power forecast error is subject to a normal distribution $N(0,{\mathrm{\text{σ}}}_t^2)$ in each time interval t. The variance series of wind power forecast error can be calculated as:

$$\begin{array}{l}{V}_0=0\hfill \\ V_1={\mathrm{\text{σ}}}_z^2\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}⋮\hfill \\ V_t={\mathrm{\text{μ}}}^2{V}_{t-1}+(1+{\nu }^2+2\mathrm{\text{μ}}\nu ){\mathrm{\text{σ}}}_z^2\hfill \end{array}$$

(1)

Correspondingly, the standard deviation series ${\sigma }_t$ is calculated as:

$${\mathrm{\text{σ}}}_t=\sqrt{V_t}$$

(2)

The parameters ${\mathrm{\text{σ}}}_Z$, μ, and $\nu$ in ARMA (1, 1) model can be identified by fitting the forecast error series as close as possible to available historical data.

In Figure 2, the time series of the forecasted wind power with its 95.4% probability interval ($\pm 2{\mathrm{\text{σ}}}_t$) are shown.

Figure 2. Forecasted wind power and its 95.4% probability interval.

Figure 2. Forecasted wind power and its 95.4% probability interval.

3.1.2. Scenario Generation

The wind power scenarios representing the future realization of uncertainties are generated using the Latin hypercube sampling (LHS) technique in this paper. LHS is a stratified sampling method which is very efficient at extracting samples from the entire distribution of random variables [25]. Compared to the traditional Monte Carlo Simulation which randomly samples the cumulative distributions to obtain variable inputs, LHS uniformly samples the cumulative distributions and approximates the required normal distribution much better [26]. The procedure for LHS can be divided into two main steps; sampling and permutation. The objective of sampling is to generate representative samples to cover probabilistic distributions of each input random variable (wind power at different time stage). The permutation aims at reducing the correlations between samples of different input random variables using the Cholesky decomposition method [27].

In our experiment, a scenario set with 3000 scenarios is generated. To alleviate the computational burden this scenario set must be reduced to contain only a finite number of scenarios which represent the original process as closely as possible.

3.1.3. Scenario Reduction

A scenario reduction technology is employed to determine a set with the requested number of scenarios and assign optimal probability to the preserved scenarios. The scenarios with very low probability will be eliminated and the close scenarios will be aggregated by measuring the distance of probability distribution between every pair of scenarios [28]. The initial scenario set with 3000 scenarios is ultimately reduced to 10 scenarios, as shown in Figure 3.

The modeling of solar power uncertainties is congruent to wind power and will not be repeated here. A scenario tree for stochastic solar power is shown in Figure 4 and the numerical results of scenario generation and reduction are shown in Appendix B.

Figure 3. Scenario tree for stochastic wind power.

Figure 3. Scenario tree for stochastic wind power.

Figure 4. Scenario tree for stochastic solar power.

Figure 4. Scenario tree for stochastic solar power.

3.2. Dynamic Power Sharing Scheme

The assignment of a droop coefficient to each DDER in the traditional method is usually based on an equitable load share in the form of m_p1P₁ = m_p2P₂ = … =m_pnP_n = constant, n_q1Q₁ = n_q2Q₂ = … = n_qnQ_n = constant, where P_i and Q_i are the rated active and reactive power output of DDER i [19]. It is remarkable that neither the economic nor the stability content is considered therein; the droop coefficients of the DDERs are scheduled based on rated power only. In fact, the operational cost and the stability of the microgrid are all closely related to the assignation of droop coefficients.

In this paper, a droop-based power sharing scheme uses the frequency and voltage of the microgrid as a common signal among the DDERs to balance the active and reactive power mismatch between generation and demand. The power injected from DDER i is mainly determined by its droop coefficient m_pi and n_qi, which has been proposed in the preliminary work published in [29]. It is shown in Appendix C.

3.3. Multi-Objective Optimization Model of Droop Coefficients

3.3.1. Pareto-Optimality Concept and Multi-Objective Optimization

A multi-objective optimization problem can be described using the following formulation:

$$\min \hspace{0.17em}y=F\left(x\right)=\left\{f_1\left(x\right),f_2\left(x\right),\mathrm{...},f_n\left(x\right)\right\}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in R^m,y\in R^n$$

(3)

while f_p(·) is the pth objective function (for 1 ≤ p ≤ n).

In most cases, there does not exist a global optimal solution to such a problem. However, the Pareto-optimality concept has been intensively used and solutions included in the Pareto-optimal set are those that cannot be improved along any dimension without simultaneously deteriorating along another dimension(s) [30].

3.3.2. Mathematical Formulation of the Proposed Problem

To minimize the total operational cost, the droop coefficient of each DDER should be selected carefully so that the DDERs with high and low fuel cost can be forced to share a relatively small and large portion of unbalanced power between generation and demand. On the other hand, the stability of a microgrid is also largely affected by the selection of droop coefficients. With inappropriate settings, a large overshoot or oscillation would appear in the transient process of the microgrid with the variation of renewable energy generation and load. In order to achieve both the economic dispatch among multiple DDERs in steady state and the satisfactory dynamic behavior of the microgrid in transient with the variation of operating conditions, a multi-objective optimization model of droop coefficients compromising between the operational cost and the ITAE criterion is proposed, which indicates the dynamic behavior of the microgrid. The multi-objective function is defined as:

$$\min \hspace{0.17em}C={\mathrm{\omega }}_1\times f_1+{\mathrm{\omega }}_2\times SC\times f_2$$

(4)

where f₁ is the operational cost and f₂ is the ITAE criterion. SC is the scaling factor, blending f₂ with f₁, and ω₁ and ω₂ are the weighting factor whose value varies between [0, 1].

The operational cost f₁ is formulated as:

$$f_1=[{\displaystyle \sum _{t=1}^{N_t}{\displaystyle \sum _{i=1}^{N_i}F{C}_i^t(P_{i,new}^t)}}\text{+}{\displaystyle \sum _{t=1}^{N_t}{\displaystyle \sum _{i=1}^{N_i}{\displaystyle \sum _{k=1}^{N_k}{\alpha }_{k}E{F}_{i,k}{P}_{i,new}^t}}}+C_{net}\times E_{net}]\hspace{0.17em}\hspace{0.17em}\forall t\in N_t,\hspace{0.17em}\forall i\in N_i,\hspace{0.17em}\forall k\in N_k$$

(5)

where $P_{i,new}^t$ is the steady state power output of DDER i at time stage t; $F{C}_i^t(P_{i,new}^t)$ is the fuel cost of DDER i at time stage t; α_k is the externality cost of emission type k (such as nitrogen oxides (NO_x), sulfur oxides (SO₂), carbon oxides (CO₂), and so on); EF_i,k is the emission factor of DDER i, emission type k; N_t, N_i, and N_k are the numbers of time stages, DDERs, and emission types, respectively; E_net is the planned exchange electric quantity between the microgrid and the main grid during the scheduling period and indicates that the microgrid purchases electricity from the main grid when the value is positive, and vice versa; and C_net is the purchase or sale of electricity price for microgrid, which can be descripted as piecewise function:

$$C_{net}=\{\begin{array}{c}{C}_b,\hspace{0.17em}{E}_{net}>0\hfill \\ C_s,\hspace{0.17em}{E}_{net}<0\hfill \end{array}\hspace{0.17em}$$

(6)

where C_b and C_s are the purchase or sale of electricity price for the microgrid, respectively.

The planned exchange electric quantity between the microgrid and the main grid during the scheduling period T can be calculated by Equation (7):

$$E_{net}={\displaystyle {\int }_0^T{P}_{net}}dt$$

(7)

where P_net is the constant exchange power between the microgrid and the main grid during the scheduling period T, which indicates that the power flows from main grid to the microgrid when the value is positive, and vice versa.

It can be seen in Equation (5) that f₁ is composed of fuel costs and environmental externality costs over the scheduled horizon. The operational costs of WTG and PV are assumed to be zero.

The fuel cost of DDER i at time stage t is modeled as a quadratic function of steady state power output of DDER i, as:

$$F{C}_i^t(P_{i,new}^t)=a_{Fi}+b_{Fi}{P}_{i,new}^t+c_{Fi}{(P_{i,new}^t)}^2\hspace{0.17em}\forall t\in N_t$$

(8)

where a_Fi, b_Fi, and c_Fi are the cost coefficients of DDER i.

The ITAE criterion f₂ is formulated as:

$$f_2={\displaystyle \sum _{k=k_0}^{T_0}\left(k-k_0\right)}\cdot W\cdot E_{abs}(k)$$

(9)

where k is the current sample time; k₀ is the starting time of wind generation and load change; T₀ is ending time of calculation; W is the weighting matrix which is an identity matrix in this paper; and E_abs(k) is the absolute error vector and it is defined as:

$$\begin{array}{cc}\hfill E_{abs}(k)=& [{\displaystyle \sum _{t=1}^{N_t}\left|f_{new}^t-f^t(k)\right|},{\displaystyle \sum _{t=1}^{N_t}{\displaystyle \sum _{j=1}^{N_j}\left|V_{j,new}^t-V_j^t(k)\right|}},\hfill \\ & {\displaystyle \sum _{t=1}^{N_t}{\displaystyle \sum _{i=1}^{N_i}\left|P_{i,new}^t-P_i^t(k)\right|}},{\displaystyle \sum _{t=1}^{N_t}{\displaystyle \sum _{i=1}^{N_i}\left|Q_{i,new}^t-Q_i^t(k)\right|}}{]}^T\hspace{0.17em}\forall t\in N_t\hfill \end{array}$$

(10)

where $f_{new}^t$ is the steady state frequency of the microgrid at time stage t; $V_{j,new}^t$ is the magnitude of steady state voltage of bus j; N_j is the number of buses; $Q_{i,new}^t$ is the steady state reactive power output of DDER i at time stage t; and $f^t(k)$, $V_j^t(k)$, $P_i^t(k)$ and $Q_i^t(k)$ are the measured frequency, voltage magnitude of bus j, power, and reactive output of DDER i at time stage t, respectively.

The objective function Equation (4) is subject to a set of operational constraints including the frequency and voltage regulation constraints and DDER capacity constraints. All these constraints are summarized as:

$$\begin{array}{l}{f}^{\text{min}}\le f^t(k)\le f^{\text{max}}\hfill \\ V_j^{\text{min}}\le V_j^t(k)\le V_j^{\text{max}}\hfill \\ \sqrt{{(P_i^t(k))}^2+{(Q_i^t(k))}^2}\le S_i^{\text{max}}\hfill \end{array}$$

(11)

where fmin and fmax are the allowed minimum and maximum frequency, $V_j^{\text{min}}$ and $V_j^{\text{max}}$ are the lower and upper limits of bus j voltage magnitude, and $S_i^{\text{max}}$ is the inverter capacity of DDER i.

The objective function Equation (4) is subjected to a set of system constraints, including the power balance constraints, generating constraints, EES constraints, operational reserve constraints, frequency, and voltage regulation constraints and scenario constraints. All the constraints are summarized as follows.

(1) Power balance constraints

$${\displaystyle \sum _{i=1}^{N_i}{P}_i^t}+P_W^t+P_P^t+P_E^t+P_{net}=P_{NL}^t+P_{IL}^t\hspace{0.17em}\hspace{0.17em}\forall t\in N_t,\hspace{0.17em}\forall i\in N_i$$

(12)

where $P_W^t$ and $P_P^t$ are the forecasted wind and solar power generation at time t; $P_E^t$ is the exchange power of EES at time t; $P_{NL}^t$ and $P_{IL}^t$ are the non-interrupted and interrupted loads at time t. Notice that the network loss is not considered here.

(2) Generation constraints

$$P_{i,\text{min}}\le P_i^t\le P_{i,\text{max}}\hspace{0.17em}\hspace{0.17em}\forall i\in N_i,\forall t\in N_t$$

(13)

where $P_{i,\text{min}}$ and $P_{i,\text{max}}$ are the lower and upper operating limits on generation of DDER i, respectively.

(3) EES constraints

The charge or discharge electric quantity of EES during period T is as:

$$E_E^t={\displaystyle {\int }_0^T{P}_E^t}dt\hspace{0.17em}\hspace{0.17em}\forall t\in N_t$$

(14)

The residual electric quantity of EES at the end of period T is as:

$$E_r^t=E_{int}+{\displaystyle \sum _{t=1}^T{E}_E(T)}\hspace{0.17em}\hspace{0.17em}\forall t\in N_t$$

(15)

where $E_{int}$ is the initial electric quantity of EES.

The electric quantity and power constraints of EES are as:

$$E_{E,\min }\le E_r^t\le E_{E,\max }\hspace{0.17em}\hspace{0.17em}\forall t\in N_t$$

(16)

$$-P_{E,min}\le P_{ESS}^t\le P_{E,\max }\hspace{0.17em}\hspace{0.17em}\forall t\in N_t$$

(17)

where $E_{E,\min }$, $E_{E,\max }$, $P_{E,min}$ and $P_{E,\max }$ are the lower and upper operating limits of the electric quantity and exchange power of EES, respectively.

(4) Operational reserve constraints:

In order to ensure the non-interrupted loads maintain continuous operation, all the power output of the microsources must be greater than the power demand of non-interrupted loads, considering the microgrid is likely to disconnect from the main grid at any time.

$${\displaystyle \sum _{i=1}^{N_i}{P}_{i,\text{max}}}+P_{E,\max }+P_W^t+P_P^t\ge P_{NL}(k)\hspace{0.17em}\hspace{0.17em}\forall t\in N_t,\hspace{0.17em}\forall i\in N_i$$

(18)

Meanwhile, in order to ensure the non-interrupted loads keep continuous operation for a period T at the end of time t, the whole electric quantity of the microgrid must be greater than the electric quantity demand of non-interrupted loads, as the energy storage constraints of EES.

$$E_r^t+{\displaystyle \sum _{i=1}^{N_i}{\displaystyle {\int }_0^T{P}_{i,max}}dt}+{\displaystyle {\int }_0^T{P}_P^t}dt+{\displaystyle {\int }_0^T{P}_W^t}dt\ge {\displaystyle {\int }_0^T{P}_{NL}^{t}dt}\hspace{0.17em}\hspace{0.17em}\forall t\in N_t,\hspace{0.17em}\forall i\in N_i$$

(19)

(5) Frequency and voltage regulation constraints

$$f^{\min }\le f^t\le f^{\max }\hspace{0.17em}\hspace{0.17em}\forall t\in N_t$$

(20)

$$V^{\min }\le V^t\le V^{\max }\hspace{0.17em}\hspace{0.17em}\forall t\in N_t$$

(21)

where $f^{\min }$, $f^{\max }$, $V^{\min }$ and $V^{\max }$ are the allowed minimum and maximum frequency and voltage, respectively, which are chosen as 49 Hz, 51 Hz, 0.95 p.u. and 1.05 p.u., respectively.

(6) Scenario constraints

For the sake of considering uncertainties of wind and solar power generation, the scenario constraints have been added to the model for the system power balance.

$$\{\begin{array}{cc}\begin{array}{l}{\displaystyle \sum _{i=1}^{N_i}{P}_i^{t,s}}+P_W^{t,s}+P_P^{t,s}+P_E^{t,s}+P_{net}=P_{NL}^t+P_{IL}^t\hfill \\ P_{i,\text{min}}\le P_i^{t,s}\le P_{i,\text{max}}\hfill \\ E_{E,\min }\le E_r^{t,s}\le E_{E,\max }\hfill \\ -P_{E,min}\le P_{ESS}^{t,s}\le P_{E,\max }\hfill \\ {\displaystyle \sum _{i=1}^{N_i}{P}_{i,\text{max}}}+P_{E,\max }+P_W^{t,s}+P_P^{t,s}\ge P_{NL}(k)\hfill \\ E_r^{t,s}+{\displaystyle \sum _{i=1}^{N_i}{\displaystyle {\int }_0^T{P}_{i,\text{max}}}dt}+{\displaystyle {\int }_0^T{P}_P^{t,s}}dt+{\displaystyle {\int }_0^T{P}_W^{t,s}}dt\ge {\displaystyle {\int }_0^T{P}_{NL}^{t,s}dt}\hfill \\ f^{\min }\le f^{t,s}\le f^{\max }\hfill \\ V^{\min }\le V^{t,s}\le V^{\max }\hfill \end{array}\hfill & \hspace{0.17em}\hspace{0.17em}\forall t\in N_t,\hspace{0.17em}\forall i\in N_i,\hspace{0.17em}\forall s\in N_s\hfill \end{array}$$

(22)

3.3.3. Fuzzy Logic Theory Involving to the Differential Evolution (DE) Algorithm

The proposed optimization problem shown above is a constrained nonlinear optimization problem. DE is used because of its characteristic of being derivative-free and, thus, effective in capturing the optimal results. It is to be mentioned that in the optimization process, after calculating the objective function, all the constraints are checked and a penalty term is added if any constraint is violated. The algorithm stops if there is no improvement in the objective function.

The values of objective functions can usually belong to different intervals with variable lengths and some objective functions may be always dominated by the others. For example, if $f_i\in \left[{10}^2,\hspace{0.17em}{10}^3\right]$ and $f_j\in \left[0,\hspace{0.17em}1\right]$ $\left(i\ne j,\hspace{0.17em}1\le i,\hspace{0.17em}j\le n\right)$, it is easy to find that the objective function $f_i$ is in a dominant position. In order to make the evaluation more efficient, fuzzy logic theory is used in the normalization of each objective function based on the following steps:

(1) Divide the interval of ω₁ and ω₂ into 10 subintervals and implement the optimization repeat from 0 to 1 with a step increment of 0.1. Two groups of non-dominated objective function f_1,i and f_2,i can be obtained and are presented in Table 1.

Table 1. Divisions of w and the objective function.

Variables	Scenarios
i	1	2	3	4	5	6	7	8	9	10	11
ω1	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
ω2	1	0.9	0.8	0.7	0.6	0.5	0.4	0.3	0.2	0.1	0
f1,i	f1,1	f1,2	f1,3	f1,4	f1,5	f1,6	f1,7	f1,8	f1,9	f1,10	f1,11
f2,i	f2,1	f2,2	f2,3	f2,4	f2,5	f2,6	f2,7	f2,8	f2,9	f2,10	f2,11

Table 1. Divisions of w and the objective function.

Variables	Scenarios
i	1	2	3	4	5	6	7	8	9	10	11
ω1	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
ω2	1	0.9	0.8	0.7	0.6	0.5	0.4	0.3	0.2	0.1	0
f1,i	f1,1	f1,2	f1,3	f1,4	f1,5	f1,6	f1,7	f1,8	f1,9	f1,10	f1,11
f2,i	f2,1	f2,2	f2,3	f2,4	f2,5	f2,6	f2,7	f2,8	f2,9	f2,10	f2,11

(2) Calculate the lower-bounds f_1,min, f_2,min and upper-bounds f_1,max, f_2,max of each objective function based on [30].

(3) For the two groups of non-dominated objective function f_1,i and f_2,i, a fuzzy membership function can be assigned as:

$$\begin{array}{l}\mu \left(f_{1,i}\right)=\{\begin{array}{c}\hfill 1,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{f}_{1,i}\le f_{1,\min }\hfill \\ \hfill \frac{f_{1,\max }-f_{1,i}}{f_{1,\max }-f_{1,\min }},f_{1,\min }\le f_{1,i}\le f_{1,\max }\hfill \\ \hfill 0,f_{1,i}\ge f_{1,\max }\hfill \end{array}\hfill \\ \mu \left(f_{2,i}\right)=\{\begin{array}{c}\hfill 1,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{f}_{2,i}\le f_{2,\min }\hfill \\ \hfill \frac{f_{2,\max }-f_{2,i}}{f_{2,\max }-f_{2,\min }},f_{2,\min }\le f_{2,i}\le f_{2,\max }\hfill \\ \hfill 0,f_{2,i}\ge f_{2,\max }\hfill \end{array}\hfill \end{array}$$

(23)

(4) By choosing an appropriate ω, we can now reduce the multi-objective optimization to the minimization of the following fitness function:

$$Fitnes{s}_i={\mathrm{\omega }}_1\cdot \mathrm{\text{μ}}(f_{1,i})+{\mathrm{\omega }}_2\cdot \mathrm{\text{μ}}(f_{2,i})$$

(24)

(5) The best compromise solution can be found when the normalized sum of the membership function Equation (23) is highest:

$$C_i=\frac{Fitnes{s}_i}{{\displaystyle \sum _{i=1}^{11}(Fitnes{s}_i)}}$$

(25)

4. Experimental Verification of Proposed Method

The exited part of the microgrid system depicted in Figure 1 is used to demonstrate the proposed approach for finding the optimal droop coefficients. The rated frequency and voltage of the microgrid is 50 Hz and 380 V, respectively. The inverter capacity of each DDER is 300 kVA. The passive filter of each microsource is r_f = 0.1 Ω, L_f = 1.35 mH, C_f = 50 μF; the coupling resistance and inductance of each microsource is r_c = 0.01 Ω, L_c = 0.02 mH. The fuel cost coefficients, output power limits and externality cost coefficients of three DDERs are listed in Table 2 and Table 3. The capacity and other coefficients of the EES are listed in Table 4.

The total load demand of the microgrid is shown in Table 5, with 50% non-interruptible load.

This paper solves the optimization model presented in Equations (1)–(21) using DE and determines the best compromise solution according to the fuzzy logic theory. The value of the objective function Equation (2) is 0.0375 and the corresponding optimal droop coefficients of DDERs are shown in Table 6. The convergence characteristic of the objective function using the DE algorithm is shown in Figure 5. It can be seen that the convergence of the objective function is fast and the DE algorithm is capable of finding the optimum of the proposed nonlinear optimization problem.

Table 2. Fuel cost coefficients of dispatchable distributed energy resources (DDERs).

Fuel Cost Coefficients	DDER 1	DDER 2	DDER 3
aFi ($)	1.150	5.717	10.193
bFi ($/kWh)	57.783	100.937	105.180
cFi ($/kW2h)	−133.0915	0	62.560
Pi,min(kW)	10	15	20
Pi,max(kW)	200	300	400

Table 2. Fuel cost coefficients of dispatchable distributed energy resources (DDERs).

Fuel Cost Coefficients	DDER 1	DDER 2	DDER 3
aFi ($)	1.150	5.717	10.193
bFi ($/kWh)	57.783	100.937	105.180
cFi ($/kW2h)	−133.0915	0	62.560
Pi,min(kW)	10	15	20
Pi,max(kW)	200	300	400

Table 3. Externality costs and emission factors for NO_x, SO₂ and CO₂.

Emission Type k	αk ($/lb)	EF1,k (lb/kWh)	EF2,k (lb/kWh)
NOx	4.200	0.440	0.400
SO2	0.990	0.008	0.006
CO2	0.014	1.596	1.378

Table 3. Externality costs and emission factors for NO_x, SO₂ and CO₂.

Emission Type k	αk ($/lb)	EF1,k (lb/kWh)	EF2,k (lb/kWh)
NOx	4.200	0.440	0.400
SO2	0.990	0.008	0.006
CO2	0.014	1.596	1.378

Table 4. The coefficients of energy storage system (EES).

$P_{ESS}^{max}$ (kW)	$E_{ini}$ (kWh)	$E_{ESS}^{max}$ (kWh)	$E_{ESS}^{min}$ (kWh)
100	10	200	10

Table 4. The coefficients of energy storage system (EES).

$P_{ESS}^{max}$ (kW)	$E_{ini}$ (kWh)	$E_{ESS}^{max}$ (kWh)	$E_{ESS}^{min}$ (kWh)
100	10	200	10

Table 5. Total load demand.

Time Stage	Total Active Load (kW)	Total Reactive Load (kVar)	Time Stage	Total Active Load (kW)	Total Reactive Load (kVar)
1	535	175.8	13	789	259.3
2	575	189.0	14	783	257.4
3	574	188.6	15	790	259.7
4	579	190.3	16	696	228.8
5	585	192.3	17	630	207.1
6	598	196.6	18	599	196.9
7	667	219.2	19	621	204.1
8	707	232.4	20	577	189.7
9	689	226.5	21	591	194.3
10	683	224.5	22	687	225.8
11	796	261.6	23	688	226.1
12	794	261.0	24	606	199.2

Table 5. Total load demand.

Time Stage	Total Active Load (kW)	Total Reactive Load (kVar)	Time Stage	Total Active Load (kW)	Total Reactive Load (kVar)
1	535	175.8	13	789	259.3
2	575	189.0	14	783	257.4
3	574	188.6	15	790	259.7
4	579	190.3	16	696	228.8
5	585	192.3	17	630	207.1
6	598	196.6	18	599	196.9
7	667	219.2	19	621	204.1
8	707	232.4	20	577	189.7
9	689	226.5	21	591	194.3
10	683	224.5	22	687	225.8
11	796	261.6	23	688	226.1
12	794	261.0	24	606	199.2

Table 6. The optimal droop coefficients.

mp1 (rad/(s·W))	mp2 (rad/(s·W))	mp3 (rad/(s·W))	nq1 (V/Var)	nq2 (V/Var)	nq3 (V/Var)
2.596 × 10−6	3.845 × 10−6	4.811 × 10−6	8.152 × 10−7	5.493 × 10−6	9.885 × 10−6

Table 6. The optimal droop coefficients.

mp1 (rad/(s·W))	mp2 (rad/(s·W))	mp3 (rad/(s·W))	nq1 (V/Var)	nq2 (V/Var)	nq3 (V/Var)
2.596 × 10−6	3.845 × 10−6	4.811 × 10−6	8.152 × 10−7	5.493 × 10−6	9.885 × 10−6

Figure 5. Convergence curve of objective function.

Figure 5. Convergence curve of objective function.

As shown in Table 7, under the control of these optimal droop coefficients, the system operation cost f₁ is $1483.3 and the ITAE criterion f₂ is 0.0229, respectively. If considering only the minimization of operational cost, which is ω₁ = 1 and is utilized in Equation (2), f₁ will be $1276.60, which is lower than $1483.3. However, f₂ will be 0.3452, which is much larger than 0.0229. The system will have a poor dynamic response characteristic. On the other hand, if only considering the minimization of ITAE criterion, which is ω₁ = 0 and is utilized in Equation (2), f₂ will be 0.0124, which is lower than 0.0229. Here f₁ will be $1846.10, which is much larger than $1483.3. The system will operate in an uneconomical power-sharing scheme. The above comparison further validates the necessity to take into account both economic and stability considerations in the design process of optimal droop coefficients.

Table 7. The optimal results of the power sharing scheme.

Weighing Factor ω1	Weighing Factor ω2	Operation Cost f1 ($)	ITAE Criterion f2 ($)
0.5	0.5	1483.3	0.0229
1	0	1276.6	0.3452
0	1	1846.1	0.0124

Table 7. The optimal results of the power sharing scheme.

Weighing Factor ω1	Weighing Factor ω2	Operation Cost f1 ($)	ITAE Criterion f2 ($)
0.5	0.5	1483.3	0.0229
1	0	1276.6	0.3452
0	1	1846.1	0.0124

To further demonstrate that the dynamic power sharing scheme will enhance the transient and dynamic stability of the system, a comparison of numerical simulations between operation curves is presented below. Figure 6, Figure 7, Figure 8 and Figure 9 shows the responses of system frequency, voltage amplitude of the main bus, active power of DDER2, and reactive power of DDER1 with the optimal droop coefficients (ω₁ = 0.5, ω₂ = 0.5) and the economy-based optimal droop coefficients (ω₁ = 1, ω₂ = 0). As shown in

Table A1. Wind power generation in 10 scenarios and the probability of each scenario.

Time Stage	Wind Power Generation in 10 Scenarios (kW) and the Probability of Each Scenario
	1	2	3	4	5	6	7	8	9	10
1	132	132	132	132	132	132	132	132	132	132
2	212	216	205	224	207	211	204	213	204	217
3	233	235	225	230	210	228	210	232	222	230
4	241	258	241	221	229	246	256	275	237	244
5	270	251	250	253	256	252	268	265	253	265
6	226	254	261	258	250	252	229	262	255	241
7	303	295	292	253	344	300	280	307	304	330
8	307	335	283	282	292	300	295	302	350	299
9	263	209	244	221	246	234	231	232	236	250
10	186	215	207	204	186	192	216	189	174	205
11	317	283	291	305	302	300	292	338	256	296
12	243	269	289	243	270	276	274	299	266	312
13	228	241	263	251	269	252	241	264	251	261
14	200	236	259	267	212	240	230	227	241	251
15	250	223	221	237	223	234	196	246	226	237
16	93	86	94	94	88	96	93	107	97	92
17	42	40	40	41	39	42	37	42	49	43
18	55	53	48	60	43	54	61	53	52	44
19	56	60	66	55	63	60	66	52	61	56
20	45	44	45	43	41	45	44	50	40	46
21	57	49	44	42	55	48	48	50	44	55
22	144	170	146	177	174	168	173	195	167	161
23	235	242	252	256	259	246	270	216	255	210
24	148	124	130	166	156	156	147	144	134	142
Probability	0.078	0.151	0.067	0.108	0.048	0.188	0.113	0.105	0.090	0.052

, scenario 6 has the largest occurrence probability for wind power and is used as an example. To mimic the variation of operating conditions in a day, the wind generation and load change every 10 s according to Table 5, Table A1 and

Table A2. Solar power generation in 10 scenarios and the probability of each scenario.

Time Stage	Solar Power Generation in 10 Scenarios (kW) and the Probability of Each Scenario
	1	2	3	4	5	6	7	8	9	10
1	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
2	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
3	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
4	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
5	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
6	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
7	10.4	10.6	9.9	10.7	10.3	10.2	9.6	11.6	10.5	10.9
8	19.7	19.3	20.0	20.6	17.6	20.7	19.9	21.4	18.7	18.9
9	19.4	19.2	18.5	19.1	22.1	20.3	19.5	19.9	18.0	18.7
10	38.0	39.5	37.9	34.7	38.4	37.0	40.4	38.7	37.5	42.3
11	75.9	72.6	74.3	79.2	76.7	78.7	73.5	77.5	85.4	75.4
12	97.0	90.5	94.9	99.2	93.4	100.1	96.7	104.8	98.0	91.7
13	98.5	90.3	104.9	94.2	94.5	99.2	101.3	98.2	103.9	96.8
14	77.2	76.4	85.6	79.0	72.8	76.2	80.5	77.8	79.8	66.6
15	100	102.1	90.8	105.9	96.8	97.8	102.9	94.6	99.6	100.4
16	78.0	86.2	79.2	73.9	78.0	80.2	73.0	81.6	76.8	75.7
17	48.0	46.6	47.5	47.3	46.0	48.5	49.2	55.7	44.1	50.0
18	18.8	18.6	18.7	18.9	19.4	19.1	19.8	18.0	21.8	17.4
19	8.1	8.9	7.5	7.6	8.4	8.0	8.1	8.2	7.8	8.4
20	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
21	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
22	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
23	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
24	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
Probability	0.063	0.133	0.136	0.198	0.107	0.061	0.053	0.082	0.06	0.107

.

Figure 6. (a) The system frequency response with the optimal droop coefficients (ω₁ = 0.5, ω₂ = 0.5); (b) The system frequency response with the economy-based optimal droop coefficients (ω₁ = 1, ω₂ = 0).

Figure 6. (a) The system frequency response with the optimal droop coefficients (ω₁ = 0.5, ω₂ = 0.5); (b) The system frequency response with the economy-based optimal droop coefficients (ω₁ = 1, ω₂ = 0).

Figure 7. (a) The voltage amplitude response of bus3 with the optimal droop coefficients (ω₁ = 0.5, ω₂ = 0.5); (b) The voltage amplitude response of bus3 with the economy-based optimal droop coefficients (ω₁ = 1, ω₂ = 0).

Figure 7. (a) The voltage amplitude response of bus3 with the optimal droop coefficients (ω₁ = 0.5, ω₂ = 0.5); (b) The voltage amplitude response of bus3 with the economy-based optimal droop coefficients (ω₁ = 1, ω₂ = 0).

Figure 8. (a) The response of DDER2 output active power with the optimal droop coefficients (ω₁ = 0.5, ω₂ = 0.5); (b) The response of DDER2 output active power with the economy-based optimal droop coefficients (ω₁ = 1, ω₂ = 0).

Figure 8. (a) The response of DDER2 output active power with the optimal droop coefficients (ω₁ = 0.5, ω₂ = 0.5); (b) The response of DDER2 output active power with the economy-based optimal droop coefficients (ω₁ = 1, ω₂ = 0).

Figure 9. (a) The response of DDER1 output reactive power with the optimal droop coefficients (ω₁ = 0.5, ω₂ = 0.5); (b) The response of DDER1 output reactive power with the economy-based optimal droop coefficients (ω₁ = 1, ω₂ = 0).

Figure 9. (a) The response of DDER1 output reactive power with the optimal droop coefficients (ω₁ = 0.5, ω₂ = 0.5); (b) The response of DDER1 output reactive power with the economy-based optimal droop coefficients (ω₁ = 1, ω₂ = 0).

It can be observed in group (a) of Figure 6, Figure 7, Figure 8 and Figure 9 that the microgrid has excellent dynamic performance with the optimal droop coefficients. When the wind generation and load change, the transient process exhibits a small overshoot, little regulation time, and rapid oscillation decay. All the operational constraints are maintained in the allowable range (the maximum frequency deviation is −0.16 Hz, the maximum voltage amplitude deviation is +0.5 V, and the maximum apparent power of DDERs is 241.1 kVA).

If only considering the minimization of operation cost, which is ω₁ = 1, ω₂ = 0, and is utilized in Equation (4), the system would probably have very poor dynamic performance with the variation of operating conditions. The system response with the economy-based optimal droop coefficients is shown in group (b) of Figure 6, Figure 7, Figure 8 and Figure 9. It can be observed that the transient process exhibits very large oscillations and the system will lose stability at time stage 19, 20 and 21.

5. Conclusions

This paper proposes a method for scheduling droop coefficients to simultaneously improve the economy and stability of a microgrid, considering the volatility of wind power generation. A multi-objective optimization model of droop coefficients compromising between the operational cost and the ITAE criterion is presented. The simulation result validates the effectiveness of the proposed method. It is emphasized in this paper that both the economy and the stability of microgrids should be taken into consideration in the design process of droop coefficients. If not, the autonomous microgrid may suffer economic hardship or stability degradation with the variation of wind generation or load demand.