A Coordinated Dispatching Model Considering Generation and Operation Reserve in Wind Power-Photovoltaic-Pumped Storage System

: Large-scale grid integration of renewable energy increases the uncertainty and volatility of power systems, which brings di ﬃ culties to output planning and reserve decision-making of power system units. In this paper, we innovatively combined the non-parametric kernel density estimation method and scenario method to describe the uncertainty of renewable energy outputs, and obtained a representative set of renewable energy output scenarios. In addition, we proposed a new method to determine the reserve capacity demand. Further, we derived the quantitative relationship between the reserve demand and the power system reliability index, which was used as the constraint condition of a day-ahead power generation dispatch. Finally, a coordinated dispatching model of power generation and reserve was established, which had the lowest penalty for curtailment of wind power and photovoltaic, as well as the lowest total operating cost for thermal power units, gas power units, and pumped storage power station. By simulating three di ﬀ erent working conditions, the proposed model was compared with the traditional deterministic model. Results showed that our proposed method signiﬁcantly improved system e ﬃ ciency while maintaining system reliability.


Introduction
In recent years, wind power (WP) and photovoltaic (PV) have achieved sustained and rapid development. However, the output of WP and PV is very unstable. Large-scale WP and PV are connected to the grid and the generation side uncertainty of the power system is aggravated [1][2][3]. In Europe, balancing the power grid is mainly realized by using gas-based power generation units and pumped storage power station (PSPS), but natural gas is imported, and the geographical location of pump storage units is greatly restricted, which brings great difficulties to power grid balancing [4]. The increasing share of intermittent renewable energy generation and changing patterns of electricity demand pose challenges not only to the balance of the grid but also to the security of the supply. In China, there are many studies on renewable energy and energy storage complementary power generation technologies, but chemical energy storage systems, such as batteries, are generally preferred. The battery has its own disadvantage that cannot be changed. The service life of the battery is affected by frequent charge and discharge, which reduces the system economy. The discarded batteries may capacity or the load percentage. This method does not consider the quantitative relationship between unit output and reserve decision, and the result makes it difficult to achieve global optimization [17]. With the large-scale integration of renewable energy such as WP to the power grid, the traditional reserve decision-making methods have been difficult to apply to the power system. In order to reduce the uncertainty caused by renewable energy, some scholars proposed a coordinated dispatch model for power generation and reserve, which described the unit output and reserve decision jointly as an optimization problem with constraints, and obtained the optimal solution while satisfying the reliability constraints. Reference [18] proposed an index of expected load not supplied ratio to quantify the minimum allowable load shedding per hour and derived the quantitative relationship between this indicator and the operation reserve. However, it failed to consider the uncertainty of WP and reserve cost, so it could not meet economic requirements. Reference [19] analyzed the distribution of wind resources and studied the impact of WP grid connection and its prediction errors on power system dispatching. However, in terms of reserve decision-making, a deterministic method was adopted. The unit output and reserve decision were scheduled in sequence, without taking into account the influence of WP volatility on reserve demand; the overly simple constraint conditions made the results unrepresentative. In reference [20], a spinning reserve acquisition model was constructed based on the opportunistic constraint programming method. However, the Monte Carlo stochastic simulation may have led to long calculation time and less practicability, and the impact of positive and negative spinning reserve on the system after WP was connected to the grid was not considered. Reference [21] aimed to achieve the optimal configuration of spinning reserve with the goal of minimum the expected loss of load and the minimum unit operating cost. However, the reserve and unit output in this method were separately optimized, which made it difficult to ensure the optimal overall result. Reference [22] realized the coordination and optimization of unit output and spinning reserve under the condition of high WP penetration, but it lacked the quantification between system reliability and reserve capacity.
In the research of renewable energy output models, there are two main types. One is to establish a deterministic model to obtain the specific value of renewable energy output. This method takes into account many factors, is difficult, and the result error is large, and thus it is challenging to provide a reliable basis for the power system. The other is to establish a probability distribution model to obtain the probability density function of a renewable energy output [23,24]. Compared with the deterministic model, the latter can better explain the uncertainty of renewable energy, with higher credibility and a wider range of applications. In establishing the probability distribution model of a renewable energy output, the parameter estimation method is generally adopted, assuming that the wind speed probability obeys the Weibull distribution and the light intensity probability obeys the Beta distribution. Then, the probability distribution functions of WP and PV output are obtained according to WP unit output function and PV panel output function [25,26]. The modeling process of this method is simple, but WP output and PV output are not only affected by wind speed and light intensity, so the results have large errors. The probability model established by the non-parametric estimation method does not require model assumptions about wind speed and light intensity, and only needs to estimate the probability model based on the historical data of renewable energy output, which can effectively reveal the statistical information hidden in the historical data and reduce the influence of uncertain factors on the probability model [27]. Reference [28] adopts the autoregressive moving average method with normal distribution of WP prediction error to model the wind speed time series and obtain the probability distribution of the WP output. Then, the Monte Carlo simulation was used to generate random samples of wind speed to generate WP output scenes. This method is inefficient and time-consuming, and requires a large number of samples to get good results. In reference [29], the authors successfully deduced the analytic expression of WP density function and the fourth-order statistics based on the historical data of the WP output, and extended the model of the WP output to a regional scale. References [30][31][32] adopted the non-parametric kernel density estimation method to calculate the probability density functions of different WP prediction errors. The obtained results had higher accuracy and better adaptability than the traditional wind Energies 2020, 13, 4834 4 of 24 speed parameter distribution method, and were applied to the field of reserve capacity demand determination and power generation scheduling.
The wind-photovoltaic-pumped storage system is mainly based on the wind farms and PV power stations, according to the local favorable terrain [33]. PSPS transfers and stores the unstable and fluctuating power supply that exist in wind farms and PV power stations during the generation process, and then converts it into stable power input. According to the characteristics of the WP output in northwest China, the wind speed at night is generally greater than that during the day [34]. Therefore, when the WP output increases, the excess energy can be stored in the PSPS. The pumping operation is performed by starting the water pump and then the excess electric energy of the wind farm is converted into the gravity potential energy of water for storage. In the same way, PV power stations convert energy through the same operation. When the grid load is high and the system power supply is difficult to balance the power demand, the PSPS generates power to release energy, which can reduce load shedding loss, improve the system's utilization of new energy, and reduce the impact of power fluctuation on the grid. The structure of wind power-photovoltaic-pumped storage system is shown in Figure 1. The wind-photovoltaic-pumped storage system is mainly based on the wind farms and PV power stations, according to the local favorable terrain [33]. PSPS transfers and stores the unstable and fluctuating power supply that exist in wind farms and PV power stations during the generation process, and then converts it into stable power input. According to the characteristics of the WP output in northwest China, the wind speed at night is generally greater than that during the day [34]. Therefore, when the WP output increases, the excess energy can be stored in the PSPS. The pumping operation is performed by starting the water pump and then the excess electric energy of the wind farm is converted into the gravity potential energy of water for storage. In the same way, PV power stations convert energy through the same operation. When the grid load is high and the system power supply is difficult to balance the power demand, the PSPS generates power to release energy, which can reduce load shedding loss, improve the system's utilization of new energy, and reduce the impact of power fluctuation on the grid. The structure of wind power-photovoltaic-pumped storage system is shown in Figure 1. In summary, many scholars have used the parameter estimation method to model the WP output and PV output. The modeling process of this method is simple, but the WP and PV outputs are not only affected by wind speed and light intensity, so the results obtained by using the parameter method have a large error. The traditional reserve decision scheme of the power system is based on the deterministic model. This method performs sequential scheduling of the unit output and reserve decisions, in turn, without considering the quantitative relationship between unit output and reserve decisions. As a result, it is difficult to achieve global optimal results and the traditional reserve decision method is difficult to apply to the power system. Aiming at the above problems, the main objectives and main innovations and contributions of this paper are as follows.
The primary aims of this paper, which are also its main novelty and contributions, are to: (a) consider that, under some special circumstances, the parameter estimation method for the prediction error of renewable energy output is invalid. Thus, this paper adopts the non-parameter kernel density estimation method to model the prediction error of the renewable energy output and obtain the probability distribution function of the prediction error of renewable energy output. (b) According to the probability distribution function obtained, the Latin hypercube sampling (LHS) method was used to sample the prediction error of renewable energy output, and a large number of renewable energy output scenes were obtained. A representative set of scenes were obtained by the simultaneous backward reduction (SBR) method. (c) A method for determining the operating reserve demand capacity based on the reliability index of power system is proposed and the quantitative relationship between the up-regulated operation reserve and the expected energy not supplied (EENS) per hour is derived. Moreover, the quantitative relationship between the down-regulated operation reserve and the expected WP and PV curtailed (EWPPC) per hour is derived. (d) A coordinated optimization model of power generation and standby is established, aiming at the minimum wind and light abandoning electricity quantity and the minimum total operating cost of thermal power units, gas units, and PSPSs. We optimized the reserve supply and optimal unit allocation scheme for each hour operation by coordinated dispatching. Finally, the validity of the proposed model is verified by comparing with the deterministic model. In summary, many scholars have used the parameter estimation method to model the WP output and PV output. The modeling process of this method is simple, but the WP and PV outputs are not only affected by wind speed and light intensity, so the results obtained by using the parameter method have a large error. The traditional reserve decision scheme of the power system is based on the deterministic model. This method performs sequential scheduling of the unit output and reserve decisions, in turn, without considering the quantitative relationship between unit output and reserve decisions. As a result, it is difficult to achieve global optimal results and the traditional reserve decision method is difficult to apply to the power system. Aiming at the above problems, the main objectives and main innovations and contributions of this paper are as follows.
The primary aims of this paper, which are also its main novelty and contributions, are to: (a) consider that, under some special circumstances, the parameter estimation method for the prediction error of renewable energy output is invalid. Thus, this paper adopts the non-parameter kernel density estimation method to model the prediction error of the renewable energy output and obtain the probability distribution function of the prediction error of renewable energy output. (b) According to the probability distribution function obtained, the Latin hypercube sampling (LHS) method was used to sample the prediction error of renewable energy output, and a large number of renewable energy output scenes were obtained. A representative set of scenes were obtained by the simultaneous backward reduction (SBR) method. (c) A method for determining the operating reserve demand capacity based on the reliability index of power system is proposed and the quantitative relationship between the up-regulated operation reserve and the expected energy not supplied (EENS) per hour is derived. Moreover, the quantitative relationship between the down-regulated operation reserve and the expected WP and PV curtailed (EWPPC) per hour is derived. (d) A coordinated optimization model of power generation and standby is established, aiming at the minimum wind and light abandoning electricity quantity and the minimum total operating cost of thermal power units, gas units, and PSPSs. We optimized the reserve supply and optimal unit allocation scheme for each hour operation by coordinated dispatching. Finally, the validity of the proposed model is verified by comparing with the deterministic model.
The paper is structured as follows. In Section 2, the probability distribution model of renewable energy output ultra-short-term prediction error is mainly studied. In Section 3, the uncertainty of the WP and PV outputs are described by the scenario method. Based on the reliability index of the power system, the quantitative relationship between the hourly reserve demand and reliability index is derived in Section 4. The coordinated dispatching model considering generation and operating reserve is proposed in Section 5, and the case study is presented in Section 6. Finally, conclusions are drawn in Section 7.

