IEEE/CAA Journal of Automatica Sinica  2017, Vol. 4 Issue(2): 186-194 PDF
A Chance Constrained Optimal Reserve Scheduling Approach for Economic Dispatch Considering Wind Penetration
Yufei Tang1, Chao Luo3, Jun Yang3, Haibo He2
1. Department of Computer & Electrical Engineering and Computer Science, and the Institute for Sensing and Embedded Network Systems Engineering, Florida Atlantic University, Boca Raton, FL 33431, USA;
2. Department of Electrical, Computer and Biomedical Engineering, University of Rhode Island, Kingston, RI 02881, USA;
3. School of Electrical Engineering, Wuhan University, Wuhan 430072, China
Abstract: The volatile wind power generation brings a full spectrum of problems to power system operation and management, ranging from transient system frequency fluctuation to steady state supply and demand balancing issue. In this paper, a novel wind integrated power system day-ahead economic dispatch model, with the consideration of generation and reserve cost is modelled and investigated. The proposed problem is first formulated as a chance constrained stochastic nonlinear programming (CCSNLP), and then transformed into a deterministic nonlinear programming (NLP). To tackle this NLP problem, a three-stage framework consists of particle swarm optimization (PSO), sequential quadratic programming (SQP) and Monte Carlo simulation (MCS) is proposed. The PSO is employed to heuristically search the line power flow limits, which are used by the SQP as constraints to solve the NLP problem. Then the solution from SQP is verified on benchmark system by using MCS. Finally, the verified results are feedback to the PSO as fitness value to update the particles. Simulation study on IEEE 30-bus system with wind power penetration is carried out, and the results demonstrate that the proposed dispatch model could be effectively solved by the proposed three-stage approach.
Key words: Chance constrained     day-ahead economic dispatch     optimal reserve scheduling     particle swarm optimization (PSO)     wind power penetration
Ⅰ. INTRODUCTION

SUSTAINABLE energy generation from renewable resources, such as wind power, are developing rapidly around the world in the past decades. Because of their highly randomness and volatility characteristics, large-scale integration of these resources brings great challenges to the system secure and economic operation, as well as ancillary services and reserves scheduling [1]. Moreover, the relative poor prediction accuracy of wind generation causes a series of stability and economic issues, ranging from short-term transient stability (e.g., frequency fluctuation [2]) to long-term generation-demand balancing problems (e.g., economic dispatch [3]). In this paper, we focus on the day-ahead economic dispatch for wind generation integrated power system considering optimal reserve scheduling, and the design of effective solution algorithm.

In the power and energy community, the research for wind generation integrated power systems could be categorized into two folds. On one hand, the wind power holds the characteristic of high temporal variations, and its large-scale integration will bring impact on the system transient stability and control [4]. Examples include impact of double-fed induction generator (DFIG) based wind turbine on power systems small signal stability [5], computational intelligence (CI) based wind farm active power and reactive power damping control [6]-[8]. On the other hand, the wind power acts as an external stochastic source in the power system economic operation, and should be carefully and dedicated addressed in the planning stage. To more strategically accommodate large-scale intergraded wind power generation in economic scheduling, significant work has been finished in the market clearing modelling [9], [10] as well as stochastic optimization techniques [11]-[15]. The related investigations include stochastic power system operation by using chance constrained day-ahead scheduling [16], [17], and risk based unit commitment (UC) for day-ahead market clearing with wind power uncertainty [18], [19].

Moreover, with the decreased wind generation cost, the operators are required to fully accommodate the power from the wind without curtailment [20]. Therefore, the system operation and management will have higher requirement for the ancillary services and system reserves [21]. Along this direction, extensive studies have been carried out in the power and energy society (PES), such as procurement for load-following reserves considering flexible demand and high wind penetration [22], security-constrained scheduling based hourly reserve allocation versus demand response [23], and game theory based multi-area spinning reserve trading with wind power uncertainty [24]. Among all the proposed methods, the scenario-based stochastic programming has attracted the attention from the researchers in recent years [25]. However, we should notice that this method relies on the past experience, and its subjective and heuristic nature leaves many academics uncomfortable [26], [27].

Based on aforementioned discussion, in this paper a novel day-ahead economic dispatch model for wind power integrated power system with optimal reserves scheduling is proposed and investigated. Considering the introduced stochastic factor from the wind, chance constrained model is proposed to address this issue, and the problem is formulated as a chance constrained stochastic nonlinear programming (CCSNLP). Based on the quantile and relaxation concepts, the CCSNLP in transformed into a deterministic nonlinear programming (NLP). Then a three-stage solution framework based on particle swarm optimization (PSO), sequential quadratic programming (SQP) and Monte Carlo simulation (MCS) is proposed to solve this model effectively. The proposed model and algorithm is tested on IEEE 30-bus system with wind power penetration, and the the results demonstrate the correctness of this dispatch model as well as the effectiveness of the three-stage solution framework.

The rest of this paper is organized as follows. Section Ⅱ formulates the problem as chance constrained day-ahead economic dispatch considering reserve scheduling. Section Ⅲ presents the methodology of the proposed three-stage algorithm, including the PSO stage, the SQP solution and the verification stages. Section Ⅳ presents the simulation results of case study on IEEE 30-bus system with wind power generation. Finally, conclusions are drawn in Section Ⅴ.

