A Chance-Constrained Economic Dispatch Model in Wind-Thermal-Energy Storage System

Abstract

As a type of renewable energy, wind energy is integrated into the power system with more and more penetration levels. It is challenging for the power system operators (PSOs) to cope with the uncertainty and variation of the wind power and its forecasts. A chance-constrained economic dispatch (ED) model for the wind-thermal-energy storage system (WTESS) is developed in this paper. An optimization model with the wind power and the energy storage system (ESS) is first established with the consideration of both the economic benefits of the system and less wind curtailments. The original wind power generation is processed by the ESS to obtain the final wind power output generation (FWPG). A Gaussian mixture model (GMM) distribution is adopted to characterize the probabilistic and cumulative distribution functions with an analytical expression. Then, a chance-constrained ED model integrated by the wind-energy storage system (W-ESS) is developed by considering both the overestimation costs and the underestimation costs of the system and solved by the sequential linear programming method. Numerical simulation results using the wind power data in four wind farms are performed on the developed ED model with the IEEE 30-bus system. It is verified that the developed ED model is effective to integrate the uncertain and variable wind power. The GMM distribution could accurately fit the actual distribution of the final wind power output, and the ESS could help effectively decrease the operation costs.

1. Introduction

The uncertainty and variability of the integration of wind power are presenting more and more challenges to the power system operations [1,2]. It is difficult for the power system to retain secure, reliable and economic operations, especially at a high level of wind power penetrations. According to the actual situation of wind power integration, power system operators (PSOs) find that the large-scale wind power can increase or decrease in a very short time period with large and rapid fluctuations [3,4]. With the recent rapid development of the energy storage system (ESS), it is becoming a hot focus to use the ESS to mitigate the uncertainty and variability of wind power for government policies, power and energy industries and academic research. The ESS has been demonstrated to play a significant role in mitigating the fluctuations for both distributed power systems and the power system operations. More and more research is focusing on blending the wind power generators with an ESS to obtain smoother wind power generation. This research background makes it more practical to study the behavior of the wind power generation based on the statistical approaches. The statistical analysis of wind power generation could help PSOs better understand the uncertainties of wind power [5] and avoid the instability of the power system frequency [6,7], load shedding [8,9] and other operation problems [10,11].

1.1. Overview of the Economic Dispatch Model

With the increasing integration of wind power, the conventional economic dispatch (ED) model is facing the challenge of the uncertainty and variation of wind power forecasts. Academic researchers have developed numerous mathematical models to manage the wind power uncertainties in the ED problems. Wang et al. [12] formulated a bi-objective economic dispatch problem considering wind penetration and utilized a modified multi-objective particle swarm optimization algorithm to solve the ED model. Liu et al. [13] implemented Benders’ decomposition method to solve the security-constrained economic dispatch considering multiple stages of rescheduling to meet different security constraints. Niknam et al. [14] proposed a teaching-learning optimization algorithm to search the global optimal solution for the reserve constrained dynamic economic dispatch. Cui et al. [15] developed a multi-timescale power system operation model integrating both the unit commitment and economic dispatch sub-models. However, none of the above research consider the stochastic characteristics of wind power forecasts.

1.2. Overview of the Energy Storage System

The energy storage system has been widely used to solve the over-estimation or under-estimation in wind power integration [16]. The ESS is capable of improving the operation flexibility of the power system and mitigating the significant and rapid fluctuations of wind power generation. It has the ability to provide flexible charging and discharging. The application of energy storage technologies is a viable solution for the modern power system operations based on the recent development and advances in the ESS. It could help the wind power participate in the electricity market and provide operation services and reserves [17] with less wind power curtailments. Gong et al. [18] used the energy storage system to control wind power ramps by establishing an optimization model considering both the output power of ESS and the wind power curtailments. Jiang et al. [19] proposed a two-time-scale coordination control method to mitigate wind power fluctuations by using the ESS. A flexible first-order low-pass filter was developed to limit the wind power fluctuations and ramping rate. Wen et al. [20] proposed a hybrid multi-objective particle swarm optimization approach to minimize the power system cost and improve the system voltage profiles by searching the economic allocation and size of the ESS.

1.3. Overview of the Gaussian Mixture Model

To depict the uncertainty and variation of wind power forecasts, some statistical models have been widely used to model the probabilistic and cumulative distributions of the wind power forecasting data. The Gaussian mixture model (GMM) is one of the statistical models used to characterize the multi-modal and irregular distributions. It has been widely utilized in the renewable energy community. Ke et al. [21] customized the GMM by three Gaussian functions and utilized the GMM to approximate the probability density function (PDF) of wind power generation with triple probability peaks. An application of the GMM was performed in the probabilistic optimal power flow model in this paper. Singh et al. [22] represented all irregular PDFs of load using GMM in various distribution system applications and utilized the GMM distribution on a 95-bus generic distribution network model. Valverde et al. [23] proposed the use of GMM to represent non-Gaussian correlated wind power output and aggregated load demands to model the probabilistic load flow. The GMM distribution was testified on the 14-bus and 57-bus systems with deterministic power flow. However, none of the aforementioned papers adopt the GMM distribution into the economic dispatch model. In addition, the cumulative distribution of the GMM is used in the numerical manner rather than in the analytical manner, which limits the utilization of the GMM distribution.

1.4. Research Motivation and Objectives

The level of uncertainty and variability in power systems tends to increase with additional renewable energy [24]. The ESS is helpful to mitigate the uncertainty and variability of the original wind power generation and obtain the more smooth final output generation. However, the quantified impacts of the ESS on the power system operations are still unclear for the power system operators. The statistical characteristics of the final output of the wind power generation with the ESS (W-ESS) are also rarely studied by academic researchers. Some questions have not been addressed yet. For instance, how will the final output of the W-ESS be analytically expressed in the PDF and CDF? What is the impact of the ESS on the economic dispatch model integrated with wind power?

This paper aims to: (i) develop a chance-constrained ED model considering the W-ESS and the analytical distribution functions of the final wind power generation; (ii) compare the impact of different distribution models on the scheduled wind power generation; and (iii) quantitatively calculate the impact of the ESS on the chance-constrained ED model.

The remainder of this paper is organized as follows: The probabilistic distribution model of the final wind power generation based on the GMM distribution is presented in Section 2. The wind-energy storage system (W-ESS) considering both the usage of the ESS and wind power curtailments is provided in Section 3.1. The developed chance-constrained ED model integrated by the wind-energy storage system (W-ESS) is discussed in Section 3.2. The sequential linear programming method is introduced in Section 3.3 to solve the chance-constrained ED model. Section 4 presents the simulation results of a case study with the actual wind power data. Conclusions and a discussion of future work are given in Section 5.

2. Probabilistic Distribution Model of Final Wind Power Generation