The Non-Parametric Kernel Density Estimation Method
Kernel density estimation is a type of non-parametric estimation, which can describe the continuous density function well. This method does not attach any assumptions to the data distribution and only studies the distribution characteristics of the data from its own characteristics [35]. In this paper, the non-parametric kernel density estimation method is used to obtain the ultra-short-term prediction error probability distribution functions of the WP output and PV output.
The historical data of the WP output prediction error and PV output prediction error are standardized, in which P w is the WP output prediction error data; (P w1 , P w2 , . . . , P wn ) is the sample space of WP prediction error data; P pv is the PV output prediction error data; (P pv1 , P pv2 , . . . , P pvn ) is the sample space of PV output prediction error data; n is the sample size. Let the probability density function of random variable x be f (x), then the probability density function of the WP output prediction error is f w (P w ), and the probability density function of PV output prediction error is f pv (P pv ). Kernel density estimation was performed based on sample data (P w1 , P w2 , . . . , P wn ) and (P pv1 , P pv2 , . . . , P pvn ).
where h is the bandwidth; K(·) is the kernel function, the most commonly used kernel functions are the Epanechikov function and Gaussian function. This paper chooses the Gaussian kernel function because f h (x) inherits the continuity and differentiability of the kernel function, so f h (x) is differentiable of any order.
If the bandwidth h is too large, the deviation is too large and the variance too small, making f h (x) too smooth; if the bandwidth h is too small, the deviation is too small and the variance increases, resulting in f h (x) under-smoothing. The bandwidth h cannot reduce the deviation and variance at the same time, so the kernel density estimation should comprehensively weigh the deviation and the variance to find a suitable bandwidth [35].
In this paper, the empirical method is used to calculate the optimal bandwidth, so as to minimize the asymptotic mean integrated squared error (AMISE). When AMISE is the smallest, the bandwidth expression is: When selecting the kernel density estimation of the Gaussian kernel function, the normal reference criterion is used to simplify Equation (3) as: where σ is the standard deviation of the sample variable. In this paper, a more robust Interquartile range (I qr ) is considered to be adopted to replace the σ in Equation (5) with Equation (6): where Φ is the standard normal cumulative distribution function. In order to achieve an accurate estimation of the probability density distribution, the coefficient is reduced to 0.9, and the optimal bandwidth is obtained as follows:

Prediction Error Probability Distribution of the Wind Power Output and PV Power Output
The WP data used in this paper was the historical WP data for the island of Ireland [36]. The data used was the measured WP of 1 January 2014 solstice on 31 December 2014. There are 35,040 data in total and the time scale was 15 min. The PV power data used was the historical PV power data of a domestic PV power station [37]. The data used was the measured PV power of 1 December 2017 solstice on 1 September 2018. There are 17,474 data in total and the time scale was 15 min.
After processing, the prediction error kernel density estimation curves of the WP output and PV output were respectively obtained and compared with the frequency histogram of the actual data, as shown in Figure 2. It can be seen that the prediction error probability density curves of the WP output and PV output obtained by kernel density estimation were basically consistent with the frequency statistical graph of the actual data.
Energies 2020, 13, x FOR PEER REVIEW 6 of 25 When selecting the kernel density estimation of the Gaussian kernel function, the normal reference criterion is used to simplify Equation (3) as: where σ is the standard deviation of the sample variable. In this paper, a more robust Interquartile range ( qr I ) is considered to be adopted to replace the σ in Equation (5) with Equation (6): where Φ is the standard normal cumulative distribution function. In order to achieve an accurate estimation of the probability density distribution, the coefficient is reduced to 0.9, and the optimal bandwidth is obtained as follows:

Prediction Error Probability Distribution of the Wind Power Output and PV Power Output
The WP data used in this paper was the historical WP data for the island of Ireland [36]. The data used was the measured WP of 1 January 2014 solstice on 31 December 2014. There are 35,040 data in total and the time scale was 15 min. The PV power data used was the historical PV power data of a domestic PV power station [37]. The data used was the measured PV power of 1 December 2017 solstice on 1 September 2018. There are 17,474 data in total and the time scale was 15 min.
After processing, the prediction error kernel density estimation curves of the WP output and PV output were respectively obtained and compared with the frequency histogram of the actual data, as shown in Figure 2. It can be seen that the prediction error probability density curves of the WP output and PV output obtained by kernel density estimation were basically consistent with the frequency statistical graph of the actual data. We performed the following integral operation on the prediction error probability density function  We performed the following integral operation on the prediction error probability density function f w (P w ) of the WP output and the prediction error probability density function f pv (P pv ) of PV power output. where f h (x) is the prediction error probability density function of the WP output and PV power output; F h (x) is the prediction error cumulative probability distribution function of the WP output and PV power output. After calculation, the cumulative probability distribution function F w (P w ) of the WP output prediction error and the cumulative probability distribution function F pv (P pv ) of the PV output prediction error are obtained, respectively, and the cumulative probability distribution curve is shown in Figure 3.

Uncertain Scenario Set Description
The combination of the WP output scenarios and PV output scenarios increase the number of formed scene sets and increase the difficulty of calculation. Therefore, this paper used the Latin hypercube sampling method to generate a large number of initial scenarios of the WP output and initial scenarios of PV output, and then used the simultaneous backward reduction method to reduce the scenarios.

Scene Generation Based on Latin Hypercube Sampling
There are two commonly used methods for sampling: (1) Monte Carlo sampling and (2) Latin hypercube sampling. Under the same sampling scale, random variable joint covering space of the Latin hypercube sampling method was larger than that of the Monte Carlo sampling method [38]. The Latin hypercube sampling method is essentially hierarchical. As shown in Figure 4, the core idea of sampling is to divide the cumulative curve into equal intervals on the cumulative probability scale [0,1]. Then, in order to ensure the sampling point covers the random distribution area of all input random variables, it is suggested to take samples randomly from each interval of the input distribution and force the sampling to represent the value of each interval.

