Spinning Reserve Capacity Optimization of a Power System When Considering Wind Speed Correlation

Usually, the optimal spinning reserve is studied by considering the balance between the economy and reliability of a power system. However, the uncertainties from the errors of load and wind power output forecasting have seldom been considered. In this paper, the optimal spinning reserve capacity of a power grid considering the wind speed correlation is investigated by Nataf transformation. According to the cost–benefit analysis method, the objective function for describing the optimal spinning reserve capacity is established, which considers the power cost, reserve cost, and expected cost of power outages. The model was solved by the quantum-behaved particle swarm optimization (QPSO) algorithm, based on stochastic simulation. Furthermore, the impact of the related factors on the optimal spinning reserve capacity is analyzed by a test system. From the simulation results, the model and algorithm are proved to be feasible. The method provided in this paper offers a useful tool for the dispatcher when increasing wind energy is integrated into power systems.


Introduction
Wind power, which is a green, clean, and renewable energy source, has developed dramatically.Unlike traditional power sources, wind energy in a power system has the following characteristics: (1) the intermittent and variable nature of wind speed causes the output power of wind farms to be stochastic; (2) as wind farms are commonly clustered in a region rich in wind resources, wind speeds and wind speed forecast errors of different wind farms are dependent.As wind power cannot be predicted with great accuracy, additional spinning reserve needs to be carried in order to guarantee the operational reliability [1].Therefore, it is of great significance that new spinning reserve dispatch methods capable of taking account into probabilistic and correlated characteristics of wind power are developed.
There have been many studies on the spinning reserve of power grids with large-scale wind power.Zhang, et al. [2] and Zhang, et al. [3] set up deterministic optimization modes for spinning reserve capacity, where the objective function is the minimum cost to buy spinning reserve.In [4], the system states are quickly selected by using the enumeration and simulation method, and a new method is used to assess the operating reserve risk of the wind power system.An optimization algorithm for the spin reserve is presented in [5] for the different demand of spin reserves under risk, using the control performance of the generator set.However, the relationship between the abundant levels of spinning reserves and the reliability and economy of the power grid is not well considered.In [6], a model of the generation and optimal coordinating dispatch reserves is built, in which the constraint condition of the quantitative relation of the ratio between the system reserves and expected loss of load is considered.Taking both system economy and reliability into consideration, Literature [7] propose a conditional value-at-risk-based optimal spinning reserve model and incorporate it into the generation scheduling model in wind power integrated systems, in order to minimize total generation cost.Considering uncertainties in wind power and load forecasts, a multi-objective energy and spinning reserve scheduling method for a wind-thermal power system is investigated in a market environment [8].Cobos, et al. [9] establishes a mathematical optimization model for the minimum cost of conventional generators and spinning reserve, and then adopts a robust linear optimization to get the optimal reserve capacity with the wind uncertainty, so as to guarantee the feasibility of the optional solutions for all the realizations of the uncertain data.In [10], based on a robust optimization approach, an energy and reserve joint dispatch model in the real-time electricity markets is presented, considering wind power generation uncertainties as well as zonal reserve constraints under both normal and N-1 contingency conditions.However, the above references have not considered the influence of wind speed correlation on the spinning reserve capacity.Li, et al. [11], Chen, et al. [12], and Xie and Xiong [13] take into account the wind speed correlation to investigate optimal power flow, reliability assessment, and dynamic economic dispatch of power systems with integrated wind power, respectively.These studies have shown wind speed correlation should not be neglected when several wind farms are connected to the same power system.
In this paper, a stochastic spinning reserve capacity optimization model is set up, taking account into the wind speed correlation and some uncertain factors.Compared with the prior works, the major contributions of this paper are summarized as: (1) The wind speed correlation is considered.The Nataf transformation and Cholesky decomposition are applied to model the correlated wind speed; (2) An optimal spinning reserve model is proposed, aiming at the minimum generator cost, spinning reserve cost, and expected outage cost, in which the other uncertain factors, such as the load forecast deviation, wind power output prediction error, and forced outage rate of the generator are all considered; (3) The model is solved by the quantum-behaved particle swarm optimization (QPSO), based on the stochastic simulation algorithm.
The rest of the paper is organized as follows.Section 2 addresses the modeling method of wind speed correlation.Section 3 presents the optimization model of spinning reserve capacity.Section 4 employs the stochastic optimization algorithm, QPSO, to solve the model.In Section 5, tests and comparisons under different wind speed correlations are demonstrated.Finally, conclusions are drawn in Section 6.