Gaussian mixture models (GMM) are composed of ${\mathrm{N}}_{\mathrm{G}}$ multivariate normal density components, where ${\mathrm{N}}_{\mathrm{G}}$ is a positive integer. Each component has a d-dimensional mean (d is a positive integer), a d-by-d covariance matrix and a mixing proportion. Mixing proportion i determines the proportion of the population composed by component i, i = 1, ..., ${\mathrm{N}}_{\mathrm{G}}$. In this section, the nonlinear least square (NLS) method is used to fit GMMs to the actual data. There are two steps to obtain the parameters of the GMM distribution. First, the weight, mean value and standard deviation parameters are estimated by the NLS method using the trust-region algorithm. Second, the optimal number of mixture components is determined by minimizing the Euclidean distance between the actual PDF and the GMM PDF.

The developed probability distribution function of the GMM in this paper is a weighted sum of a finite number of Gaussian functions. It is characterized by the number of mixture components, standard deviations, mean values and weights of each component, given by:

$$f\left(x\right|{\mathrm{N}}_{\mathrm{G}};\Gamma )=\sum _{i=1}^{{\mathrm{N}}_{\mathrm{G}}}{\omega }_i{g}_i\left(x\right|{\mu }_i,{\sigma }_i),\forall x\in \mathcal{X},\forall i\in \mathcal{I}$$

(1)

where $\mathcal{X}$ is the dataset of a random variable, x, with a total number of ${\mathrm{N}}_x$. $\mathcal{I}$ is the set of mixture components with a total number of ${\mathrm{N}}_{\mathrm{G}}$. Γ is the overall parameter matrix. Each vector of Γ (${\gamma }_i$) defines a mixture component of the GMM, i.e., $\Gamma ={\{{\gamma }_i:{\gamma }_i=\{{\omega }_i,{\mu }_i,{\sigma }_i\}\}}_{i=1}^{{\mathrm{N}}_{\mathrm{G}}}$. ω is the weight. μ is the mean value. σ is the standard deviation. The function of each component $g\left(x\right|\mu ,\sigma )$ conforms to a Gaussian function:

$$g\left(x\right|\mu ,\sigma )=e^{-\frac{{(x-\mu )}^2}{2{\sigma }^2}}$$

(2)

The NLS method is adopted to estimate all of the parameters (ω, μ, and σ) of the mixture components of the GMM. The number of these mixture components is fixed during each step of parameter estimation. Thus, in this section, the PDF is simplified as $f\left(x\right|{\mathrm{N}}_{\mathrm{G}};\Gamma )\stackrel{{\mathrm{N}}_{\mathrm{G}}}{\to }f\left(x\right|\Gamma )$. Suppose a set of data points for ramping features (x) and their actual probabilities (p), i.e., ($x_1,p_1$), ..., ($x_j,p_j$), ..., ($x_{{\mathrm{N}}_x},p_{{\mathrm{N}}_x}$); the objective function, $S_2$, of the NLS is to minimize the sum of the squares of fitting residuals, given by:

$$S_2=\sum _{j=1}^{\mathrm{N}x}{r}_j^2=\sum _{j=1}^{\mathrm{N}x}{\left[p_j-f\left(x_j\right|\Gamma )\right]}^2$$

(3)

where $r_j$ is the residual of the j-th data point.

Taking the mean value, μ, as an example, the minimum value of $S_2$ occurs when the gradient is zero, given by:

$$\frac{\partial S_2}{\partial {\mu }_i}=2\sum _{j=1}^{{\mathrm{N}}_x}{r}_j\frac{\partial r_j}{\partial {\mu }_i}=-2\sum _{j=1}^{{\mathrm{N}}_x}{r}_j\frac{\partial f\left(x_j\right|\Gamma )}{\partial {\mu }_i}=0$$

(4)

Since the model contains $3\times {\mathrm{N}}_{\mathrm{G}}$ parameters, there are $3\times {\mathrm{N}}_{\mathrm{G}}$ gradient equations. Then, each mean value of the GMM, ${\mu }_i$, is refined iteratively by the successive approximation:

$${\mu }_i\approx {\mu }_i^{k+1}={\mu }_i^k+\Delta {\mu }_i$$

(5)

where k is an iteration number and the increment, $\Delta {\mu }_i$, is known as the step size. At each iteration, the GMM is linearized by approximating to a first-order Taylor series expansion:

$$f\left(x\right|\Gamma )\approx f\left(x\right|{\Gamma }^k)+\sum _{i=1}^{{\mathrm{N}}_{\mathrm{G}}}\frac{\partial f\left(x\right|{\Gamma }^k)}{\partial {\omega }_i}\Delta {\omega }_i+\sum _{i=1}^{{\mathrm{N}}_{\mathrm{G}}}\frac{\partial f\left(x\right|{\Gamma }^k)}{\partial {\mu }_i}\Delta {\mu }_i+\sum _{i=1}^{{\mathrm{N}}_{\mathrm{G}}}\frac{\partial f\left(x\right|{\Gamma }^k)}{\partial {\sigma }_i}\Delta {\sigma }_i$$

(6)

where the derivatives $\partial f\left(x\right|{\Gamma }^k)/\partial {\omega }_i$, $\partial f\left(x\right|{\Gamma }^k)/\partial {\mu }_i$, and $\partial f\left(x\right|{\Gamma }^k)/\partial {\sigma }_i$ consist of the Jacobian matrix, $\mathbf{J}$. Each derivative is analytically deduced by:

$$\begin{array}{c}\hfill J_{ij}^{\omega }=\frac{\partial f\left(x_j\right|{\Gamma }^k)}{\partial {\omega }_i}=e^{-\frac{{(x_j-{\mu }_i^k)}^2}{2{\left({\sigma }_i^k\right)}^2}}\end{array}$$

(7a)

$$\begin{array}{c}\hfill J_{ij}^{\mu }=\frac{\partial f\left(x_j\right|{\Gamma }^k)}{\partial {\mu }_i}=\frac{{\omega }_i^k(x_j-{\mu }_i^k)}{{\left({\sigma }_i^k\right)}^2}{e}^{-\frac{{(x_j-{\mu }_i^k)}^2}{2{\left({\sigma }_i^k\right)}^2}}\end{array}$$

(7b)

$$\begin{array}{c}\hfill J_{ij}^{\sigma }=\frac{\partial f\left(x_j\right|{\Gamma }^k)}{\partial {\sigma }_i}=\frac{{\omega }_i^k{(x_j-{\mu }_i^k)}^2}{{\left({\sigma }_i^k\right)}^3}{e}^{-\frac{{(x_j-{\mu }_i^k)}^2}{2{\left({\sigma }_i^k\right)}^2}}\end{array}$$

(7c)

$$\begin{array}{c}\hfill \mathbf{J}={\left[{\mathbf{J}}^{\omega }{\mathbf{J}}^{\mu }{\mathbf{J}}^{\sigma }\right]}^{\mathbf{T}},\forall J_{ij}^{\omega }\in {\mathbf{J}}^{\omega },\forall J_{ij}^{\mu }\in {\mathbf{J}}^{\mu },\forall J_{ij}^{\sigma }\in {\mathbf{J}}^{\sigma }\end{array}$$

(7d)

Since the iterative residuals are given by: $\Delta p_j=p_j-f\left(x_j\right|{\Gamma }^k)$, the original residuals are rearranged by:

$$\begin{array}{cc}\hfill r_j& =[p_j-f\left(x_j\right|{\Gamma }^k)]+[f\left(x_j\right|{\Gamma }^k)-f\left(x_j\right|\Gamma )]\hfill \\ & =\Delta p_j-\sum _{s=1}^{{\mathrm{N}}_{\mathrm{G}}}{J}_{sj}^{\omega }\Delta {\omega }_s-\sum _{s=1}^{{\mathrm{N}}_{\mathrm{G}}}{J}_{sj}^{\mu }\Delta {\mu }_s-\sum _{s=1}^{{\mathrm{N}}_{\mathrm{G}}}{J}_{sj}^{\sigma }\Delta {\sigma }_s\hfill \end{array}$$

(8)

Then, substituting these expressions in (7) and (8) into the gradient equations in (4), we can rearrange and get the normal equations:

$$\begin{array}{c}\hfill \sum _{j=1}^{{\mathrm{N}}_x}\sum _{s=1}^{{\mathrm{N}}_{\mathrm{G}}}{J}_{ij}^{\mu }(J_{sj}^{\omega }\Delta {\omega }_s+J_{sj}^{\mu }\Delta {\mu }_s+J_{sj}^{\sigma }\Delta {\sigma }_s)=\sum _{j=1}^{{\mathrm{N}}_x}{J}_{ij}^{\mu }\Delta p_j\end{array}$$

(9a)

$$\begin{array}{c}\hfill \sum _{j=1}^{{\mathrm{N}}_x}\sum _{s=1}^{{\mathrm{N}}_{\mathrm{G}}}{J}_{ij}^{\omega }(J_{sj}^{\omega }\Delta {\omega }_s+J_{sj}^{\mu }\Delta {\mu }_s+J_{sj}^{\sigma }\Delta {\sigma }_s)=\sum _{j=1}^{{\mathrm{N}}_x}{J}_{ij}^{\omega }\Delta p_j\end{array}$$

(9b)

$$\begin{array}{c}\hfill \sum _{j=1}^{{\mathrm{N}}_x}\sum _{s=1}^{{\mathrm{N}}_{\mathrm{G}}}{J}_{ij}^{\sigma }(J_{sj}^{\omega }\Delta {\omega }_s+J_{sj}^{\mu }\Delta {\mu }_s+J_{sj}^{\sigma }\Delta {\sigma }_s)=\sum _{j=1}^{{\mathrm{N}}_x}{J}_{ij}^{\sigma }\Delta p_j\end{array}$$

(9c)

The normal equations are written in the matrix notation as:

$${\mathbf{J}}^{\mathbf{T}}\mathbf{J}\left[\begin{array}{c}\Delta \mathit{\omega }\hfill \\ \Delta \mathit{\mu }\hfill \\ \Delta \mathit{\sigma }\hfill \end{array}\right]={\mathbf{J}}^{\mathbf{T}}\Delta \mathbf{p}$$

(10)

Since the estimated initial parameters may be far from the optimum, Equation (10) is improved by using the trust-region algorithm [25], given by:

$$[{\mathbf{J}}^{\mathbf{T}}\mathbf{J}+\lambda diag\left({\mathbf{J}}^{\mathbf{T}}\mathbf{J}\right)]\left[\begin{array}{c}\Delta \mathit{\omega }\hfill \\ \Delta \mathit{\mu }\hfill \\ \Delta \mathit{\sigma }\hfill \end{array}\right]={\mathbf{J}}^{\mathbf{T}}\Delta \mathbf{p}$$

(11)

where $diag\left({\mathbf{J}}^{\mathbf{T}}\mathbf{J}\right)$ is the diagonal matrix with the same diagonal as ${\mathbf{J}}^{\mathbf{T}}\mathbf{J}$ and λ is used to control the trust-region size.

After estimating the parameters of each mixture component of the GMM, the next step is to determine the optimal number of mixture components, ${\mathrm{N}}_{\mathrm{G},\mathrm{opt}}$, by minimizing the Euclidean distance between the actual PDF, ${\mathrm{PDF}}_{\mathrm{A}}$, and the PDF of GMM, f, i.e., $f\left(x\right|{\mathrm{N}}_{\mathrm{G}})\stackrel{{\mathrm{N}}_{\mathrm{G},\mathrm{opt}}}{\to }f\left(x\right)$. The objective function is formulated as:

$$min\sqrt{\sum _{x\in \mathcal{X}}{[f\left(x\right|i)-{\mathrm{PDF}}_{\mathrm{A}}]}^2},i=1,2,...,{\mathrm{N}}_{\mathrm{G}}$$

(12)

The cumulative distribution function (CDF), F, is another essential statistic metric to analyze wind power forecasts due to its monotonicity, which is analytically expressed as:

$$\begin{array}{cc}\hfill F\left(x\right|{\mathrm{N}}_{\mathrm{G}};\Gamma )& ={\int }_{-\infty }^x\sum _{i=1}^{{\mathrm{N}}_{\mathrm{G}}}{\omega }_i{e}^{-\frac{1}{2}{\left[\frac{t-{\mu }_i}{{\sigma }_i}\right]}^2}dt\hfill \\ & =\sum _{i=1}^{{\mathrm{N}}_{\mathrm{G}}}[\frac{\sqrt{\pi }}{2}{\omega }_i{\sigma }_i\mathrm{erf}(\frac{{\mu }_i-x}{{\sigma }_i})]+\mathrm{C}\hfill \end{array}$$

(13)

where $\mathrm{erf}$ is the Gaussian error function and defined as:

$$\mathrm{erf}\left(x\right)=\frac{2}{\sqrt{\pi }}{\int }_0^x{e}^{-t^2}dt$$

(14)

Equation (13) is an indefinite integral with a constant $\mathrm{C}$, which can be solved by (15). Since the wind power data are normalized into the range [0, 1], it can be derived that $F(x<0)=0$ and $F(x>1)=1$. Hence, we use $x=-0.1$ (i.e., $F(-0.1)=0$) to obtain the constant, $\mathrm{C}$, given by:

$$\mathrm{C}=F(-0.1)-\sum _{i=1}^{{\mathrm{N}}_{\mathrm{G}}}[\frac{\sqrt{\pi }}{2}{\omega }_i{\sigma }_i\mathrm{erf}(\frac{-0.1-{\mu }_i}{{\sigma }_i})]$$

(15)

In addition to the PDF (f) and CDF (F) of the GMM, the inverse function of the CDF ($F^{-1}$) is deduced in this paper, formulated as:

$$\widehat{x}=F^{-1}(\sum _{i=1}^{{\mathrm{N}}_{\mathrm{G}}}[\frac{\sqrt{\pi }}{2}{\omega }_i{\sigma }_i\mathrm{erf}(\frac{{\mu }_i-x}{{\sigma }_i})]+\mathrm{C})$$

(16)

However, the inverse function in (16) cannot be analytically deduced. Instead of obtaining the analytical solution, we use the Newton–Raphson method to obtain a numerical solution. The process repeats until the range is less than the stopping threshold ε, given by:

$$\left|x_k-x_{k-1}\right|<\epsilon$$

(17)

where ε is set as $1\times {10}^{-8}$.