Uncertain Scenario Set Description
The combination of the WP output scenarios and PV output scenarios increase the number of formed scene sets and increase the difficulty of calculation. Therefore, this paper used the Latin hypercube sampling method to generate a large number of initial scenarios of the WP output and initial scenarios of PV output, and then used the simultaneous backward reduction method to reduce the scenarios.

Scene Generation Based on Latin Hypercube Sampling
There are two commonly used methods for sampling: (1) Monte Carlo sampling and (2) Latin hypercube sampling. Under the same sampling scale, random variable joint covering space of the Latin hypercube sampling method was larger than that of the Monte Carlo sampling method [38]. The Latin hypercube sampling method is essentially hierarchical. As shown in Figure 4, the core idea of sampling is to divide the cumulative curve into equal intervals on the cumulative probability scale [0,1]. Then, in order to ensure the sampling point covers the random distribution area of all input random variables, it is suggested to take samples randomly from each interval of the input distribution and force the sampling to represent the value of each interval.
The steps are as follows [39]: Step 1: Let x 1 , x 2 , . . . , x T be T independent random variables, whose cumulative probability distribution function is: Step 2: Suppose M is the sampling size and the vertical axis of the cumulative probability distribution function curve Φ t = F h (x t ) is divided into M equal intervals, non-overlapping, and with a width of 1 M , whose size is [ n−1 M , n M ], n = 1, 2, . . . , M.
Energies 2020, 13, 4834 8 of 24 Step 3: The authors select the midpoint of each interval as the sampling value of Φ t and calculate the sampling value of x t by taking the inverse function of the cumulative distribution function Φ t = F h (x t ), i.e., the m-th sampling value of x t .
All sampling values x tm form a T × M initial sampling matrix X.
Energies 2020, 13, x FOR PEER REVIEW 8 of 25 The steps are as follows [39]: Step 1：Let 1 2 , ,..., T x x x be T independent random variables, whose cumulative probability distribution function is: Step 2：Suppose M is the sampling size and the vertical axis of the cumulative probability distribution function curve Step 3：The authors select the midpoint of each interval as the sampling value of Φ t and calculate the sampling value of t x by taking the inverse function of the cumulative distribution All sampling values tm x form a × T M initial sampling matrix X .

Scenario Reduction Based on Simultaneous Backward Reduction Method
A large number of time series scenarios with the same probability were obtained by sampling methods, which described the uncertainty of the WP output and PV output more accurately. However, the combination of the WP and PV output scenarios made the number of scenes increase sharply, which created low calculation efficiency. Therefore, it was possible to use the simultaneous backward reduction method to reduce the WP and PV scenarios. Lastly, we obtained the probability of each scenario. The steps are as follows [40]: Step 1：The authors first get the large-scale scenarios X X X X and set the number of scenarios to be deleted as K .
Step 2：In this step, the authors calculate the Kantorovich distance for each pair of scenarios: T Figure 4. The structure of Latin hypercube sampling.

Scenario Reduction Based on Simultaneous Backward Reduction Method
A large number of time series scenarios with the same probability were obtained by sampling methods, which described the uncertainty of the WP output and PV output more accurately. However, the combination of the WP and PV output scenarios made the number of scenes increase sharply, which created low calculation efficiency. Therefore, it was possible to use the simultaneous backward reduction method to reduce the WP and PV scenarios. Lastly, we obtained the probability of each scenario. The steps are as follows [40]: Step 1: The authors first get the large-scale scenarios X = {X 1 , X 2 , . . . , X M } and set the number of scenarios to be deleted as K.
Step 2: In this step, the authors calculate the Kantorovich distance for each pair of scenarios: where X i is the i-th scenario, and x t,i is the t-th element of scenario i.
Step 3: For any scenario, X i , the authors compared the distance of the matching scene pair and then found scenarios closest to scenario, X i . These scenarios are deleted according to the principle of scene reduction. Finally, the authors accumulate the probability of deleted scenes to the closest scene.
Step 4: Repeat step 3 until the number of deleted scenes are reached K. Finally, the authors obtain the reduced WP output scenarios, PV output scenarios, and the corresponding scenario probability.

Reserve Capacity Quantitative Model Based on the Reliability Index
Based on the scenarios established above, the system reliability was expressed by the expected energy not supplied (EENS) and the expected WP and PV curtailed (EWPPC). Based on reliability indexes and considering the uncertain factors of power system, the quantitative relationship between the reserve demand and reliability indexes was deduced.

The Expected Energy Not Supplied
Operation reserve refers to the rapid active power response capacity reserved to meet the reliable and continuous power supply for the load. It can cope with the load and the renewable energy output fluctuation, generator outage fault, and so on [15]. If the prediction error of the renewable energy output in the system is too large or the generator is faulty, it leads to a power shortage in the system. At the same time, if the up-regulated reserve cannot meet the power shortage, then the loss of load occurs.
If the reserve capacity provided by the system at time t under scenario k is less than the power shortage, then the required shear load is: (12) where R k up,t is provided by the thermal power units, gas power units, and pumped storage unit in the normal operation.
ga,j,t = min{r up ga, j ∆t, P max ga,j − P k ga,j,t } R +,k ps,t = min{P max ps,out , (W k ps,t − W min ps )/∆t} (14) Assuming that E t 0 is the minimum load shedding capacity at time t, according to the reliability requirements, then: Equation (16) is the relationship between up-regulated reserve and reliability index of the system at time t.
j P k ga,j,t − P k ps,out,t + P k ps,in,t − R k up,t ) (16) where the minimum R k up,t is the minimum up-regulated reserve required by the system at time t under scenario k:

Expected Wind Power and PV Curtailed
After the large-scale renewable energy is connected, the system peak shaving problem is prominent. Especially for systems with cogeneration units, the peak shaving depth of the system is limited due to the operation mode in winter, which intensifies wind and PV curtailment [41]. Therefore, the reasonable arrangement of the system's operation reserve is beneficial to reduce wind and solar curtailment.
According to the above scene set, the expected WP and PV curtailment at time t under scenario k is: pv,t + ξ k pv,t + i P k th,i,t + j P k ga,j,t + P k ps,out,t − P k ps,in,t − P pre l,t − ξ k l,t − R k dn,t , 0} (18) where R k dn,t is provided by the thermal power units, gas power units, and pumped storage unit in normal operation.
ga,j,t = min{r dn ga,j ∆t, P k ga, j,t − P min ga,j } R −,k ps,t = min{P max ps,in , (W max ps − W k ps,t )/∆t} (20) Assuming that C t 0 is the minimum wind and PV curtailment at time t, according to the reliability requirements, then we can find that: Equation (22) is the relationship between the down-regulated reserve and reliability index of the system at time t.
where the minimum R k dn,t is the minimum down-regulated reserve required by the system at time t, we can find that:

The Coordination and Optimization Model of Power Generation and Reserve
In this paper, the system reserve capacity demand was divided into up-regulated reserve and down-regulated reserve, which were used to deal with load shortage and renewable energy output climbing events, respectively. Considering the coordination between reserve dispatching and generation dispatching, the system reserve plan can be determined while units were committed and the output plan was formulated. Therefore, the model obtained the system units generation plan, reserve plan, and minimum reserve demand of the system at each moment, simultaneously. Based on the scenario set, the following coordinated optimization model of power generation and reserve was established, which considered the reliability and economic efficiency of the system.

Objective Function
In order to realize the maximum absorption of the renewable energy and economic efficiency of the system's operation, the objective function of the model comprehensively considered the minimum amount of WP and PV power curtailment, as well as the lowest total operation cost of thermal power units, gas power units and PSPS.
The specific calculation formula of operating cost of various units and equipment is shown in Equation (25).
[a ga,j (P k ga,j,t ) 2 + b ga,j P k ga,j,t + c ga,j ] C k ps,t = w out P k ps,out,t U out,t + w in P k ps,in,t U in,t Energies 2020, 13, 4834 11 of 24

Operation Constraints of Thermal Power Units
The upper and lower limits of output constraints and the climb rate constraints of thermal power units are mainly considered in the model. The specific constraints are as follows: − r dn th,i ·∆t ≤ P k th,i,t+1 − P k th,i,t ≤ r up th,i ·∆t Large thermal power units are the main generators of the power system. It takes 1-2 days for the units from the cold reserve state of the boiler to the grid connection, which may take longer or longer due to the impact of scheduling instructions [42]. Therefore, the model in this paper did not consider unit commitment. Once the thermal power units are determined to operate, the intra-day change did not occur.

Other constraints
The constraints of gas power unit also includes the upper and lower limits of output constraints and the climb rate constraints, which are similar to the constraints of thermal power unit, so it will not be repeated here.

State constraints of generation/pumping
In this paper, a virtual generator and motor are used to represent the two working states of the power generation and pumping, and they can only be in one working mode at the same time.

Power constraints of generation/pumping
The power generation of PSPS can change continuously, but the pumping power is usually a constant value.
       P k ps,in,t = P max ps,in 0 ≤ P k ps,out,t ≤ P max ps,out

Reliability constraints
3. Network security constraints P d,min ≤ P k d,t ≤ P d,max

Calculating Procedures
In this paper, we established a coordinated dispatching model of the WP-PV-pumped storage system with WP, PV power, thermal power units, and gas power units. The ultra-short term prediction error probability density distributions of the WP output and PV output were analyzed by the nonparametric kernel density estimation method and the uncertainties were described by a scenario set. Considering the operation constraints of various power sources and the minimum amount of WP and PV power curtailment, as well as the lowest total operation cost, the coordinated dispatching model of power generation and reserve was established. The unit operating cost was transformed into a linear function by piecewise linearization and the reliability constraint was transformed into the inequality constraint by derivation. Therefore, the whole problem was transformed into a mixed integer linear programming problem, which was solved by using MATLAB 2016a programming and the Yalmip toolbox to call the solver CPLEX12.8.

Parameters of the Calculation Example
The examples in this paper included a 200 MW wind farm, a 40 MW PV power station, four thermal power units with a total installed capacity of 1290 MW, two gas power units with a total installed capacity of 120 MW, and a pumped storage unit with a maximum capacity of 150 MW·h, a maximum pumping power of 30 MW, and a maximum generation power of 30 MW. The specific parameters of the units are shown in Table A1.
Examples in this paper verified the practicability of the model from three working conditions. Due to the small installed capacity of PV power station, the impact caused by its uncertainty was far less than that caused by WP, so the impact of PV output uncertainty was not considered. The three working conditions were as follows: ordinary operating conditions, the working conditions with continuous large wind power output, and the working conditions with continuous a small wind power output. The day-ahead WP output forecast data under different working conditions is shown in Figure 5. The day-ahead load forecast data and the day-ahead PV output forecast data are shown in Figure 6. less than that caused by WP, so the impact of PV output uncertainty was not considered. The three working conditions were as follows: ordinary operating conditions, the working conditions with continuous large wind power output, and the working conditions with continuous a small wind power output. The day-ahead WP output forecast data under different working conditions is shown in Figure 5. The day-ahead load forecast data and the day-ahead PV output forecast data are shown in Figure 6.   less than that caused by WP, so the impact of PV output uncertainty was not considered. The three working conditions were as follows: ordinary operating conditions, the working conditions with continuous large wind power output, and the working conditions with continuous a small wind power output. The day-ahead WP output forecast data under different working conditions is shown in Figure 5. The day-ahead load forecast data and the day-ahead PV output forecast data are shown in Figure 6.