Nataf Transformation
The Nataf transformation and Cholesky decomposition are used to convert non-normal relevant variables to independent standard normal ones.
For the input speed wind vector V = [v 1 , v 2 , . . ., v m ] T , whose probability density function (PDF) f i (v i ) and marginal cumulative distribution function (CDF) F i (v i ) are known, using marginal transformations, Liu and Der [14] obtained the jointly standard normal variate vector X = [x 1 , x 2 , . . ., x m ] T , which is expressed as: where Φ(•) is the CDF of standard normal variable (SNV), and Φ −1 (•) is the inverse CDF of SNV.
The relationship between correlation coefficient ρ ij of vector V and correlation coefficient ρ 0,ij of vector X can be expressed as where ρ ij of ρ 0,ij of are, respectively, the elements of correlation matrix ρ of vector V and the ρ 0 of vector X, µ i and σ i are the mean and standard deviation of wind speed v i , respectively; φ 2 (x i , x j , ρ 0,ij ) is the bi-dimensional standard normal PDF of zero means, unit standard deviations, and correlation coefficient ρ 0,ij .
In most engineering applications, ρ 0 is a positive definite, so it can be decomposed by Cholesky decomposition.
where L 0 is an inferior triangular matrix.Then, the correlated standard norm vector X is transformed into an independent standard normal vector U = [u 1 , u 2 , . . ., u m ] by using L 0 .
Through Equations ( 1)-( 4) above, the correlated non-normal wind speed vector V is transformed to the independent standard normal vector U, and this is the positive process of Nataf transformation.

The Solution for the Correlation Coefficient of the Wind Speed
For the complexity of calculating Equation (2), Dagang [15] employed the integral space transformation method to transform Equation (2) into Formula (5), as follows: where u i , u j is the ith and jth element of the vector U, respectively, and φ(•) is the PDF.
Then, the bi-dimensional Nataf transform and Gauss-Hermite integral method are applied to calculate the integral in Formula (5), ρ 0,ij can be achieved by solving the nonlinear equation as follows: where m is the number of integral nodes; w l and w k are the weight values of the integral nodes l and k, respectively; and (v il , v jk ) T is derived from Equations ( 7) and ( 8), as follows: where L 0 is a bi-dimensional inferior triangular decomposition matrix, which is derived from Nataf transformation; (z il , z jk ) T = √ 2(u il , u jk ) T ; and the typical weight values w l and w k and the node values z il and z jk can be obtained from the literature [16].

Generation of the Random Numbers of Correlated Wind Speeds
When CDFs F(V) and correlation matrix ρ of correlated wind speeds are available, the random numbers of wind speeds can be generated by inverse Nataf transformation.The basic steps are as follows: (1) Generate the random numbers U s of vector U; (2) Achieve the correlation matrix ρ 0 of the correlated standard norm vector X by Equations ( 6)-( 8), then ρ 0 is decomposed by Cholesky decomposition to get the matrix L 0 ; (3) Generate the random numbers X s of vector X by Equation ( 4); (4) Generate the correlated random numbers V s of wind speed vector V by marginal transformation.

Objective Function
According to the cost-benefit analysis method, an objective function is built as the following, which is the minimum cost containing the generation cost, reserve cost, and expected outage cost.
where N is the number of thermal power units; T is the number of the scheduling interval; f (p i,j ) is the cost of thermal power generator unit i at the time t, which includes the operating cost and start-up fee S i,t ; a i , b i , and c i are the cost coefficients of the generator unit i; P i,t is the output of thermal power unit i at the time t; g(r i,j )is the reserve cost; α i,u is the up spinning reserve price of the unit; α i,d is the down spinning reserve price of the unit i; r u,i,t and r d,i,t are the up and down spinning reserve capacities, respectively; O t is the expected outage cost; γ is the power loss cost, which can be obtained from statistical results; E t is the expected energy not serve (EENS) of the power system at the time t, which can be calculated by where y is the electricity demand deviation of the whole power system, M r is the up spinning reserve capacity at the time t, ∆P g is the dispatching difference of the thermal power units due to outage, and f (y) is PDF of the electricity demand deviation y.