The developed GMM distribution could be used to fit the actual histogram distribution of the FWPG of the wind-energy storage system. The detailed information about the W-ESS is discussed in the following section.

3. Economic Dispatch Model in the Wind-Thermal-Storage System

This section develops an economic dispatch model in the wind-thermal-storage system. First, a linear programming optimization model is established for the wind-storage system considering both less usage of the ESS and less curtailment of wind generation; second, the W-ESS model is integrated into a chance-constrained economic dispatch model to construct the wind-thermal-ESS system. This chance-constrained ED model considers the operation costs of conventional thermal generators, the operation costs of wind power plants, the underestimation costs of W-ESS plants and the overestimation costs of W-ESS plants. The sequential linear programming method is used to convert the chance-constrained problem to a deterministic and linear problem considering the uncertainties of the W-ESS outputs. The nonlinear spinning reserve constraints are substituted by the linear constraints using the inverse function of the GMM CDF for all W-ESS plants.

3.1. Optimization Model of the Wind-Energy-Storage System

The uncertainty and variability of wind power generation impact the penetration of wind in the power system operations, especially for the significant wind power ramps [26,27]. These extreme cases could cause a large amount of wind power curtailments. In order to take full advantage of wind power without any curtailment, PSOs have proposed to use the ESS to mitigate the high variation. However, the massive usage of the ESS could increase the operation cost of the system. Hence, the wind-energy storage system is established with the consideration of both the economic benefits of the system and less wind curtailments. The W-ESS structure with a wind farm and an integrated ESS is shown in Figure 1.

A linear programming (LP) optimization model is established to control the integration of the ESS and wind power curtailment [18]. In order to reduce the repeating usage of the ESS and the curtailment of wind generation, the objective function of the W-ESS optimization model is formulated by:

$$\underset{p_{\mathrm{S}}^t,p_{\mathrm{C}}^t}{min}{\lambda }_{\mathrm{C}}{p}_{\mathrm{C}}^t+{\lambda }_{\mathrm{S}}\left|p_{\mathrm{S}}^t\right|$$

(18)

where $p_{\mathrm{C}}^t$ represents the non-negative curtailment of wind power generation at time t. $p_{\mathrm{S}}^t$ represents the power output of the ESS at time t. A positive value of $p_{\mathrm{S}}$ means the discharging state of the ESS, and a negative value of $p_{\mathrm{S}}$ means the charging state of the ESS, respectively. ${\lambda }_{\mathrm{C}}$ and ${\lambda }_{\mathrm{S}}$ represent different penalty multipliers. Generally, the power system balancing authorities prefer to make full use of the ESS with less wind power curtailment. Thus, the penalty multipliers of ${\lambda }_{\mathrm{C}}$ and ${\lambda }_{\mathrm{S}}$ are set as 0.9 and 0.1, respectively.

As shown in Figure 1, the FWPG, $p_{\mathrm{W}}$, of the W-ESS conforms to the equality constraints at each time point t, given by:

$$p_{\mathrm{W}}^t={\tilde{p}}_{\mathrm{W}}^t+p_{\mathrm{S}}^t-p_{\mathrm{C}}^t$$

(19)

where $p_{\mathrm{W}}^t$ represents the FWPG of the W-ESS at time t. ${\tilde{p}}_{\mathrm{W}}^t$ represents the original wind power generated by the doubly-fed induction generator (DFIG) of the wind farm or plant at time t.

The energy balance for charging and discharging the storage system is given by:

$$\begin{array}{c}\hfill S^t=\left(1-d_{\mathrm{s}}\right)S^{t-1}+{\eta }_{\mathrm{c}}{L}_{\mathrm{S}}^t\end{array}$$

(20a)

$$\begin{array}{c}\hfill S^t=\left(1-d_{\mathrm{s}}\right)S^{t-1}-{\eta }_{\mathrm{d}}{G}_{\mathrm{S}}^t\end{array}$$

(20b)

where $S^t$ and $S^{t-1}$ are the energy stored in the storage system at time t and $t-1$, respectively. $d_{\mathrm{s}}$ is the self-discharge rate for the ESS. ${\eta }_{\mathrm{c}}$ and ${\eta }_{\mathrm{d}}$ represent the charging and discharging efficiencies for the ESS, respectively. $L_{\mathrm{S}}^t$ is the storage loading capacity at time t. $G_{\mathrm{S}}^t$ is the storage generating capacity at time t.

The state of charge (SOC) of the ESS at time t is constrained by the minimum storage capacity requirement, $S_{\min }$, and maximum storage capacity requirement, $S_{\max }$, given by:

$$S_{min}\le S^t\le S_{max}$$

(21)

At each time point t, the ESS power, $p_{\mathrm{S}}^t$ must not exceed the ESS rating, $P_{\mathrm{S},max}$, as follows:

$$\left|p_{\mathrm{S}}^t\right|\le P_{\mathrm{S},max}$$

(22)

Upward and downward ramping constraints for the storage are formulated by:

$$\begin{array}{c}\hfill G_{\mathrm{S}}^t-G_{\mathrm{S}}^{t-1}\le R{U}_{\mathrm{S}}\end{array}$$

(23a)

$$\begin{array}{c}\hfill G_{\mathrm{S}}^{t-1}-G_{\mathrm{S}}^t\le R{D}_{\mathrm{S}}\end{array}$$

(23b)

where $R{U}_{\mathrm{S}}$ and $R{D}_{\mathrm{S}}$ are the upward and downward ramping rate for the ESS, respectively.

The FWPG is constrained by the upper and lower ramping power limits, $P_{\mathrm{W},max}$ and $P_{\mathrm{W},min}$, given by:

$$P_{\mathrm{W},min}\le {\tilde{p}}_{\mathrm{W}}^t+p_{\mathrm{S}}^t-p_{\mathrm{C}}^t\le P_{\mathrm{W},max}$$

(24)

where the curtailment of wind power, $p_{\mathrm{C}}^t$, satisfies the constraint of the maximum original wind power generation, given by:

$$0\le p_{\mathrm{C}}^t\le {\tilde{p}}_{\mathrm{W}}^t$$

(25)

When $p_{\mathrm{C}}^t=0$, it means there is no curtailment at all. When $p_{\mathrm{C}}^t={\tilde{p}}_{\mathrm{W}}^t$, it means all of the original wind power is curtailed.

To analyze the impact of the uncertainty of wind power, the W-ESS with wind farm and energy storage is modeled in the chance-constrained ED model for the power system analysis.

3.2. Chance-Constrained Economic Dispatch Model in the Wind-Thermal-Energy-Storage System

A typical chance-constrained ED model integrated by the W-ESS considering the chance constraints is modeled as a wind-thermal-energy-storage system (WTESS) in this section, as shown in (26). In this model, only spinning reserve is considered. For the sake of simplicity, there are no regulation, replacement and non-spinning reserve schedules in this study. The piece-wise linear approximation of the cost curves of generators is utilized. The objective function consists of the operation costs of conventional thermal generators, the operation costs of wind power plants, the underestimation costs of W-ESS plants and the overestimation costs of W-ESS plants. The underestimation costs mean that the actual FWPG is not completely utilized. The overestimation costs mean that the scheduled FWPG fails to be produced by the W-ESS plants. The objective function is formulated as:

$$min\sum _{i=1}^{NI}{C}_i\left(p_{\mathrm{T},i}\right)+\sum _{w=1}^{NW}{C}_w\left(p_{\mathrm{W},w}\right)+\sum _{w=1}^{NW}{C}_{\mathrm{U},w}(p_{\mathrm{A},w}-p_{\mathrm{W},w})+\sum _{w=1}^{NW}{C}_{\mathrm{O},w}(p_{\mathrm{W},w}-p_{\mathrm{A},w})$$

(26)

where i and w are the index of thermal generators and W-ESS plants. $NI$ and $NW$ are the number of thermal generators and W-ESS plants. $p_{\mathrm{T},i}$ and $p_{\mathrm{W},w}$ are the scheduled power of thermal generator i and W-ESS plant w. $p_{\mathrm{A},w}$ is the actual generation power of W-ESS plant w. $C_i$ is the cost function of thermal generator i, formulated by (27). $C_w$ is the cost function of W-ESS plant w, formulated by:

$$C_i\left(p_{\mathrm{T},i}\right)=a_i{p}_{\mathrm{T},i}^2+b_i{p}_{\mathrm{T},i}+c_i$$

(27)

$$C_w\left(p_{\mathrm{W},w}\right)=d_w{p}_{\mathrm{W},w}$$

(28)

where $a_i$, $b_i$, and $c_i$ are the fuel cost coefficients of thermal generator i. $d_w$ is the cost coefficient of W-ESS plant w. $C_{\mathrm{U},w}$ and $C_{\mathrm{O},w}$ are the underestimation and overestimation cost functions of W-ESS plant w, respectively. In order to consider the uncertainties of actual FWPG w, the probabilistic underestimation and overestimation costs, $C_{\mathrm{U},w}$ and $C_{\mathrm{O},w}$, can be expressed as the forms of expectation values developed by [28,29] and shown as below:

$$C_{\mathrm{U},w}(p_{\mathrm{A},w}-p_{\mathrm{W},w})=k_{\mathrm{U}}(p_{\mathrm{A},w}-p_{\mathrm{W},w})=k_{\mathrm{U}}{\int }_{p_{\mathrm{W},w}}^{p_{max,w}}(t-p_{\mathrm{W},w})g_w\left(t\right)dt$$

(29)

$$C_{\mathrm{O},w}(p_{\mathrm{W},w}-p_{\mathrm{A},w})=k_{\mathrm{O}}(p_{\mathrm{W},w}-p_{\mathrm{A},w})=k_{\mathrm{O}}{\int }_0^{p_{\mathrm{W},w}}(p_{\mathrm{W},w}-t)g_w\left(t\right)dt$$

(30)

where $k_{\mathrm{U}}$ and $k_{\mathrm{O}}$ are the underestimation and overestimation cost coefficients. $g_w(·)$ is the PDF of the FWPG of W-ESS plant w. In this paper, we use the GMM in Section 2 to model the PDF of the FWPG. By using the aforementioned formulations, the underestimation cost is equivalent to the expected cost of FWPG surplus, and the overestimation cost is equivalent to the expected cost of FWPG deficit. Wind plant operators would be penalized according to the electricity market rules when the FWPG deviations are not covered.

The aforementioned objective function is constrained by the following constraints. The power balance constraint is given by:

$$\sum _{i=1}^{NI}{p}_{\mathrm{T},i}+\sum _{w=1}^{NW}{p}_{\mathrm{W},w}=L$$

(31)

where L is the system demand as constant.

The power generation limit constraints of thermal generators and W-ESS are given by:

$$p_{min,i}\le p_{\mathrm{T},i}\le p_{max,i}$$

(32)

$$0\le p_{\mathrm{W},w}\le p_{max,w}$$

(33)

where $p_{min,i}$ and $p_{max,i}$ are the minimum and maximum output of thermal generator i. $p_{max,w}$ is the maximum output of W-ESS plant w.

The spinning reserve constraints are given by:

$$0\le r_{\mathrm{u},i}\le min\{p_{max,i}-p_{\mathrm{T},i},r_{\mathrm{u},i}^{max}\}$$

(34)

$$0\le r_{\mathrm{d},i}\le min\{p_{\mathrm{T},i}-p_{min,i},r_{\mathrm{d},i}^{max}\}$$

(35)

$$\sum _{i=1}^{NI}{r}_{\mathrm{u},i}\ge \sum _{w=1}^{NW}(p_{\mathrm{W},w}-p_{\mathrm{A},w})$$

(36)

$$\sum _{i=1}^{NI}{r}_{\mathrm{d},i}\ge \sum _{w=1}^{NW}(p_{\mathrm{A},w}-p_{\mathrm{W},w})$$

(37)

where $r_{\mathrm{u},i}$ and $r_{\mathrm{d},i}$ are the scheduled amount of up- and down-spinning reserve services provided by thermal generator i. $r_{\mathrm{u},i}^{max}$ and $r_{\mathrm{d},i}^{max}$ are the maximum amount of up- and down-spinning reserve services provided by thermal generator i within a specific time interval. Considering the uncertainties of actual FWPG, $p_{\mathrm{A},w}$, the probabilistic constraints in (36) and (37) are replaced by the chance-constrained formulations [30,31], shown as below:

$$\mathrm{Pr}\left\{\sum _{i=1}^{NI}{r}_{\mathrm{u},i}\ge \sum _{w=1}^{NW}(p_{\mathrm{W},w}-p_{\mathrm{A},w})\right\}\ge C{L}_{\mathrm{u}}$$

(38)

$$\mathrm{Pr}\left\{\sum _{i=1}^{NI}{r}_{\mathrm{d},i}\ge \sum _{w=1}^{NW}(p_{\mathrm{A},w}-p_{\mathrm{W},w})\right\}\ge C{L}_{\mathrm{d}}$$

(39)

where $C{L}_{\mathrm{u}}$ and $C{L}_{\mathrm{d}}$ are the confidence levels to prepare sufficient up- and down-spinning reserve services for the system operation. $\mathrm{Pr}\left\{·\right\}$ represents the probability of a specific chance-constraint. The probabilities should be greater than the confidence levels of $C{L}_{\mathrm{u}}$ and $C{L}_{\mathrm{d}}$, respectively.

3.3. Approach for Solving the Chance-Constrained Economic Dispatch Problem

The chance-constrained ED model is established in Section 3.2. The conventional approaches have difficulty solving this stochastic nonlinear ED model. Hence, we use the sequential linear programming method to convert the chance-constrained problem to a deterministic and linear problem considering the FWPG uncertainties. The detailed solution procedure of the sequential linear programming method can be seen in [32]. Assuming an operation point in the iteration k of the solution procedure, $(p_{\mathrm{W},1}^k,p_{\mathrm{W},2}^k,...,p_{\mathrm{W},NW}^k,p_{\mathrm{T},1}^k,p_{\mathrm{T},2}^k,...,p_{\mathrm{T},NI}^k)$, the underestimation and overestimation costs of the objective function in (26) can be analytically linearized by using the Taylor series expansion, given by:

$$\frac{\partial C_{\mathrm{U},w}}{\partial p_{\mathrm{W},w}}=k_{\mathrm{U}}[G_w\left(p_{\mathrm{W},w}\right)-G_w\left(p_{\mathrm{A},w}\right)]$$

(40)

$$\frac{\partial C_{\mathrm{O},w}}{\partial p_{\mathrm{W},w}}=k_{\mathrm{O}}[G_w\left(p_{\mathrm{W},w}\right)-G_w\left(0\right)]$$

(41)

where $G_w(·)$ represents the CDF of GMM, which is discussed in Section 2 in detail.

It is difficult to directly use the nonlinear constraints in (38) and (39) to solve the chance-constrained ED problem. The linear constraints are adopted to substitute the nonlinear ones, which are shown as below:

$$\sum _{w=1}^{NW}{p}_{W,w}-\sum _{i=1}^{NI}{r}_{\mathrm{u},i}\le G_{\Sigma }^{-1}(1-C{L}_{\mathrm{u}})$$

(42)

$$\sum _{w=1}^{NW}{p}_{W,w}+\sum _{i=1}^{NI}{r}_{\mathrm{d},i}\ge G_{\Sigma }^{-1}\left(C{L}_{\mathrm{d}}\right)$$

(43)

where $G_{\Sigma }^{-1}(·)$ represents the inverse function of the GMM CDF of the FWPG for all W-ESS plants, which is also discussed in Section 2 in detail.

Based on the aforementioned formulations, the nonlinear chance-constrained ED problem can be converted to a deterministic linear problem by using the sequential linear programming, given by:

$$min\Lambda \Delta L+\sum _{i=1}^{NI}{\lambda }_{\mathrm{T},i}^k(p_{\mathrm{T},i}-p_{\mathrm{T},i}^k)+\sum _{w=1}^{NW}{\lambda }_{\mathrm{W},w}^k(p_{\mathrm{W},w}-p_{\mathrm{W},w}^k)$$

(44)

where $\Delta L$ represents the amount of shedding load and is used to ensure the existence of a resolution in the iteration process of sequential linear programming. Λ represents the penalty factor defined by the user and satisfies $\Lambda \gg max\{{\lambda }_{\mathrm{T},i}^k,{\lambda }_{\mathrm{W},w}^k\}$. After removing the constant part in (44), this objective function can be further simplified as:

$$min\lambda \Delta L+\sum _{i=1}^{NI}{\lambda }_i^k{p}_{\mathrm{T},i}+\sum _{w=1}^{NW}{\lambda }_{\mathrm{W},w}^k{p}_{\mathrm{W},w}$$

(45)

Except constraints in (32)–(35) and (42)–(43), this linear problem is additionally subject to:

$$\left\{\begin{array}{c}{\displaystyle \sum _{i=1}^{NI}}{p}_{\mathrm{T},i}+{\displaystyle \sum _{w=1}^{NW}}{p}_{\mathrm{W},w}=L-\Delta L\hfill \\ 0\le \Delta L\le L\hfill \\ \left|p_{\mathrm{T},i}-p_{\mathrm{T},i}^k\right|\le \Delta p_{\mathrm{T},i}^k\hfill \\ \left|p_{\mathrm{W},w}-p_{\mathrm{W},w}^k\right|\le \Delta p_{\mathrm{W},w}^k\hfill \end{array}\right.$$

(46)

where $\Delta p_{\mathrm{T},i}^k$ is the step size of $p_{\mathrm{T},i}^k$ in iteration k. $\Delta p_{\mathrm{W},w}^k$ is the step size of $p_{\mathrm{W},w}^k$ in iteration k. Both step sizes are defined by the user. The constant coefficient ${\lambda }_{\mathrm{T},i}^k$ can be deduced by the cost function of thermal generators in (27), given by:

$${\lambda }_{\mathrm{T},i}^k={\left.\frac{\partial C_{\mathrm{T},i}}{\partial p_{\mathrm{T},i}}\right|}_{p_{\mathrm{T},i}=p_{\mathrm{T},i}^k}=2a_i{p}_{\mathrm{T},i}^k+b_i$$

(47)

The constant coefficient ${\lambda }_{\mathrm{W},w}^k$ can be deduced by the cost function of W-ESS plants in (28), the underestimation cost in (40) and the overestimation cost in (41), given by:

$$\begin{array}{cc}\hfill {\lambda }_{\mathrm{W},w}^k& ={\left.\left(\frac{\partial C_{\mathrm{W},w}}{\partial p_{\mathrm{W},w}}+\frac{\partial C_{\mathrm{U},w}}{\partial p_{\mathrm{W},w}}+\frac{\partial C_{\mathrm{O},w}}{\partial p_{\mathrm{W},w}}\right)\right|}_{p_{\mathrm{W},w}=p_{\mathrm{W},w}^k}\hfill \\ & =d_w-k_{\mathrm{U}}{G}_w\left(p_{\mathrm{A},w}\right)-k_{\mathrm{O}}{G}_w\left(0\right)+(k_{\mathrm{U}}+k_{\mathrm{O}})G_w\left(p_{\mathrm{W},w}^k\right)\hfill \end{array}$$

(48)

Based on the aforementioned formulations, the sequential linear programming can be used to solve the linear ED problem in (44) and obtain the resolution of $(p_{\mathrm{W},1}^{k+1},p_{\mathrm{W},2}^{k+1},...,p_{\mathrm{W},NW}^{k+1},p_{\mathrm{T},1}^{k+1},p_{\mathrm{T},2}^{k+1},...,p_{\mathrm{T},NI}^{k+1})$ in the iteration $(k+1)$, which is near the new operation point. This calculation process repeats until the range is less than the convergence criterion parameter, ϵ, given by:

$$\left|(p_{\mathrm{W},1}^{k+1},p_{\mathrm{W},2}^{k+1},...,p_{\mathrm{W},NW}^{k+1},p_{\mathrm{T},1}^{k+1},p_{\mathrm{T},2}^{k+1},...,p_{\mathrm{T},NI}^{k+1})-(p_{\mathrm{W},1}^k,p_{\mathrm{W},2}^k,...,p_{\mathrm{W},NW}^k,p_{\mathrm{T},1}^k,p_{\mathrm{T},2}^k,...,p_{\mathrm{T},NI}^k)\right|\le ϵ$$

(49)

The flowchart of the proposed approach of the WTESS model is shown in Figure 2.

4. Case Studies and Results

The developed WTESS model is used on the IEEE 30-bus, system to verify its effectiveness. In this system, there are 6 thermal generators, 1 W-ESS plant and 41 transmission lines. Detailed information about this system can be seen in [33]. Line flow capacities of the IEEE 30-bus model test system are listed in

Table 1. Line flow capacities of the IEEE 30-bus model test system.

Line No.	Cap. (MW)	Line No.	Cap. (MW)	Line No.	Cap. (MW)	Line No.	Cap. (MW)	Line No.	Cap. (MW)	Line No.	Cap. (MW)
1	20.0	8	15.0	15	70.0	22	50.0	29	20.0	36	50.0
2	20.0	9	15.0	16	50.0	23	50.0	30	20.0	37	50.0
3	20.0	10	55.0	17	70.0	24	50.0	31	20.0	38	50.0
4	20.0	11	50.0	18	50.0	25	20.0	32	50.5	39	50.0
5	20.0	12	50.0	19	50.0	26	20.0	33	50.0	40	50.0
6	20.0	13	30.0	20	50.0	27	20.0	34	50.0	41	50.0
7	10.0	14	50.0	21	50.0	28	20.0	35	15.0

, where the maximum line flow capacities are given by 110% of the corresponding standard value. Parameters of the conventional thermal generators are listed in

Table 2. Parameters of the conventional thermal generators.

Gen. No.	Fuel Cost Coefficients ($/h)			Output Limits (p.u.)		Reserve Limits (p.u.)
	${\mathit{a}}_{\mathit{i}}$	${\mathit{b}}_{\mathit{i}}$	${\mathit{c}}_{\mathit{i}}$	${\mathit{p}}_{\mathbf{min},\mathit{i}}$	${\mathit{p}}_{\mathbf{max},\mathit{i}}$	${\mathit{r}}_{\mathit{u},max,\mathit{i}}$	${\mathit{r}}_{\mathit{d},max,\mathit{i}}$
1	100	200	10	0.4	1	0.2	0.2
2	120	150	10	0.4	1	0.2	0.2
3	40	180	20	0.4	1	0.2	0.2
4	60	100	10	0.4	1	0.2	0.2
5	40	180	20	0.4	1	0.2	0.2
6	40	180	20	0.4	1	0.2	0.2

. The wind power data used for statistical analysis are collected from four wind plants in China. The simulation data represent wind power generation and forecasts spanning from 1 January 2010–31 December 2012 with a 1-min time resolution. Four regions are selected for the case study, denoted as WF-I, WF-II, WF-III and WF-IV, respectively. The rated capacity with approximately 150 MW is chosen. The total number of samples is 6,312,960, which is sufficiently large for statistical analysis of probabilistic wind power generation. The total load demand of the IEEE 30-bus system to be satisfied is approximately 283.4 MW. The solver used in the optimization model of the WTESS is the IBM ILOG CPLEX Optimizer 12.6. The computation is run in MATLAB R2015a on a Core-i5 2.60-GHz laptop.

4.1. Statistical Comparisons of the PDF and CDF of the FWPG

In order to compare the performance of different distributions with the GMM distribution developed in this paper, three distributions are used for comparison, including: normal, logistic and versatile distributions. Normal and logistic distributions are widely used for statistical analysis. The versatile distribution is relatively new for statistical analysis and has shown good performance in fitting the distribution of wind power [29]. In the following case studies, we use the aforementioned three distributions to verify the effectiveness of the GMM distribution. The FWPG data of the WTESS are used for simulation and comparison.

Based on the wind power data at 1-min time resolution, we calculate the 5-min time resolution data in an average manner to compare the probability distribution of four wind farms. The effectiveness of the PDF and the CDF of the GMM distribution is verified based on the 5-min time resolution data. As shown in Figure 3, the normal and logistic distributions fail to track the actual probability values, especially for the peak values. Though the versatile distribution could fit the actual probability values accurately, it performs worse than the GMM distribution with lower fitting accuracy. To quantitatively compare the fitting performance, three metrics are adopted to assess the distribution accuracy, including: the mean absolute error (MAE), chi-square goodness-of-fit (GOF) and the root mean square error (RMSE). All three metrics are calculated between the actual PDF and the PDFs of fitting distributions. The mathematical expressions of the metrics are shown in (50)–(52). Numerical results are listed in

Table 3. Metrics of four distributions estimated for the PDF.

Dist.	WF-I			WF-II			WF-III			WF-IV
	MAE	GOF	RMSE	MAE	GOF	RMSE	MAE	GOF	RMSE	MAE	GOF	RMSE
Norm.	0.11	8.11	0.15	0.17	6.44	0.18	0.16	7.48	0.21	0.17	5.32	0.19
Log.	0.09	6.52	0.12	0.08	3.18	0.14	0.08	4.32	0.18	0.12	4.18	0.15
Vers.	0.04	0.44	0.06	0.05	0.84	0.07	0.05	1.01	0.11	0.08	1.08	0.08
GMM	0.02	0.11	0.03	0.02	0.23	0.03	0.01	0.19	0.05	0.01	0.23	0.04

to show the metrics of four distributions estimated for the PDF.

$$MAE=\frac{1}{N}\sum _{n=1}^N\left|F_{act}\left(x_n\right)-F_{sim}\left(x_n\right)\right|$$

(50)

$$GOF=\sum _{n=1}^N{\left[F_{act}\left(x_n\right)-F_{sim}\left(x_n\right)\right]}^2/F_{sim}\left(x_n\right)$$

(51)

$$RMSE=\sqrt{\frac{1}{N}\sum _{n=1}^N{\left[F_{act}\left(x_n\right)-F_{sim}\left(x_n\right)\right]}^2}$$

(52)

As shown in Table 3, the metrics of the GMM distribution are the smallest, compared with those of the normal, logistic and versatile distributions. Taking the MAE metrics as an example, the normal distribution presents the metrics values with the range 0.11∼0.17; the logistic distribution presents the metrics values with the range 0.08∼0.12; the versatile distribution presents the metrics values with the range 0.04∼0.08; and the GMM distribution presents the metrics values with the smallest range 0.01∼0.03.

As shown in Figure 4, the normal, logistic and versatile distributions fail to track the actual CDF values. Similarly to the comparisons of the PDF, it is also shown that the GMM follows the actual CDF curve significantly better than other distributions (normal, logistic and versatile). Numerical results are listed in

Table 4. Metrics of four distributions estimated for the CDF.

Dist.	WF-I			WF-II			WF-III			WF-IV
	MAE	GOF	RMSE	MAE	GOF	RMSE	MAE	GOF	RMSE	MAE	GOF	RMSE
Norm.	0.15	5.47	0.21	0.21	8.77	0.25	0.21	5.33	0.17	0.22	7.99	0.16
Log.	0.11	4.35	0.17	0.14	5.69	0.19	0.17	3.88	0.12	0.15	5.84	0.12
Vers.	0.07	1.38	0.09	0.09	1.39	0.12	0.15	1.11	0.09	0.11	3.25	0.05
GMM	0.04	0.49	0.05	0.04	0.41	0.07	0.09	0.56	0.03	0.07	0.99	0.02

to show the metrics of four distributions estimated for the CDF. The metrics of the GMM distribution are also the smallest, compared with those of the normal, logistic and versatile distributions. Taking the GOF metrics as an example, the normal distribution presents the metrics values with the range 5.32∼8.11; the logistic distribution presents the metrics values with the range 3.18∼6.52; the versatile distribution presents the metrics values with the range 0.44∼1.08; and the GMM distribution presents the metrics values with the smallest range 0.11∼0.23. Based on the results of Table 3 and Figure 4, the GMM distribution could fit the actual distribution most accurately for all four wind farms.

4.2. Results of the Chance-Constrained ED Model

Each FWPG case is performed on the IEEE 30-bus system by using the sequential linear programming. For the simplicity of comparison, the actual FWPG data of the W-ESS, normal, logistic, versatile and GMM distributions are studied for comparison based on the numerical computation. As shown in Figure 5, the solutions of the chance-constrained ED model with scheduled FWPG are compared. It is observed that the solutions obtained by using the GMM distribution are approximately the same as those obtained from the distribution of the actual FWPG data. However, for the normal, logistic and versatile distributions, the accuracy of the chance-constrained ED model of the WTESS system is sacrificed due to the worse fitting performance of each distribution model. However, the solutions obtained by using the GMM distribution have the smallest fitting errors.

Figure 6 compares the operation costs of four cases (i.e., actual, normal, versatile and GMM) with six wind power penetrations, which consist of 8.66%, 12.99%, 17.32%, 21.65%, 25.98% and 30.31%, respectively. The wind penetration level is calculated by the percentage of the total wind power accounting for the total load of the system. For all wind power penetrations in Figure 6, the case using the actual FWPG data obtains the least operation cost; the case using the normal distribution obtains the most operation cost; and the case using the GMM distribution obtains less operation costs, comparing to the normal and versatile distributions. This is because the GMM distribution could fit the actual FWPG data of the W-ESS most accurately. Comparing to the cost using the normal distribution, the saving cost using the GMM distribution is the most, $8,336, with the saving percentage 25.05% when the wind power penetration level is 25.98%. Comparing to the cost using the versatile distribution, the saving cost using the GMM distribution is the most, $4,449, with the saving percentage 14.51% when the wind power penetration level is 21.65%.

In addition, the operation costs decrease with the increasing penetration level of wind power. When the wind power penetration levels increase from 8.66%–30.31%, the average operation costs decrease approximately 36.29%, 31.21%, 34.86% and 38.53%, respectively. This finding is meaningful for the integration of renewable resources, such as wind power.

4.3. Comparisons with and without ESS

Figure 7 shows the curves of the original wind power generation (OWPG), ${\tilde{p}}_{\mathrm{W}}$, generated by the DFIG of the wind farm or plant, and the FWPG, $p_{\mathrm{W}}$, of the W-ESS in (19), respectively. The simulation time is spanning from 0:00–22:00. It is observed that the FWPG is significantly smoother than the OWPG for the whole simulation time window. Figure 8 shows the power output of the ESS model $p_{\mathrm{S}}$ with simultaneous variations of charging and discharging power. The positive value above the zero line (in red) means the discharging state of the ESS. The negative value below the zero line (in red) means the charging state of the ESS.

To compare the impact of the ESS on the chance-constrained ED of the WTESS model based on the GMM distribution, the wind power data of four wind farms are input into the WTESS model with and without considering the ESS. Figure 9 shows the operation costs comparisons with and without (w/o) the ESS of four wind farms considering different wind power penetration levels. Six wind power penetrations are utilized for statistical analysis, which consist of 8.66%, 12.99%, 17.32%, 21.65%, 25.98% and 30.31%, respectively. As can be seen in Figure 9, all of the operation costs decrease with the wind power penetration level. The usage of the ESS could help effectively decrease the operation costs to different degrees.

Numerical comparison results are illustrated in

Table 5. Saving costs and percentages using the ESS in the WTESS with six wind power penetration levels on four W-ESS plants.

W-ESS Plant	Saving Metrics	Wind Power Penetration Levels
		8.66%	12.99%	17.32%	21.65%	25.98%	30.31%
WF-I	Cost ($)	20,533	21,609	15,170	17,423	15,708	16,072
	Percentage (%)	35.55	38.29	31.21	36.43	34.61	36.58
WF-II	Cost ($)	26,854	20,166	21,533	22,075	16,733	16,817
	Percentage (%)	38.75	32.32	35.45	37.32	32.76	35.64
WF-III	Cost ($)	25,621	26,009	15,287	15,216	17,881	17,976
	Percentage (%)	38.86	39.61	28.03	29.10	34.46	36.11
WF-IV	Cost ($)	21,216	17,315	18,776	23,151	19,752	18,963
	Percentage (%)	33.19	29.97	33.68	41.71	38.71	38.29

to show the saving costs and percentages using the ESS in the WTESS with six wind power penetration levels on four W-ESS plants. In W-ESS plant WF-I, the most saving cost is $21,609 with the saving percentage 38.29% when the penetration level is 12.99%. In W-ESS plant WF-II, the most saving cost is $26,854 with the saving percentage 38.75% when the penetration level is 8.66%. In W-ESS plant WF-III, the most saving cost is $26,009 with the saving percentage 39.61% when the penetration level is 12.99%. In W-ESS plant WF-IV, the most saving cost is $23,151 with the saving percentage 41.71% when the penetration level is 21.65%. Figure 9 and Table 5 verify the effectiveness of the usage of the ESS in the WTESS model.

Table 6. Comparison of the consuming time with four wind farms.

Wind Farms	Consuming Time (s)
WF-I	635
WF-II	689
WF-III	701
WF-IV	653

compares the consuming time of the proposed approach with four wind farms. It is seen that the average consuming time using the proposed approach is approximately 670 s. The simulation on the W-ESS plant WF-III takes the maximum consuming time, approximately 701 s. The simulation on the W-ESS plant WF-I takes the minimum consuming time, approximately 635 s. It is verified that the proposed approach is convenient to be used in the IEEE 30-bus system with relatively less consuming time.

5. Conclusions

The paper developed a chance-constrained economic dispatch (ED) model for the wind-thermal-energy storage system (WTESS). An optimization model with the wind power and the energy storage system (W-ESS) was first established with the consideration of both economic benefits of the system and less wind power curtailments. The original wind power generation (OWPG) was processed by the ESS to obtain the final wind power generation (FWPG). A Gaussian mixture model (GMM) distribution was adopted to characterize the probabilistic distribution function (PDF) with an analytical expression. The cumulative distribution function (CDF) was deduced by using the PDF of the GMM distribution. Then, a chance-constrained ED model integrated by the W-ESS was developed by considering both the overestimation costs and the underestimation costs of the system and solved by the sequential linear programming method. The actual wind power data in four W-ESS plants were performed on the developed chance-constrained ED model to verify the effectiveness. Numerical simulations results on the IEEE 30-bus system showed that:

(i)

The GMM distribution could fit the actual probabilistic and cumulative distributions of the FWPG of the W-ESS most accurately for all four W-ESS plants, comparing to other distributions.

(ii)

The operation costs decreased with the increasing penetration level of wind power. When the wind power penetration levels increased from 8.66%–30.31%, the average operation costs decreased approximately 31.21%∼38.53%.

(iii)

The case using the GMM distribution obtained the least operation costs, comparing to the normal and versatile distributions. The saving cost using the GMM distribution was the most, $4449, with the saving percentage 14.51% when the wind power penetration level was 21.65%.

(iv)

The usage of the ESS could help effectively decrease the operation costs to different degrees. The effectiveness of the usage of the ESS in the WTESS model was verified.