Scenario Generation and Reduction
The authors used a multi-scene model to simulate the uncertainty of the WP output and PV output. In this section, the WP output under ordinary operating conditions was used as an example. The LHS method was adopted to obtain 200 WP output scenes and PV output scenes, respectively. It can be seen from Figure 7 that the variation trend of the scene simulated by LHS method was roughly the same, except that the output in each period appeared with corresponding deviation within a certain confidence interval.
After simultaneous backward reduction, three WP output curves and three PV output curves were obtained, as shown in Figure 8. The scenario probability corresponded to the reduced scenarios, as shown in Table 1.
The authors used a multi-scene model to simulate the uncertainty of the WP output and PV output. In this section, the WP output under ordinary operating conditions was used as an example. The LHS method was adopted to obtain 200 WP output scenes and PV output scenes, respectively. It can be seen from Figure 7 that the variation trend of the scene simulated by LHS method was roughly the same, except that the output in each period appeared with corresponding deviation within a certain confidence interval. After simultaneous backward reduction, three WP output curves and three PV output curves were obtained, as shown in Figure 8. The scenario probability corresponded to the reduced scenarios, as shown in Table 1.   the same, except that the output in each period appeared with corresponding deviation within a certain confidence interval. After simultaneous backward reduction, three WP output curves and three PV output curves were obtained, as shown in Figure 8. The scenario probability corresponded to the reduced scenarios, as shown in Table 1.    Similarly, the WP output scenarios and its corresponding scenario probability can be obtained under the other working conditions.

The Dispatching Test Results under Different Working Conditions
In this section, the authors analyzed the output status of each power source, the working status of PSPS, and the reserve condition from the perspective of three working conditions. In order to verify the validity of the model in this paper, two different alternative decision making schemes were selected for comparative analysis under the same calculation example and parameters. Scheme 1 was the operation reserve capacity scheme proposed in this paper. Scheme 2 was a deterministic reserve capacity scheme, whose operating reserve is shown in the following formula.

Ordinary Operating Conditions
The output of each power source under ordinary operating conditions is shown in Figure 9. Under ordinary operating conditions, the WP had a strong anti-peak-regulating property, so its output was large at night and low in the morning. However, the night load level was low, which seriously compresses the generating space of thermal power units, so the output of thermal power units showed a downward trend. At this time, the down-regulated reserve provided by the thermal power unit was lower. If the down-regulated reserve was not supplemented, there was wind curtailment. The pumped storage unit provided a major down-regulated reserve capacity through pumping operations, converting wind curtailment electricity into potential energy and storing it. During 9:00-12:00, the first load peak occurred and a large number of up-regulated reserves were consumed. Thermal power units all operated at a high output level and, at this time, up-regulated reserve units were mainly provided by gas power units and pumped storage units. When the gas power units ran at a lower output level, they provided a large number of up-regulated reserves to cope with load surge events. Similarly, during 19:00-22:00, the second load peak occurred, and thermal power units operated at the maximum output level. At this time, gas power units participated in dispatching to provide up-regulated reserve, while pumped storage unit participates in the dispatching to reduce the load level.
the validity of the model in this paper, two different alternative decision making schemes were selected for comparative analysis under the same calculation example and parameters. Scheme 1 was the operation reserve capacity scheme proposed in this paper. Scheme 2 was a deterministic reserve capacity scheme, whose operating reserve is shown in the following formula.

Ordinary Operating Conditions
The output of each power source under ordinary operating conditions is shown in Figure 9. Under ordinary operating conditions, the WP had a strong anti-peak-regulating property, so its output was large at night and low in the morning. However, the night load level was low, which seriously compresses the generating space of thermal power units, so the output of thermal power units showed a downward trend. At this time, the down-regulated reserve provided by the thermal power unit was lower. If the down-regulated reserve was not supplemented, there was wind curtailment. The pumped storage unit provided a major down-regulated reserve capacity through pumping operations, converting wind curtailment electricity into potential energy and storing it. During 9:00-12:00, the first load peak occurred and a large number of up-regulated reserves were consumed. Thermal power units all operated at a high output level and, at this time, up-regulated reserve units were mainly provided by gas power units and pumped storage units. When the gas power units ran at a lower output level, they provided a large number of up-regulated reserves to cope with load surge events. Similarly, during 19:00-22:00, the second load peak occurred, and thermal power units operated at the maximum output level. At this time, gas power units participated in dispatching to provide up-regulated reserve, while pumped storage unit participates in the dispatching to reduce the load level.  It can be seen from Figure 10 that the PSPS increased the load level by pumping water during the period of abundant WP at night, while in peak load periods, such as the period of 9:00-12:00 and 19:00-22:00, the PSPS reduced the load level through generation. Therefore, the PSPS had a significant effect of peak load cutting, which provided great flexibility for the operation and scheduling of power system.
The comparison between the up-regulated reserve supply and demand curves and the down-regulated reserve supply and demand curves at each moment is shown in Figure 11. By comparing the results of the two kinds of reserve capacity decision-making schemes, it can be seen that the reserve capacity demand obtained by the traditional deterministic scheme was lower than that obtained by the optimized scheme in this paper. The reserve capacity supply obtained by the traditional scheme was also smaller than the result of the scheme in this paper. Because the scheme proposed in this paper met the system power balance and reserve capacity constraints under all scenarios, there were many units involved in the dispatching, which led to the opening of some small capacity units with high operating costs. Therefore, compared with the traditional scheme, the optimized scheme proposed in this paper increased the total operating cost. It can be seen from Figure 10 that the PSPS increased the load level by pumping water during the period of abundant WP at night, while in peak load periods, such as the period of 9:00-12:00 and 19:00-22:00, the PSPS reduced the load level through generation. Therefore, the PSPS had a significant effect of peak load cutting, which provided great flexibility for the operation and scheduling of power system. The comparison between the up-regulated reserve supply and demand curves and the downregulated reserve supply and demand curves at each moment is shown in Figure 11. By comparing the results of the two kinds of reserve capacity decision-making schemes, it can be seen that the reserve capacity demand obtained by the traditional deterministic scheme was lower than that obtained by the optimized scheme in this paper. The reserve capacity supply obtained by the traditional scheme was also smaller than the result of the scheme in this paper. Because the scheme proposed in this paper met the system power balance and reserve capacity constraints under all scenarios, there were many units involved in the dispatching, which led to the opening of some small capacity units with high operating costs. Therefore, compared with the traditional scheme, the optimized scheme proposed in this paper increased the total operating cost.  The comparison between the up-regulated reserve supply and demand curves and the downregulated reserve supply and demand curves at each moment is shown in Figure 11. By comparing the results of the two kinds of reserve capacity decision-making schemes, it can be seen that the reserve capacity demand obtained by the traditional deterministic scheme was lower than that obtained by the optimized scheme in this paper. The reserve capacity supply obtained by the traditional scheme was also smaller than the result of the scheme in this paper. Because the scheme proposed in this paper met the system power balance and reserve capacity constraints under all scenarios, there were many units involved in the dispatching, which led to the opening of some small capacity units with high operating costs. Therefore, compared with the traditional scheme, the optimized scheme proposed in this paper increased the total operating cost.