Ⅱ. PROBLEM FORMULATION A. Generator Output Analysis

The power grid is a real-time system requiring the plants produce the right amount of electricity at the right time to consistently and reliably meet the load demand, such that the system frequency is maintained at the specified value. To fulfill this task, a certain amount of active power called control reserve is stored in the system [28], and the related control schemes could be categorized as primary, secondary and tertiary control.

The automatic generation control (AGC) in the power system consists of the primary and the secondary control, and the related reserve is called AGC reserve or frequency-response reserve. The primary control is provided by participated spinning generators, which response to disturbances-caused frequency deviations from the nominal value according to their speed droop characteristics. Meanwhile, the objective of the secondary control is to help the primary control to clear the frequency error, and bring the frequency back to its nominal value as soon as possible. Moreover, in the inter-connected power system, the secondary control is also responsible for maintaining the power on the tie-line to the pre-defined values. The third frequency control scheme is called tertiary control, and is usually activated manually such that the used primary and secondary control reserves are released after a large disturbance (e.g., $N-$ 1 contingency). The reserve related to this control is called contingency reserve or replacement reserve. This control scheme occurs 10-15 min after a serious system contingency and its goal is to set-up new post-contingency operating points. In tertiary control, non-spinning generation reserves can also be used [29], such as energy storage device and small hydropower plant, since it allows for a longer time lag before their deployment.

In this paper, we are focusing on optimal reserve scheduling and we only consider the reserves available in the generators. Fig. 1 shows the AGC-participated generator output power and the reserves in the normal and $N-$ 1 conditions [30]. The following two constraints should be satisfied in the normal condition:

 Download: larger image Fig. 1 Generator output with reserve scheduling at time period t.
 $$$\begin{cases} P_{Gi, t} + U_{Gi, t} \le \min ( P_{Gi, \max }, \, A_{Gi, \max } ) \\ P_{Gi, t} - D_{Gi, t} \ge \max ( P_{Gi, \min }, \, A_{Gi, \min } ) \end{cases}$$$ (1)

where the meaning of the symbols are shown in the figure. For the $N-$ 1 condition (i.e., one generator is disconnected from the main grid), the following constraint should still be satisfied for system operation:

 $$$P_{Gi, t} + U_{Gi, t} + R_{Gi, t} \le \min (P_{Gi, \max }, \, A_{Gi, \max }).$$$ (2)

We should also notice that, the deployment of system reserves subject to the generator ramping up and ramping down capability, and could be modelled as:

 $$$\begin{cases} U_{Gi, t} \le T_{1i} r_{ui} \\ D_{Gi, t} \le T_{1i} r_{di} \\ R_{Gi, t} \le T_{2i} r_{di} \end{cases}$$$ (3)

where $r_{ui}$ and $r_{di}$ are the ramping up and ramping down rate for the $i$th generator, respectively; and $T_{1i}$ and $T_{2i}$ are the ramping up and ramping down time for the $i$th generator. In this paper, we are using $T_{1}=$ 5 min and $T_{2}=$ 15 min for all the generators.

B. AGC Based Reserves Representation

Before we formulate the optimal AGC based reserve scheduling, we will first introduce the stochastic wind generation model and the load demand model. We assume the wind speed is composed as:

 $$$v_t = v_{f, t} + e_{W, t}$$$ (4)

where $v_{f, t}$ is the forecasted wind speed at time step $t$, and $e_{W, t}$ is the wind speed forecast error, which is represented as a Gaussian distribution of $\mathcal{N}(0, \sigma_{W, t})$. Based on this wind speed model, the wind farm output power is calculated as:

 $$$P_{W, t} = \begin{cases} {P_{Wr}, }&{v_r < v_t \le v_{\rm out} } \\ {\displaystyle\frac{{v_t^3 - v_{in}^3 }}{{v_r^3 - v_{in}^3 }} \cdot P_{Wr}, }&{v_{\rm in} < v_t \le v_r } \\ {0, }&\mbox{otherwise} \end{cases}$$$ (5)

where $v_{\rm in}$, $v_{r}$ and $v_{\rm out}$ are the cut-in, rated, and cut-out wind speeds, respectively, and $P_{Wr}$ is the rated power of the wind farm.

Similar as the stochastic wind power model, we formulated the total load demand as:

 $$$P_{L, t} = P_{Lf, t} + e_{L, t}$$$ (6)

where $P_{Lf, t}$ is the forecasted load demand at time step $t$, and $e_{L, t}$ is the load forecast error, which is also represented as a Gaussian distribution of $\mathcal{N}(0, \sigma_{L, t})$. For each bus in the system, we assume the total load is distributed to each bus according to the following:

 $$$P_{Ls, t} = d_s \cdot P_{L, t}, \quad 1 \le s \le N_s$$$ (7)

where $d_{s}$ is the load coefficient of bus number $s$, and $N_{s}$ is the total bus number.

The load forecasting and wind power forecasting error will be the sources of the system frequency instability, and should be compensated by the AGC system, which is shown in Fig. 2. As we have mentioned above, there are multiple generators in the system participating in the AGC task. The power balancing task are sharing among all the generators according to the following equation:

 Download: larger image Fig. 2 Overall architecture of the AGC based system power balancing.
 $$$P_{Gi, t}' = P_{Gi, t} + k_{i, t} \left( {P_{L, t} - P_{W, t} - \sum\limits_{j = 1}^{N_G } {P_{Gj, t} } } \right)$$$ (8)