Constraints
Constraint conditions mainly include the power balance (Equation ( 15)), the generation limits (Equation ( 16)), minimum running time (Equation ( 17)), outage time constraints (Equation ( 18)), and ramping constraints (Equation ( 19)), which are defined as follows: where S is the number of wind farms, P w j,t is the power output of the jth wind farm at the time t; P L,t is the load power at t; P min i,t and P max i,t are the minimum and maximum available power output of thermal power unit i at t, respectively; T on i,t and T off i,t are the operation and shutdown time, respectively, of the unit i at t; T on i,min and T off i,min are the minimum values of operation time and shutdown time, respectively; R i,Damp and R i,Uamp are the ramp-down and the ramp-up limit of the thermal power generator i, respectively.
Due to the random and fluctuation of the wind power outputs and loads, some additional constraint conditions should be considered.In this paper, the probability intervals of the wind power and loads are both set as 95%, which are respectively defined as The up and down spinning reserve constraints of the thermal power unit are as follows: where r i,u and r i,d are the up and down spinning reserve capacity supplied by the thermal power generator unit i, respectively; β is the confidence level of the system; d i,t is the state variable of the thermal power generator unit i, and d i,t = 0 when the unit i is not scheduled to run in the day-ahead dispatch.Otherwise, d i,t is determined by the forced outage rate q i of the unit i by Monte Carlo stochastic simulation: randomly generate a pseudo-random number δ i following the uniform distribution [0,1], if δ i ≤ q i , d i,t = 0, otherwise, d i,t = 1.

The Prediction Deviation
In a power system that includes wind farms, the prediction deviation is mainly from the load forecast error and the prediction error of the wind farm output.Firstly, the load forecast deviation is as follows: where P L,t is the forecast load and ∆P L,t is the load forecast deviation, which follows the normal distribution N(0, σ 2 L ).The wind farm output prediction deviation is where P W,t and P W,t are the actual and forecast values of wind farm output, respectively, and ∆P W,t is the wind power output prediction deviation, which follows the normal distribution N(0, σ 2 w ).The standard prediction deviation at the time t can be calculated as where K is the factor of wind power prediction error and W I is the total installed capacity of the wind farms.
The actual electric demand of the power system can be calculated as where P Q,t is the electric demand forecast of the whole power system, and ∆P Q,t is the electric demand deviation of the whole power system, which is assumed to follow the normal distribution N(0, σ 2 Q ).The standard deviation σ Q,t at the time t is According to the Wang, et al. [17], the continuous probability density function of the electric demand deviation is defined as

Solution Algorithm
The QPSO has many advantages [18], such as global convergence, faster convergence speed, fewer control parameters, and powerful search abilitities, so the QPSO based on the stochastic simulation algorithm is applied to solve the optimization model in this paper.

Stochastic Simulation
The procedure of the stochastic simulation is as follows: (1) Set the counter to N' = 0; (2) Randomly generate the variable samples of the thermal power output and reserve capacity, then substitute them into Formulas (23) and (24) to test the feasibility of these samples-if they satisfy (23) and (24), then N = N + 1; (3) Repeat Step 2 above for N times, until N /N > β.