The Working Conditions with Continuous Large Wind Power Output
The output of each power source under the second working condition is shown in Figure 12. It can be seen from Figure 12 that the WP output was large throughout the day. Therefore, in order to absorb more WP, the generation space of thermal power units had to be compressed. However, the load at night was low and thermal power units were at the minimum output level. In order to maintain the operation of thermal power units and meet the reserve capacity demand, WP was reduced. PSPS stored energy by pumping water during 1:00-4:00, which absorbed excess WP to a certain extent. However, due to capacity and pumping power rate restrictions, it was difficult to fully absorb WP. Due to the high WP output, thermal power units were enough to bear the change of load, so there was no gas power unit scheduling throughout the day. In the second peak load, the pumped storage unit generated electricity instead of the gas power unit, which not only ensured the reliability of the power system but also reduced the generation cost of the gas power unit and ensured the economic efficiency of the system. certain extent. However, due to capacity and pumping power rate restrictions, it was difficult to fully absorb WP. Due to the high WP output, thermal power units were enough to bear the change of load, so there was no gas power unit scheduling throughout the day. In the second peak load, the pumped storage unit generated electricity instead of the gas power unit, which not only ensured the reliability of the power system but also reduced the generation cost of the gas power unit and ensured the economic efficiency of the system. The changes of pumping, power generation, and storage capacity of the PSPS under this working condition are shown in Figure 13. It can be seen that the PSPS absorbs part of the WP through pumping operation during the night wind curtailment period. However, the WP output throughout the day is relatively high, so during 14:00-17:00, the PSPS reached the maximum storage capacity, and there is no space to store excess WP, thus causing certain wind curtailment. At 24:00, the storage capacity of the pumping and storage plant reached its maximum, which may have had a negative impact on WP consumption the next day. If the WP output was still at a high level, it caused more wind curtailment. The changes of pumping, power generation, and storage capacity of the PSPS under this working condition are shown in Figure 13. It can be seen that the PSPS absorbs part of the WP through pumping operation during the night wind curtailment period. However, the WP output throughout the day is relatively high, so during 14:00-17:00, the PSPS reached the maximum storage capacity, and there is no space to store excess WP, thus causing certain wind curtailment. At 24:00, the storage capacity of the pumping and storage plant reached its maximum, which may have had a negative impact on WP consumption the next day. If the WP output was still at a high level, it caused more wind curtailment. The comparison between the up-regulated reserve supply and demand curves and the downregulated reserve supply and demand curves at each moment is shown in Figure 14. Due to the large WP output, the net load was lower, and the thermal power unit output was at a low level, so the upregulated reserve was higher, and the down-regulated reserve supply capacity was lower. High WP output increased the uncertainty of the system and the reserve demand. However, the reserve results obtained by the traditional deterministic scheme were too conservative, which may have led to the shortage of reserve in the actual operation. The comparison between the up-regulated reserve supply and demand curves and the down-regulated reserve supply and demand curves at each moment is shown in Figure 14. Due to the large WP output, the net load was lower, and the thermal power unit output was at a low level, so the up-regulated reserve was higher, and the down-regulated reserve supply capacity was lower. High WP output increased the uncertainty of the system and the reserve demand. However, the reserve results obtained by the traditional deterministic scheme were too conservative, which may have led to the shortage of reserve in the actual operation. The comparison between the up-regulated reserve supply and demand curves and the downregulated reserve supply and demand curves at each moment is shown in Figure 14. Due to the large WP output, the net load was lower, and the thermal power unit output was at a low level, so the upregulated reserve was higher, and the down-regulated reserve supply capacity was lower. High WP output increased the uncertainty of the system and the reserve demand. However, the reserve results obtained by the traditional deterministic scheme were too conservative, which may have led to the shortage of reserve in the actual operation.

The Working Conditions with Continuous Small Wind Power Output
The output of each power source under the third working condition is shown in Figure 15. In this working condition, the all-day output of WP was at a low level, which brought great pressure to the generation scheduling of thermal power units. At night, WP output was low and the thermal power units increased output to balance the load demand. During the daytime, during the two peak load periods and thermal power units were at a high output level. In order to ensure the reserve capacity, the pumped storage unit and gas power units participated in the dispatching. The pumped storage unit reduced the load level and the gas power units provided reserve capacity for the power system and increased the flexibility of the power system. WP was fully consumed by the coordinated dispatching of power generation and reserve capacity without affecting the reliability of the power system.