where $P_{Gi, t}'$ is the $i$th generator actual output with AGC adjustment in time period $t$, $N_{G}$ is the total number of generators in the system, and $k_{i, t}$ is the distribution coefficient. Usually, the distribution vector $\pmb{k}_{t}$ satisfies the following condition:

 $$$k_{1, t} + k_{2, t} + \cdots + k_{N_{\rm G}, t} = 1, \quad 0 \le k_{i, t} \le 1.$$$ (9)

When the system is in $N-$ 1 condition, the contingency reserves in the rest of the generators will be activated by the tertiary control, and the actual output of each generator will be updated as:

 $\begin{array}{l} {P_{Gi,{t^\prime }}} = {P_{Gi,t}} + \frac{{{r_{i,t}}}}{{1 - {r_{l,t}}}}{P_{Gl,t}}\\ \;\;\; + \frac{{{k_{i,t}}}}{{1 - {k_{l,t}}}}\left( {{P_{L,t}} - {P_{W,t}} - \sum\limits_{\begin{array}{*{20}{c}} {j = 1}\\ {j \ne l} \end{array}}^{{N_G}} {{P_{Gj,t}}} } \right) \end{array}$ (10)

where the $l$th generator is out-of-service, and the lost generation $P_{Gl, t}$ will be distributed among all the other generators according to the distribution vector $\pmb{r}_{t}$. Similar as $\pmb{k}_{t}$, the $\pmb{r}_{t}$ is also designed to satisfy the following condition:

 $$$r_{1, t} + r_{2, t} +\cdots + r_{N_{\rm G}, t} = 1, \quad 0 \le r_{i, t} \le 1.$$$ (11)

In current ancillary service market, we assume the transmission system operator (TSO) purchases the AGC reserve and contingency reserve according to the distribution vector $\pmb{k}_{t}$ and $\pmb{r}_{t}$ for each generator. Therefore, in each time period $t$, the following reserves constraints are imposed:

 $$$\begin{cases} {U_{Gi, t} = k_{i, t} U_t } \\ {D_{Gi, t} = k_{i, t} D_t } \\ {R_{Gi, t} = r_{i, t} R_t } \end{cases}$$$ (12)

where $U_{t}$, $D_{t}$ and $R_{t}$ are the total AGC ramping up, ramping down, and contingency reserves at time period $t$, respectively.

C. Dynamic Economic Dispatch Model

In this paper, we consider day-ahead economic dispatch with optimal reserve scheduling, therefore the optimization horizon is set as $T=24$ with hourly steps. The objective function is formulated as:

 $$$\begin{split} \min&\sum_{t = 1}^{N_{\rm T} } {\sum_{i = 1}^{N_{\rm G} } {[a_i (P_{{\rm G}i, t} )^2 + b_i P_{{\rm G}i, t} + c_i]} } \\ & + \sum_{t = 1}^{N_{\rm T} } {\sum_{i = 1}^{N_{\rm G} } {[\alpha _i U_{Gi, t} + \beta _i D_{Gi, t} + \gamma _i R_{Gi, t}]} } \end{split}$$$ (13)

where $a_{i}$, $b_{i}$ and $c_{i}$ are the generation cost coefficients of unit $i$, $\alpha_{i}$, $\beta_{i}$ and $\gamma_{i}$ are the reserve cost coefficients of unit $i$, and $N_{T}$ is the total hourly steps. As can be observed from this objective function, it considers the commonly used generation cost in the first part (i.e., the part with coefficients $a_{i}$, $b_{i}$ and $c_{i}$), and also tries to minimize the reserve scheduling cost in the second part (i.e., the part with coefficients $\alpha_{i}$, $\beta_{i}$ and $\gamma_{i}$). To simplify the problem, we assume all the AGC generators are participating in the reserve adjustment.

The operational constraints considered in this paper are as follows:

1) System Power Balance:

 $$$\sum\limits_{i = 1}^{N_G } {P_{Gi, t} + P_{Wf, t} = P_{Lf, t} }.$$$ (14)

2) Regular Generation Limits:

 $$$P_{Gi, \min } \le P_{Gi, t} \le P_{Gi, \max }.$$$ (15)

3) Regular Generation Ramping Limits:

 $$$- r_{di} T_{di} \le P_{Gi, t} - P_{Gi, t - 1} \le r_{ui} T_{ui}$$$ (16)

where $T_{di}$ and $T_{ui}$ are the minimum down and up time of unit $i$, respectively. In this paper, these two numbers are set as the same as $T_{di} = T_{ui} = 60$ min.

4) System Reserves Constraints: the system reserves (i.e., AGC and contingency reserves) limits have been introduction in Section \ref{sec:2.1} and \ref{sec:2.2} as shown in (1)-(3) and (8)-(12). Moreover, the reserves limits should also consider the $N-$ 1 contingency situation, and result in the following chance constraints:

 $$$\begin{cases} {\Pr \left\{ {\sum\limits_{\scriptstyle i = 1 \atop \scriptstyle i \ne j}^{N_G } {\left( {U_{Gi, t} + P_{Gi, t} } \right) + } P_{W, t} \ge P_{L, t} } \right\} \ge \eta } \\ {\Pr \left\{ {\sum\limits_{\scriptstyle i = 1 \atop \scriptstyle i \ne j}^{N_G } {\left( {D_{Gi, t} - P_{Gi, t} } \right) - } P_{W, t} \ge - P_{L, t} } \right\} \ge \eta } \\ {\sum\limits_{\scriptstyle i = 1 \atop \scriptstyle i \ne j}^{N_G } {R_{i, t} \ge P_{Gj, t} } } \end{cases}$$$ (17)