Quantum-Behaved Particle Swarm Optimization
The QPSO is a probability search algorithm that introduces quantum mechanics into particle swarm optimization (PSO).The average best position C(k) is introduced into QPSO to calculate variables in the following iteration: where M is the number of the swarm, k is the current number of iterations, and p i is the local best position of the swarm i.
The particle position of decision variables x i is updated through the Monte Carlo stochastic simulation, as follows: where u is uniform random numbers in the interval [0,1], and λ is the contraction-expansion coefficient, which can be applied to control the convergence speed of PSO algorithm.The contraction-expansion coefficient is defined as where M axiter is the maximum number of the inner iterations.The flow chart of the proposed QPSO based on Equations (31)-( 33) is shown as Figure 1.The detailed procedure is as follows: (1) Obtain the correlated wind speed data of wind farms, according to Section 2.3; (2) Read the system data, such as the prediction values of the loads and wind power output and the probability distribution of the prediction deviation.In addition, set the input the parameters of the quantum-behaved particle swarm optimization algorithm, such as maximum iteration number and particle swarm size; (3) Initialize the population.The active power and reserve capacity of each thermal power unit are randomly generated to form the population, which is tested according to Formulas ( 15)-(24); if the population is not feasible, it will be regenerated until Formulas ( 15)-( 24) are satisfied; (4) Calculate the average best position of the particle swarm, according to Formula (31), and then calculate the fitness function value of particles at the current location according to Formula (10); (5) Update the position of the particles.The position of each particle is updated according to Formula (32), and the limit is verified.If the decision variables are beyond their limits, the particles are renewed.Besides, a random simulation is also employed to verify whether the particles are satisfied with the predetermined confidence level, and if they do not satisfy the confidence level, the particles will be renewed.Then the fitness function values of each particle are calculated.
If they are superior to the extreme value of the current particles, the individual extremum is updated.If the individual extreme value of the population is better than the current global extreme value, then update the global optimum; (6) Determine whether the convergence condition |F(k + 1) − F(k)| ≤ ε is satisfied, or if the number of iterations is reached.If it is not satisfied, then go back to Step 4-otherwise, output the best particle as the optimal solution.

Parameter Configuration
In this paper, the test system consists of 16 thermal power units and four wind farms.The basic data are listed in Tables A1 and A2 in Appendix A. The relationship between the wind farm power and the wind speed is as follows: where P r is the rated power of a wind farm; υ ci , υ r , and υ out are the cut-in speed, rated speed, and cut-out speed, which are 3, 12, and 20 m/s, respectively for all wind farms.Load forecast deviation obeys the normal distribution N(0, 20 2 ).The wind power prediction error factor is 0.02.The loss value for outage is 500 USD/MW•h.In QPSO, the number of swarms is 40, and the maximum number of iterations is 500, the parameter ε= 1 × 10 −4 .

Nataf Transformation
The wind speed distribution for every wind farm is assumed to follow the Weibull distribution, in which the scale and shape parameters are 8 m/s and 2.2, respectively.All four wind farms are correlated with correlation coefficient ρ, and the correlation coefficient matrix is defined as For the strong positive correlation, moderate correlation, low correlation, and negative correlation, the correlation coefficients are set as 0.9, 0.5, 0.1, and −0.5, respectively.In this paper, five integral nodes are selected to solve the correlation coefficient ρ 0 , which are listed it in Table 1.
Table 1.The correlation coefficients of wind speed.

Correlation Coefficient ρ
Correlation Coefficient ρ 0 0.9 0.9230 0.5 0.5125 0.1 0.1022 −0.5 −0.5118 Then the method proposed in Section 2.3 is applied to generate the random numbers of the correlated wind speeds, in order to get the wind speed time series of all wind farms, which are shown in Appendix B, Figure A1, and subsequently transformed them into power production to obtain the random numbers of wind power.

Optimal Spinning Reserve Capacity with Different Wind Speed Correlation
The installed capacity of the wind generator is 180 MW.The confidence level β is 0.9.The spinning reserve capacities of different wind speed correlations from 15:00 to 24:00 in one day are shown in Table 2. From Table 2, it can be seen that the up and down spinning reserve capacities rise with the increase of wind speed correlations.Wind speed with positive correlation will strengthen the synchronization (simultaneous increase and decrease) of different wind turbines' power output and increase the fluctuation of total wind power output.Therefore, compared with no correlation, the power system needs more up and down spinning reserve capacities.However, for the wind farms with negative wind speed correlation, their outputs can be complementary, so that the wind power output of the whole system becomes smooth.Therefore, the up and down spinning reserve capacity is smaller than in the case without considering wind speed correlation.Hence, the wind speed correlation has an important impact on the spinning reserve capacity of the power system with high wind power penetration, so it should not be ignored.