The Working Conditions with Continuous Small Wind Power Output
The output of each power source under the third working condition is shown in Figure 15. In this working condition, the all-day output of WP was at a low level, which brought great pressure to the generation scheduling of thermal power units. At night, WP output was low and the thermal power units increased output to balance the load demand. During the daytime, during the two peak load periods and thermal power units were at a high output level. In order to ensure the reserve capacity, the pumped storage unit and gas power units participated in the dispatching. The pumped storage unit reduced the load level and the gas power units provided reserve capacity for the power system and increased the flexibility of the power system. WP was fully consumed by the coordinated dispatching of power generation and reserve capacity without affecting the reliability of the power system. As shown in Figure 16, due to the small WP output level, although the PSPS carried out energy storage, the storage capacity did not reach the maximum value, and the PSPS was underutilized, which reduced the generation capacity available at the peak load in the daytime. 150 150 Water pumping power  As shown in Figure 16, due to the small WP output level, although the PSPS carried out energy storage, the storage capacity did not reach the maximum value, and the PSPS was underutilized, which reduced the generation capacity available at the peak load in the daytime. As shown in Figure 16, due to the small WP output level, although the PSPS carried out energy storage, the storage capacity did not reach the maximum value, and the PSPS was underutilized, which reduced the generation capacity available at the peak load in the daytime. The comparison between the up-regulated reserve supply and demand curves, and the downregulated reserve supply and demand curves at each moment is shown in Figure 17. Due to the low WP output, the net load was relatively high and the thermal power unit was output at a higher level. Therefore, the down-regulated reserve supply capacity was larger, while the up-regulated reserve supply capacity was lower. The scheme in this paper met the reserve capacity constraints in all scenarios. Although the traditional deterministic scheme added additional reserve requirements, it provided less reserve capacity, which means there may have been an insufficient reserve in the actual operation. The comparison between the up-regulated reserve supply and demand curves, and the down-regulated reserve supply and demand curves at each moment is shown in Figure 17. Due to the low WP output, the net load was relatively high and the thermal power unit was output at a higher level. Therefore, the down-regulated reserve supply capacity was larger, while the up-regulated reserve supply capacity was lower. The scheme in this paper met the reserve capacity constraints in all scenarios. Although the traditional deterministic scheme added additional reserve requirements, it provided less reserve capacity, which means there may have been an insufficient reserve in the actual operation.

Conclusions
In order to solve the dispatching problem of unit outputs and reserve capacity decisions caused by large-scale WP and PV power in the power system, this paper innovatively combined the nonparametric kernel density estimation method and scenario method to describe the uncertainty of renewable energy output, based on the multi-scenario analysis method. In addition, a new method to determine the reserve capacity demand was proposed to derive the quantitative relationship between the reserve demand and the reliability index of power system. Through example analysis, we can draw the following conclusions: • The non-parametric kernel density estimation method did not need to assume the distribution model of variables and had little limitation on the model. According to Figure 2, the modeling

Conclusions
In order to solve the dispatching problem of unit outputs and reserve capacity decisions caused by large-scale WP and PV power in the power system, this paper innovatively combined the non-parametric kernel density estimation method and scenario method to describe the uncertainty of renewable energy output, based on the multi-scenario analysis method. In addition, a new method to determine the Energies 2020, 13, 4834 20 of 24 reserve capacity demand was proposed to derive the quantitative relationship between the reserve demand and the reliability index of power system. Through example analysis, we can draw the following conclusions: • The non-parametric kernel density estimation method did not need to assume the distribution model of variables and had little limitation on the model. According to Figure 2, the modeling process of the non-parametric method was simple and was subject to little interference from external factors. Compared with the parametric method, the results obtained by the non-parametric method had a small error and was highly practical. • According to Figures 9, 12 and 15, it can be seen that the PSPS achieved peak load cutting and valley load reduction, as well as reduced the peak-valley difference of load. The gas power unit had flexible adjustment ability and provided a large amount of reserve capacity. Therefore, globally, the number of gas power units and PSPSs were appropriately increased so that they can participate in power grid dispatching, which effectively relieved peak pressure of thermal power units and further reduced renewable energy waste and load cutting accidents.

•
Compared with the traditional deterministic alternative decision method, the method proposed in this paper effectively solved its blindness. By comparing Figure 11, Figure 14, and Figure 17, we can clearly find that the method proposed in this paper met the system power balance constraint and reserve capacity constraint in all scenarios, and the results of the reserve demand and supply were higher than those of the traditional deterministic reserve decision method. Therefore, there was no shortage of reserve in the operation process. The method proposed in this paper is of great practical value for active power dispatching of power systems with large-scale renewable energy sources.
At present, the complementary power generation technology of new energy system and energy storage system is not fully mature. Moreover, due to the complexity of problems, the actual WP-PV-storage complementary power generation technology has other factors that need to be considered. The depth and breadth of the research are limited and further studies can be carried out in the follow-up work as follows: • The correlation between WP output and PV output was not considered. WP output and PV output have certain complementarity, i.e., negative correlation. In the follow-up research, different Copula functions can be introduced to form the WP-PV joint probability distribution model according to the situation, which can better fit with the engineering practice.

•
The actual power system parameters should be used for verification and more reference examples should be added, so that the analysis of the problem are more convincing and targeted.

•
The coordinated scheduling problem proposed in this paper is an unusually complex mixed integer nonlinear programming problem. When the system size is large, it is very difficult to directly optimize and solve the problem. Therefore, it is urgent to seek an optimization algorithm to solve this problem.