where the $j$th generator is out-of-service, and $\eta$ is the confidence level.

5) Network Security Constraints:

 $$$\begin{cases} {\Pr \left\{ {\pmb{A}_j \left( {P_{G, t}^{'} - P_{L, t} + P_{W, t} } \right) \le \overline P _{Linej} } \right\} \ge \eta } \\ {\Pr \left\{ {\pmb{A}_j \left( {P_{G, t}^{'} - P_{L, t} + P_{W, t} } \right) \ge - \overline P _{Linej} } \right\} \ge \eta } \end{cases}$$$ (18)

where $A_{j}$ is the injection correlation matrix for the $j$th line, $P_{G, t}^{'}$ is the generation injection vector, and $\overline P_{Linej}$ is the power flow limit for the $j$th line. The generation injection vector is shown in (8) with normal condition, and in (10) with $N-$ 1 condition.

In summary, (1)-(3) and (8)-(18) is the proposed CCSNLP economic dispatch model considering the generation and reserve cost. The stochastic source comes from the load and wind forecasting error and the model is given in (4)-(7). The variables need to be optimized are the planned generator output $\pmb{P}_{G}$, and the distribution vector $\pmb{k}$ and $\pmb{r}$ for all the generators in the whole time span.

Ⅲ. THREE-STAGE SOLUTION FRAMEWORK A. Chance Constraints Transformation

We usually transform the probabilistic chance constraints to deterministic constraints based on the following principles:

1) If the chance constraint has the following decoupled form:

 $$$\Pr \left\{ {f\left( x \right) \ge h\left( \xi \right)} \right\} \ge \eta$$$ (19)

where $x$ and $\xi$ are decoupled, then this chance constraint could be transformed to a deterministic constraint based on the quantile of the random variable as follows:

 $$$f\left( x \right) \ge Q_\eta \left( {h\left( \xi \right)} \right).$$$ (20)

Obviously, the chance constraints in (17) could be transformed based on this situation.

2) If the chance constraint has the following coupled form:

 $$$\Pr \left\{ {f\left( x \right) \ge h\left( \xi, x \right)} \right\} \ge \eta$$$ (21)