The Impact of Wind Speed Correlation on Spinning Reserve under Different Wind Power Capacities
Figure 2a-d shows the up and down spinning reserve capacities with different wind speed correlations, in which the wind power installed capacities are 180 and 360 MW, and the confidence level β is 0.9.
Figure 2 illustrates that the wind speed correlation has a greater influence on the spinning reserve capacity with the increase of the wind power installed capacity.In particular, the wind power forecast error will increase with time, and the wind output with positive correlation is more fluctuating; therefore, the difference in the spinning reserve capacity between the positive and the negative correlation becomes greater with time.Meanwhile, the up spinning reserve is more easily affected by the wind speed correlation than the down spinning reserve.Therefore, it is necessary to consider the influence of wind speed correlation on spinning reserve in the power system with the large-scale wind power integration.

The Effect of Wind Speed Correlation on Expected Energy Not Served
Table 3 lists the influence of the wind speed correlation on EENS, when the confidence level β is 0.9.The wind power installed capacities are 90, 180, 270, and 360 MW, respectively.As shown in Table 3, with the increase of the wind speed correlation, the EENS increases.When the wind speed correlation goes up, the wind farm output is more volatile, especially in the case of a high positive correlation; in that situation, the prediction deviation of wind power becomes larger, leading to a greater standard deviation of electricity demand in the whole power system, which gives rise to the increase of EENS.Meanwhile, EENS goes up as the wind power capacity rises.

Optimization Results under Different Confidence Levels
As shown above in Table 4, where the wind power installed capacity is 180 MW and the wind speed correlation is moderate, with an increasing confidence level the EENS decreases; however, the spinning reserve capacity and the total cost rise.Therefore, it is important to choose an appropriate confidence level to provide the trade-off between security and economy.

Conclusions
This paper proposes an optimal model of spinning reserve capacity with wind speed correlation, in which the load forecast error, wind power output prediction error, and unit forced outage rate of the uncertainty factors are considered.At the same time, the model is solved by the QPSO algorithm based on stochastic simulation.The main conclusions of this paper are as follows: • With an increase of wind speed correlation, the fluctuation of the wind farm output is greater, so the spinning reserve capacity should be increased.However, when the wind speed correlation is negative, the outputs of wind farms can be complementary, so the spinning reserve capacity should be reduced in comparison to non-correlation;

•
With the increase of wind power installed capacity, the wind speed correlation has a greater effect on the spinning reserve capacity;

•
The relationship between the total cost and confidence level of the test system was analyzed so that the results can provide decision support for dispatchers in the balance between reliability and economy.
Author Contributions: This paper is a collaborative effort among the authors.J.Z. conceived the main ideas presented in the paper and conducted the main analysis; H.Z. performed simulations and wrote part of the first draft.L.Z. and J.G. provided an in-depth review and revised the final version thoroughly together.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A

Figure 2 .
Figure 2. Spinning reserve comparisons under wind speed correlation: (a) up spinning reserve capacity with wind speed correlation under a wind power installed capacity of 180 MW; (b) down spinning reserve capacity with wind speed correlation under a wind power installed capacity of 180 MW; (c) up spinning reserve capacity with wind speed correlation under a wind power installed capacity of 360 MW; (d) down spinning reserve capacity with wind speed correlation under wind power installed capacity of 360 MW.

Funding:
This research was funded by Department of Science and Technology of Sichuan Province, China (2017GZ0263) and Department of Science and Technology of Cheng Du Province, China (2016-HM01-00275-SF).

Table 2 .
Results of the spinning reserve under different wind speed correlations.

Table 3 .
Effect of wind speed correlation on expected energy not served (EENS).

Table 4 .
Optimal dispatch results with different confidence levels.
Time/h Load/MW Time/h Load/MW Time/h Load/MW Time/h Load/MW Time/h Load/MW Time/h Load/MW

Table A2 .
Data of thermal power generators.