where $f$ and $h$ are coupled, then this chance constraint should be relaxed to a deterministic constraint as follows:

 $$$\left| {A_j \left( {P_{G, t}' - P_{L, t} + P_{W, t} } \right)} \right| \le \overline P _{Linej}$$$ (22)

and the solution based on this relaxation need to be verified by using Monte Carlo simulation (MCS). Obviously, the network security chance constraints in (18) fit the form in (21). Because the transmission line power flow limit $\overline P_{Linej}$ is adjustable in this paper, which will affect the economic dispatch result as well as the injection correlation matrix $A_{j}$.

To calculate the quantile of the random variable $h(\xi)$, the Bootstrap sampling technique [31] is used in this paper, and the procedure is as follows:

Step 1: Set $k=1$;

Step 2: Carry out the $k$th sampling and obtain $N$ samples from the load and wind profile as $\vec \xi = \left[{\xi _{1, k}, \xi _{2, k}, \ldots, \xi _{N, k} } \right]$ in ascending order;

Step 3: Then calculate $s_k = \left( {\xi _{\left\lfloor {\beta N} \right\rfloor, k} + \xi _{\left\lceil {\beta N} \right\rceil, k} } \right)/2$ as the estimated value for this sampling;

Step 4: Calculate the mean value $\bar s_{k} = \frac{1}{k}\sum_{i = 1}^k {s{}_i}$ and the corresponding deviation $\sigma^{2}_{k}$;

Step 5: If the deviation $\sigma^{2}_{k} \le \sigma^{2}$, where the $\sigma$ is a pre-defined small value, then output $\bar s_{k}$, and end the sampling;

Step 6: Otherwise, $k=k+1$, go to Step 2.

In the literature, there are also many other methods could be used to find the quantile, such as the kernel density estimation based method [32] and the convolution based method [33]. However, we should notice that, the function $h(\xi)$ in this paper is mixed with load demand variable and the wind power variable. The load demand variable is a continuous random variable, while the wind power variable is a mixed random variable. The adopted Bootstrap sampling technique is easy to implement and also shows smaller error.

Based on the aforementioned technique, the original CCSNLP is transformed to a NLP. Then a three-stage solution framework is proposed and the flow chart is shown in Fig. 3. The three stages are PSO optimization stage to search the line power flow limits $\overline P_{Linej}$, the solution stage to solve the NLP by using SQP, and the verification stage to randomly verify the solution by using MCS.

 Download: larger image Fig. 3 Flow chart of the proposed three-stage solution framework.
B. PSO Optimization Stage

Before we present the detailed design of the PSO stage, we will first introduce this optimization method briefly. PSO belongs to the nature-inspired CI firstly brought forward by Kennedy and Eberhart [34]. It represents the swarm's behavior by mimicking the moving organisms in a fish school or a bird flock for food or resource searching. The main idea is to judge the adaptability of each particle in each generation based on the fitness function. The global and local best particles are selected as the examples to guide the direction and velocity updating of the rest particles. The searching capability of this method is highly related to the swarm size, generation number and the fitness function design [35], [36]. The procedure can be generalized in the following four steps.

Step 1: Initialization.

The position dimension of each particle in the swarm corresponding to the number of the transmission line limits in the whole time span $T$, denoted by $N_{\rm Line}$ in this paper. The original speed value of each particle is formed by experience. The particle number is set as 10, the iteration times (i.e., generation) is set as 20, and the searching range of each particle is set as [90%$\sim$100%] line rated power flow limits. The velocities for the position updating are initialized as follows [37]:

 $$$\begin{cases} {v_{i, \max } = 0.1(x_{i, \max } - x_{i, \min } )} \\ {v_{i, \min } = - v_{i, \max }, \quad i = 1, 2, \ldots, N_{\rm Line}\times T} \end{cases}$$$ (23)

where $x_{i, \max}$ and $x_{i, \min}$ are the upper and lower bounds of the $i$th particle, respectively. The initial velocities are generated randomly between $[v_{i, \min} \sim v_{i, \max}]$, and the other parameters, such as the accelerating constants and the initial and final learning rates, are set as: $c_{1}=c_{2}=$ 2.05, $w_{i}=$ 0.9, and $w_{f}=$ 0.4.

Step 2: Updating the local and global best.

After the PSO initialization, the particles are send to NLP solution stage and the MCS verification stage, and then the following fitness function is feedback:

 \begin{align} \nonumber Fitness =& \sum\limits_{t = 1}^{N_{\rm T} } {\sum_{i = 1}^{N_{\rm G} } {[a_i (P_{{\rm G}i, t} )^2 + b_i P_{{\rm G}i, t} + c_i]} } \\ \nonumber &+ \sum_{t = 1}^{N_{\rm T} } {\sum_{i = 1}^{N_{\rm G} } {[\alpha _i U_{Gi, t} + \beta _i D_{Gi, t} + \gamma _i R_{Gi, t}]} } \\ & + m\sum_{t = 1}^{N_{\rm T} } {\sum_{j = 1}^{N_{\rm Line}' } {(\eta - \eta _j)} } \end{align} (24)

where $\eta_{j}$ is the calculated confidence for each particle in the the MCS verification stage, and $N_{\rm Line}'$ is the number of the lines which are not satisfied the confidence interval. As can be observed from the fitness function design, we consider to minimize the generation and reserves cost, meanwhile, we also want the MCS could satisfy the specified confidence interval $\eta$. Therefore, we add this part as a punishment in the fitness function design. In general, the coefficient $m$ is a relative large value and should be tuned according to the capacity of the system [38], [39]. The local and global best particles will be determined according to the smallest fitness value of each particle.

Step 3: Position updating.

In the third step, the velocity and position of each particle is updated as:

 $$$\left\{ \begin{split} v_{i, j} (t + 1) =\, &w(t)v_{i, j} (t)+ c_1 r_1 (x^L - x_{i, j} (t)) \\ & + c_2 r_2 (x^G - x_{i, j} (t)) \\ x_{i, j} (t + 1) =\,&x_{i, j} (t) + v_{i, j} (t + 1) \end{split}\right.$$$ (25)

where $r_{1}$ and $r_{2}$ are uniformly distributed numbers between [0$\sim$1], and $x^{L}$ and $x^{G}$ are the local and global best in the current generation.

Step 4: Determining whether to finish procedure.

The process of optimization will be finished once one of the following termination criteria is satisfied:

1) The maximum generation reached.

2) The fitness function converged to a pre-defined small value.

C. NLP Solution Stage

As can be observed from the flow chart, the NLP solution stage consists of two parts. The fist part is the initialization to find proper initial value. This is carried out by randomly giving the coefficient vector $\pmb{k}$ and $\pmb{r}$ for this model, and then employ quadratic programming (QP) to solve the degenerated problem, and to find feasible initial value for the NLP. The benefit for this part is that it could speed up the whole algorithm and facilitate the convergence.

The second part is the solution part based on the initial value provided by the first part. This is carried out by using sequential quadratic programming (SQP), which uses a serials of quadratic programming to approximate the solution for Kuhn-Tucker equations [40]. The converged solution will send to the third stage for verification.

D. MCS Verification Stage

As we have mentioned above, the chance constraints have been relaxed to deterministic constraints, therefore the Monte Carlo simulation (MCS) could be carried out to verify the solution obtained in stage two. The procedure is summarized as follows:

Step 1: Set $t=1$;

Step 2: In time period $t$, sampling the wind speed and load demand for $N$ samples, and calculate the power flow for all the samples. The confidence interval for each transmission line is calculated as $\eta_{j}=n/N$, where $n$ is the number of power flow result for line $j$ which satisfies the limit. Then for all the unsatisfied results, denoted as $N_{\rm Line}'$, the calculation ${\sum_{j = 1}^{N_{\rm Line}' } {(\eta - \eta _j)} }$ is saved;

Step 3: If $t > 24$, output the result $\sum_{t = 1}^{N_{\rm T} } {\sum_{j = 1}^{N_{\rm Line}' } {(\eta - \eta _j)} }$ and feedback to the PSO;

Step 4: Otherwise, $t=t+1$, and go to Step 2.

In this three-stage framework, the NLP solution and MCS verification stages are implemented as functions called by the PSO. We should notice that, in the literature, the NLP and verification stages are usually calculated iteratively. After the verification, the line limits will be adjusted manually (e.g., decrease the limits by a mandatory value), and then calculate the NLP again with the updated limits. The proposed PSO stage in this paper could provide a heuristic technique to facilitate the line limits adjustment, therefore speed up the algorithm convergence.

Ⅳ. NUMERICAL STUDY A. Simulation Setup

In this section, numerical simulation is carried out on the modified IEEE 30-bus system with wind power integration. The system parameters are shown in Table Ⅰ and the system structure is shown in Fig. 4. As can be observed from this figure, a wind farm with rated capacity of 300 MW is integrated to the system from bus 15. The cut-in, rated, and cut-out wind speeds are set as $v_{\rm in} =$ 3.5 m/s, $v_{r} =$ 13.5 m/s and $v_{\rm out} =$ 25 m/s, respectively. The rated power flow limits for line 4-12, 9-10, 9-11, and 12-13 are set as 176 MW, 310 MW, 350 MW and 250 MW, respectively. The load forecasting error is set as 2 % of the forecasted value, and the wind speed forecasting error is set as linearly increased from 5 % to 16.5 % in the whole time span (i.e., from time period 1 to 24). The punishment coefficient is set as $m=1e^{5}$, the confidence interval is set as $\eta =$ 0.95, and the total sampling number is set as 2500. The code is implemented in MATLAB R2014b, and performed on an Intel Core i7 Processor (3.40 GHz) with 8 GB RAM.

Table Ⅰ
PARAMETERS OF THE IEEE 30-BUS SYSTEM
 Download: larger image Fig. 4 One-line diagram of the IEEE 30-bus system with wind generation.

The forecasted load demand and wind power generation for the day-ahead dispatch are shown in Fig. 5. As can be observed from this figure, the wind power is abundant around 8:00 am and decreased to a valley at 19:00 pm. However, the pattern of the load changing is just opposite to the wind, where the demand peak is between 12:00 pm to 20:00 pm.

 Download: larger image Fig. 5 The forecasted values of the wind power and load demand.
B. Numerical Results

Based on these forecasted values, the three-stage algorithm is carried out on the benchmark system, and the convergence curves with 20 independent trials are shown in Fig. 6 and 7. Specifically, the convergence of the punishment part and generation and reserve cost part are shown in Fig. 6, and the mean fitness value is shown in Fig. 7. It is interesting to notice that as the punishment decreasing (i.e., more lines are satisfied the confidence interval in the verification stage), the generation and reserve cost increasing. This result is in consistent with our intuition that the higher the system security, the greater the operating cost.

 Download: larger image Fig. 6 The punishment part, generation and reserve cost part during the optimization.
 Download: larger image Fig. 7 The mean fitness value on 20 independent trials.

The results of generation dispatch, distribution coefficient $\pmb{k}$ and $\pmb{r}$ are shown in Tables Ⅱ-, respectively. It can be observed that, $G$1 and $G$5 are scheduled with highest output for their large generation capability. However, since their ramping up/down rates are relative small, thus they are less responsible for the AGC reserve and contingency reserve task (i.e., with smaller $\pmb{k}$ and $\pmb{r}$). Meanwhile, the dispatch results for $G$3 and $G$6 are totally different in this case. Their scheduled output are relative small in the generation dispatch, but with larger $\pmb{k}$ and $\pmb{r}$ for the reserve adjustment.

Table Ⅱ
GENERATION DISPATCH (MW)
Table Ⅲ
DISTRIBUTION COEFFICIENT k
Table Ⅳ
DISTRIBUTION COEFFICIENT r

The total AGC reserve and contingency reserve purchased by the system operator is shown in Fig. 8. To accommodate the increased load demand with changing wind power generation, the AGC ramping up reserve is increasing all the way to the end. However, the AGC ramping down reserve is decreasing during 7:00-12:00 am and 21:00-24:00 pm. The reason for this result is that the wind power is abundant during these two time periods, and the wind generator could be operated at the rated power with relative higher confidence. Moreover, we should notice that the contingency reserve has the same trend with the load demand changing. This is because once a $N-$ 1 occurs in the system with a generator disconnected from the grid, the contingency reserve should compensate this generation lost to satisfy the generation-demand balancing all the time.

 Download: larger image Fig. 8 The total AGC and contingency reserves.
Ⅴ. CONCLUSIONS

In this paper, dynamic economic dispatch model for wind power penetrated system with optimal reserve scheduling was proposed and investigated. The problem was formulated as a CCSNLP problem, and transformed into a deterministic NLP problem by using Bootstrap based sampling to find the quantile and the relaxation methodology. Then a three-stage algorithm framework, with PSO optimization, SQP solution and MCS verification, was developed to tackle this problem. Simulation results on IEEE 30-bus system with wind power penetration demonstrated the convergence characteristics and the effectiveness of the proposed three-stage algorithm.

References
 [1] M. Shahidehpour, H. Yamin, Z. Y. Li, Market Operations in Electric Power Systems: Forecasting, Scheduling, and Risk Management. New York: John Wiley & Sons, Inc., 2002 . [2] Y. F. Tang, J. Yang, J. Yan, and H. B. He, "Intelligent load frequency controller using GrADP for island smart grid with electric vehicles and renewable resources, " Neurocomputing, vol. 170, pp. 406-416, Dec. 2015. [3] J. Wang, M. Shahidehpour, and Z. Li, "Contingency-constrained reserve requirements in joint energy and ancillary services auction, " IEEE Trans. Power Syst. , vol. 24, no. 3, pp. 1457-1468, Aug. 2009. [4] F. Blaabjerg, R. Teodorescu, M. Liserre, and A. V. Timbus, "Overview of control and grid synchronization for distributed power generation systems, " IEEE Trans. Ind. Electron. , vol. 53, no. 5, pp. 1398-1409, Oct. 2006. [5] D. Gautam, V. Vittal, and T. Harbour, "Impact of increased penetration of DFIG-based wind turbine generators on transient and small signal stability of power systems, " IEEE Trans. Power Syst. vol. 24, no. 3, pp. 1426-1434, Aug. 2009. [6] Y. F. Tang, H. B. He, Z. Ni, J. Y. Wen, and X. C. Sui, "Reactive power control of grid-connected wind farm based on adaptive dynamic programming, " Neurocomputing, vol. 125, pp. 125-133, Feb. 2014. [7] Y. F. Tang, H. B. He, J. Y. Wen, and J. Liu, "Power system stability control for a wind farm based on adaptive dynamic programming, " IEEE Trans. Smart Grid, vol. 6, no. 1, pp. 166-177, Jan. 2015. [8] Y. Tang, H. He, Z. Ni, J. Wen, and T. Huang, "Adaptive modulation for DFIG and STATCOM with high-voltage direct current transmission, " IEEE Trans. Neural Netw. Learn. Syst. , vol. 27, no. 8, pp. 1762-1772, Aug. 2016. [9] C. Huang, F. X. Li, T. Ding, Z. Q. Jin, and X. Ma, "Second-order cone programming-based optimal control strategy for wind energy conversion systems over complete operating regions, " IEEE Trans. Sustain. Energy, vol. 6, no. 1, pp. 263-271, Jan. 2015. [10] J. Yang, X. Feng, Y. F. Tang, J. Yan, H. B. He, and C. Luo, "A power system optimal dispatch strategy considering the flow of carbon emissions and large consumers, " Energies, vol. 8, no. 9, pp. 9087-9106, Aug. 2015. [11] Q. L. Wei, D. R. Liu, and G. Shi, "A novel dual iterative Q-learning method for optimal battery management in smart residential environments, " IEEE Trans. Ind. Electron. , vol. 62, no. 4, pp. 2509-2518, Apr. 2015. [12] Q. L. Wei, D. R. Liu, G. Shi, and Y. Liu, "Multibattery optimal coordination control for home energy management systems via distributed iterative adaptive dynamic programming, " IEEE Trans. Ind. Electron. , vol. 62, no. 7, pp. 4203-4214, Jul. 2015. [13] Y. F. Tang, H. B. He, Z. Ni, X. N. Zong, D. B. Zhao, and X. Xu, "Fuzzy-based goal representation adaptive dynamic programming, " IEEE Trans. Fuzzy Syst. , vol. 24, no. 5, pp. 1159-1175, Oct. 2016. [14] C. Huang, F. X. Li, T. Ding, Z. Q. Jin, and X. Ma, "Second-order cone programming-based optimal control strategy for wind energy conversion systems over complete operating regions, " IEEE Trans. Sustain. Energy, vol. 6, no. 1, pp. 263-271, Jan. 2015. [15] C. Huang, F. X. Li, and Z. Q. Jin, "Maximum power point tracking strategy for large-scale wind generation systems considering wind turbine dynamics, " IEEE Trans. Ind. Electron. , vol. 62, no. 4, pp. 2530-2539, Apr. 2015. [16] H. Y. Wu, M. Shahidehpour, Z. Y. Li, and W. Tian, "Chance-constrained day-ahead scheduling in stochastic power system operation, " IEEE Trans. Power Syst. , vol. 29, no. 4, pp. 1583-1591, Jul. 2014. [17] H. P. Li, C. Z. Zang, P. Zeng, H. B. Yu, and Z. W. Li, "A stochastic programming strategy in microgrid cyber physical energy system for energy optimal operation, " IEEE/CAA J. Automat. Sinica, vol. 2, no. 3, pp. 296-303, Jul. 2015. [18] N. Zhang, C. Q. Kang, Q. Xia, Y. Ding, Y. H. Huang, R. F. Sun, J. H. Huang, J. H. Bai, "A convex model of risk-based unit commitment for day-ahead market clearing considering wind power uncertainty". IEEE Trans. Power Syst. , 2015, 30 (3) :1582–1592. DOI:10.1109/TPWRS.2014.2357816 [19] T. Ding, R. Bo, F. X. Li, and H. B. Sun, "A bi-level branch and bound method for economic dispatch with disjoint prohibited zones considering network losses, " IEEE Trans. Power Syst. , vol. 30, no. 6, pp. 2841-2855, Nov. 2015. [20] F. Bouffard, F. D. Galiana, "Stochastic security for operations planning with significant wind power generation". IEEE Trans. Power Syst. , 2008, 23 (2) :306–316. DOI:10.1109/TPWRS.2008.919318 [21] J. M. Morales, A. J. Conejo, J. Perez-Ruiz, "Economic valuation of reserves in power systems with high penetration of wind power". IEEE Trans. Power Syst. , 2009, 24 (2) :900–910. DOI:10.1109/TPWRS.2009.2016598 [22] N. G. Paterakis, O. Erdinc, A. G. Bakirtzis, J. P. S. Catalao, "Load-following reserves procurement considering flexible demand-side resources under high wind power penetration". IEEE Trans. Power Syst. , 2015, 30 (3) :1337–1350. DOI:10.1109/TPWRS.2014.2347242 [23] C. Sahin, M. Shahidehpour, and I. Erkmen, "Allocation of hourly reserve versus demand response for security-constrained scheduling of stochastic wind energy, " IEEE Trans. Sustain. Energy, vol. 4, no. 1, pp. 219-228, Jan. 2013. [24] Q. Y. Xu, N. Zhang, C. Q. Kang, Q. Xia, D. W. He, C. Liu, Y. H. Huang, L. Cheng, and J. H. Bai, "A game theoretical pricing mechanism for multi-area spinning reserve trading considering wind power uncertainty, " IEEE Trans. Power Syst. , vol. 31, no. 2, pp. 1084-1095, Mar. 2016. [25] A. Papavasiliou, S. S. Oren, and R. P. O'Neill, "Reserve requirements for wind power integration: a scenario-based stochastic programming framework, " IEEE Trans. Power Syst. , vol. 26, no. 4, pp. 2197-2206, Nov. 2011. [26] M. Vrakopoulou, K. Margellos, J. Lygeros, and G. Andersson, "A probabilistic framework for reserve scheduling and N-1 security assessment of systems with high wind power penetration, " IEEE Trans. Power Syst. , vol. 28, no. 4, pp. 3885-3896, Nov. 2013. [27] H. B. He and E. A. Garcia, "Learning from imbalanced data, " IEEE Trans. Knowl. Data Eng. , vol. 21, no. 9, pp. 1263-1284, Sep. 2009. [28] P. Kundur, Power System Stability and Control. New York, USA: McGraw-Hill, 1994 . [29] D. S. Kirschen, G. Strbac, Fundamentals of Power System Economics. Chichester, UK: Wiley, 2004 . [30] C. Luo, J. Yang, Y. Z. Sun, F. Lin, and M. J. Cui, "Dynamic economic dispatch of wind integrated power system considering optimal scheduling of reserve capacity, " Proc. CSEE, vol. 34, no. 34, pp. 6109-6118, Dec. 2014. [31] X. Chen, Z. Y. Dong, K. Meng, Y. Xu, K. P. Wong, and H. W. Ngan, "Electricity price forecasting with extreme learning machine and bootstrapping, " IEEE Trans. Power Syst. , vol. 27, no. 4, pp. 2055-2062, Nov. 2012. [32] J. Shi, W. J. Lee, Y. Q. Liu, Y. P. Yang, and P. Wang, "Forecasting power output of photovoltaic systems based on weather classification and support vector machines, " IEEE Trans. Ind. Appl. vol. 48, no. 3, pp. 1064-1069, May-Jun. 2012. [33] H. Wen, Z. S. Teng, Y. Wang, and X. G. Hu, "Spectral correction approach based on desirable sidelobe window for harmonic analysis of industrial power system, " IEEE Trans. Ind. Electron. , vol. 60, no. 3, pp. 1001-1010, Mar. 2013. [34] J. Kennedy and R. Eberhart, "Particle swarm optimization, " in Proc. IEEE Int. Conf. Neural Networks, Perth, WA, Australia, 1995, pp. 1942-1948. [35] Y. Tang, H. He, and J. Wen, "Optimized control of DFIG based wind generation using swarm intelligence, " in Proc. 2013 IEEE Power and Energy Society General Meeting (PES), Vancouver, BC, Canada, 2013, pp. 1-5. [36] H. B. Duan, P. Li, and Y. X. Yu, "A predator-prey particle swarm optimization approach to multiple UCAV air combat modeled by dynamic game theory, " IEEE/CAA J. Automat. Sinica, vol. 2, no. 1, pp. 11-18, Jan. 2015. [37] Y. F. Tang, P. Ju, H. B. He, C. Qin, and F. Wu, "Optimized control of DFIG-based wind generation using sensitivity analysis and particle swarm optimization, " IEEE Trans. Smart Grid, vol. 4, no. 1, pp. 509-520, Mar. 2013. [38] X. C. Sui, Y. F. Tang, H. B. He, and J. Y. Wen, "Energy-storage-based low-frequency oscillation damping control using particle swarm optimization and heuristic dynamic programming, " IEEE Trans. Power Syst. , vol. 29, no. 5, pp. 2539-2548, Sep. 2014. [39] G. K. Venayagamoorthy, R. K. Sharma, P. K. Gautam, and A. Ahmadi, "Dynamic energy management system for a smart microgrid, " IEEE Trans. Neural Netw. Learn. Syst. , vol. 27, no. 8, pp. 1643-1656, Aug. 2016. [40] D. P. Bertsekas, Nonlinear Programming (Second edition). Belmount, Mass: Athena Scientific, 1999 .