Exploring the Regulation Reliability of a Pumped Storage Power Plant in a Wind–Solar Hybrid Power Generation System

: In the coming decades, the proportion of wind–solar energy in power system signiﬁcantly increases, resulting to uncertainties of power ﬂuctuation in abundant wind–solar energy regions. The ﬂexibility operation of Pumped Storage Power Plants (PSPPs) has already been widely recognized to regulate wind–solar power ﬂuctuations; however, less is known about the regulation reliability of the PSPP affected by them. It is a challenge, since various uncertainties exist during this regulation process. Here, a mathematical model with a solar–wind–hydro hybrid power generation system is adopted to investigate the regulation reliability of PSPP. The uncertainties and limitations of model parameters are considered during this process. Five regulation indexes, i.e., rise time, settling time, peak value, peak time and overshoot of the reactive power generator terminal voltage, guide vane opening and angular velocity, are extracted to evaluate the PSSP’s regulation quality. Finally, the PSPP reliability probability affected by parametric uncertainties is presented. The obtained results show that the inertia coefﬁcient is the most sensitivity parameters for the settling time, peak value and peak time with sensitivity index 33.7%, 72.55% and 71.59%, respectively. The corresponding total contribution rate of the top 10 sensitive parameters are 74.45%, 93.45% and 87.15%, respectively. Despite some types of uncertainties not being considered, the results of this research are important for the regulation reliability evaluation of PSPPs in suppressing power ﬂuctuations of wind and solar generation.


Introduction
Wind-Solar-Hydro (WSH) hybrid power stations in China are under construction to keep sustainable growth in the economies, for it is a widely acknowledged fact that the use of renewable energy sources plays an important role in achieving the Paris Agreement of 2015 [1]. The flexibility characteristics of PSPPs, such as flexible adjustment [2], rapidity and reliability [3], increase the utilization and the proportion of wind-solar power in the power system [4]. However, the intermittent and random characteristics of wind-solar power easily lead to low frequency oscillation of active power and fatigue failures of the guide vane in the regulation process. Moreover, the structure of the power system connected with large scales of wind-solar farms becomes more complex and its dynamic characteristics are affected by various uncertainties, such as climate uncertainty, load demand uncertainty, etc. [5]. These uncertainties bring a great burden on peak shaving and frequency modulation, posing a challenge for regulation reliability of PSPPs in suppressing power fluctuations of wind and solar generation.
Historically, studies of this challenge have been divided into three research directions, namely, complementary uncertainty, regulation performance, and reliability. As for complementarity uncertainty, Han et al., proposed a complementarity evaluation method for the WSH system by thoroughly examining the fluctuation of the independent and combined power generation, and the results showed that the best complementarity level is obtained by changing the proportion of wind and photovoltaic power [6]. François et al., pointed out that the uncertainty of the independent power generation is reduced by integrating the hydropower to the wind/solar mix, thus improving the complementary [7]. Cantao et al., applied the hydro-wind correlation maps to the representative hydropower plants of Brazilian basin, and concluded that the method is a useful tool to analyze the complementarity uncertainty and similarity [8]. With regard to regulation performance during the complementarity process, Tang et al., studied the regulation characteristics of pumped-storage plants integrated with wind power and estimated the regulation quality using time domain simulation [9]. Martinez-Lucas et al., took an isolated island with generation 100% wind and hydro as a research case to investigate the regulation performance with the commonly used PI controller adjusting rules. The results showed that the regulation performance is not good and it is improved by a new proposed adjustment [10]. Martinez-Lucas et al., used pumped storage hydro plants with a PI governor tuning criterion to regulate frequency of a wind-solar isolated system and evaluated the regulation performance of the La Palma power system [11]. Parida and Chatterjee proposed an energy conversion mechanism of a wind−solar hybrid system with an improved control strategy, and the strategy was verified to regulate the system stability well [12]. With respect to reliability, Qin et al., proposed the Monte Carlo method to evaluate the generation and transmission system reliability of the coupled large-scale wind and photovoltaic power, then verified the effectiveness of the method in reliability evaluation based on the IEEE system with two PV power stations and wind farms [13]. Hu et al., evaluated the reliability of a power system with wind power and energy storage, and demonstrated the influence of wind energy dispatch restrictions (WEDR), wind farm location, etc., on the reliability benefits. One of the results showed that the maximum reliability benefits occurs at WEDR around 6% system load [14]. Zheng et al., established a complex uncertainty model considering the cost uncertainty of renewable energy power generation, then further proposed an improved model of integrated resource strategic planning. Finally, reliability is discussed by applying this model to a case study in China [15]. Billinton and Karki established the wind power model using the wind data and evaluated the generating capacity using the Monte Carlo simulation approach [16]. Hashemi-Dezaki et al., studied the influence of direct cyber-power interdependencies (DCPIs) on the reliability of smart grid reliability with micro turbine−wind-solar distributed generations, and the results showed that the DCPIs impacts gradually increase with the distributed generation penetration increases [17]. Li et al., investigated the impact of wind speed types, price volatilities and wind turbine number on the reliability benefits, and the results showed that wind speed types have a significant influence on benefits in a small−scale wind system [18]. Reliability is directly related to the economy and security operation of power system. There are lots of published papers on the reliability of the wind farm, wind−solar system, wind−hydro system or solar−hybrid system, etc. However, related studies on the regulation reliability of a PSPP in a wind−solar hybrid power generation system are scarce. Therefore, it is of significance to study the regulation reliability of the PSPP in a WSH hybrid system. Motivated by the above analyses, there are three advantages which make our study attractive compared with the prior work.
First, a mathematical model with a solar−wind−hydro hybrid power generation system is adopted to investigate the regulation reliability of PSPP. Second, the uncertainties and limitations of the model parameters are investigated to quantify the effect of regulation reliability. Third, five regulation indexes, i.e., rise time, settling time, peak value, peak time and overshoot of the reactive power, generator terminal voltage, guide vane opening and angular velocity, are extracted to evaluate the PSSP's regulation quality. Finally, the failure probabilities of power supply for the PSPP are calculated.
The rest of this study is structured as follows. The model of the wind−solar−hydro hybrid power generation system and the methods are presented in Section 2. Dynamic performances, uncertainty and sensitivity analysis are discussed in Section 3. The reliability analysis is studied in Section 4. Finally, conclusions are given in Section 5.

Model and Method
As shown in Figure 1, the hybrid solar-wind system with pumped storage system in this study is equipped with a Photovoltaic (PV) array, wind turbine, pumped storage power plant, an end-user (load) and a control station [19]. The mathematical model for each component of the hybrid system with pumped storage system is proposed in this section. Meanwhile, the uncertainty analysis and sensitivity analysis method are also described.  According to the IEEE group, the traveling wave transfer function between head and flow rate is written as [20] h where T 0 = L/α and Z 0 = αQ r /A i gH r . h q is the relative value of head change caused by flow change. T 0 is the elastic time of the equivalent penstock. α is the water hammer wave speed. L is the length of penstock. Z 0 is the surge impedance in per unit of the equivalent penstock. Q r and H r are the rated flow and head, respectively. A i is the section dimension of the penstock. g is the acceleration of gravity. q represents the relative value of flow. s is the Laplace operator.

Hydraulic Speed Regulation System
The servomotor is used to amplify the control signals and provide power to regulate the guide vane. The transfer function of the hydraulic servo system can be shown as [22] G(s) = 1 1 + T y s where T y is the engager relay time constant. A PID controller is commonly used. The transfer function of a PID can be expressed as [23] where k p , k i and k d denote the proportional, integral and differential adjustment coefficient, respectively. Substituting Equation (4) into Equation (3), one can get where ω and δ are the relative value of the generator rotor speed and the relative value of the rotor angle, respectively. y is the relative value of the guide vane opening.

Turbine
The nonlinearity of the turbine has a great influence on the stability calculation of the power system. Thus, the nonlinear model is considered in this section. Due to the fact that the efficient of the turbine is not 100%, the algebra equation of the IEEE Working Group Model can be used to calculate the output of the turbine [20]: where p m stands for the power output of the hydro turbine per unit. A t and q n1 denote the gain coefficient of the turbine and the no-loading flow per unit. D t and ∆ω represent the mechanical damping coefficient of the hydro turbine and the difference of angular velocity, respectively. The block diagram of the turbine is shown in Figure 2. The block diagram of the hydro turbine. A t is the hydro-turbine gain. h and q are the deviation of the water head and flow of the hydro turbine, respectively. fp is the head loss coefficient p m is the relative value of the output mechanical power. q n1 is the no-load flow deviation. ∆ω is the deviation value of the angular velocity of the unit. D t is the damping factor. hfc is the relative value of the pipe friction head loss. y is the guide vane opening. h q denotes the variation of the water head of the hydro turbine caused by the flow change of the penstock.

Excitation System
The stable operation of the power system is threatened due to the transient characteristics, such as sudden short-circuit faults [24]. The excitation system has the ability for additional damping to deal with these transients, which supplies a powerful guarantee for the safe and economical operation of the power system. Therefore, the excitation system is considered in this paper. A typical excitation system configuration is shown in Figure 3. The excitation system model of each unit is shown in Table 1. Figure 3. A typical excitation system configuration. U t is the generator terminal voltage. U ref is the reference voltage. U R is the output of the voltage regulator. E f is the excitation voltage. r f is the excitation winding resistance of generator. x ad is the inductance coefficient of the d-axis armature reaction. U f is the output of the excitation system stabilizer. U s is the output of the power system stabilizer. PSS stands for the power system stabilizer.

. Generator
A generator is the key component of the power system. First−order, three−order, five−order, and even seven-order generator models are commonly used to conduct research. However, it is a fact that the higher the model order, the higher the cost in terms of data requirements and calculation time [25]. Generally speaking, a third-order model has the advantage of simple structure and also considers the excitation system, which is widely used in dynamic analysis of the power system. Therefore, a third-order synchronous generator model is used in this paper, which is written as where δ, ω, and ω B are the rotor angle of the generator, the deviation of the relative angular speed and the nominal generator rotor speed, respectively. P m and P G represent the hydro turbine output power and the generator magnetic power, respectively. T j , D t and T d0 denote the inertia time constant of the generator, the generator damping coefficient and the generator time constant, respectively. X d∑ , X d∑ , U s and E f stand for the d-axis synchronous reactance, the d-axis transient reactance, the bus voltage and the controller output, respectively.

Pumped Storage Power Plant Model
According to the above analysis, the block diagram of the pumped storage power plant model is shown as Figure 4. In Figure 4, a step disturbance is set to investigate the dynamic response of output variables of PSPP.

Wind Power Generation System (WPGS)
The wind power generation system consists of the wind turbine, drive train, generator, pitch blade servo system and AC-DC-AV converter. The Double Fed Induction Generator (DFIG) is widely used in the wind farm due to its excellent operating performance; that is, the WPGS with DFIG attains lower requirements for power converter capacity, and flexible regulation of active and reactive power. [26]. The DFIG is connected directly to the power grid of 50 Hz via its stator, while the rotor is connected to the power grid through a power converter. A wind turbine is a complex non-linear mechanical device [27], and it ensures the conversion of the wind energy into mechanical energy by coupling generator with wind turbine [28].
The power versus wind speed curves are written as [17] where P WT and P rated are the power output of a wind turbine and the rated electrical power, respectively. v ci represents the cut−in wind speed. v r is the rated wind speed. v co stand for the cut−off wind speed. A, B and C are the intermediate variables, i.e., and respectively. Here, the wind farm cluster consists of 6 wind farms, each with a 1.5 MW wind turbine and a DFIGURE. It is connected to the external 25 KV AC grid through a 30 km and 25 KV transmission cable with the nominal power rating of 9 MW. The main purpose of this paper is to investigate the parameter uncertainty on the regulation performance of PSPP in a hybrid system. In order not to increase plant cost, the number of wind turbines and wind speed types remains unchanged. The WPGS model based on MATLAB/Simulink is shown in Figure 5.

Photovoltaic Power Generation System (PPGS)
Being environmentally friendly, solar energy has become one of the most suitable renewable energies [29], which is used in various ways and has the potential to be an alternative to conventional energy resources [30]. The most common application form of solar energy is photovoltaic power generation (PPG) [31]. A typical solar photovoltaic generator system consists of solar panel, DC/DC and DC/AC converter and the associated control. The key specifications of the solar panels are presented in Table 2.
where I ph is the photo current. I 0 is the diode saturation current. R s is the series resistance. R p is the shunt/parallel resistance. V t is the diode thermal voltage. The power output from a PV array is where N S and N p are the number of PV cells in a series for the studied array and the number of PV module in parallel, respectively. P A , I A and V A are the power output, current and voltage of a PV array, respectively. Due to the nonlinear photovoltaic cell array, the related controller is considered to maintain the photovoltaic power generation efficiency. Based on the above consideration, a generic model of PPGS developed in MATLAB/Simulink software is shown in Figure 6.

Model of the Wind−Solar−Hydro Hybrid System
According to the above analysis, the proposed WSH hybrid system schematic is demonstrated in Figure 7. Please see Ma et al., for more detailed information about the operating principle of the hybrid system [19]. This section supplies an interface to connect the wind, solar and hydro system model, which poses a potential reference to the multi−energy complementary system to better use clean energy.

Uncertainty Analysis
The mathematical model is crucial to simulate, study and predict the output of the WSH hybrid power generation system. The hybrid system model is complex, with nonlinear relations and many parameters. The important characteristic of these parameters is uncertainty, which has the ability to lead to the output uncertainty of the system [32]. Therefore, it is of importance to identify the relevance of these uncertainty parameters for the model.
One should note that uncertainty analysis often involves a larger number of parameters and data [33]. However, data of most parameters are generally very scarce and cannot be measured directly, which significantly affects the analysis results. Therefore, it is necessary to generate a large amount of random data within an acceptable range using a random sampling method. The Monte Carlo method has become a practical random sampling since it is simple, model-independent and generally applicable [32]. Please refer to Iooss and Le Gratiet (2019) for detailed information on the Monte Carlo method [34].

Sensitivity Analysis
The sensitivity analysis method is a valuable tool for building and using numerical simulation models [35]. It is used to study the influence of parameters uncertainty on the system output by setting the variation of model parameters in the corresponding design space [36]. Through sensitivity analysis, the model correction can set the parameters with low sensitivity to fixed values, and only calibrate the parameters that have great influence on the output variables, thus effectively simplifying the model, and improving the calibration accuracy of the model and saving time [37]. The first historical approach to sensitivity analysis is known as the local sensitivity where the effect of small input disturbance occurred around nominal values on the model output is studied [38]. To overcome the limitations of local methods (linearity and local variations), global sensitivity analysis considering the whole variation range of the inputs was developed in the late 1980s [35]. The extended Fourier amplitude sensitivity test (EFAST) is one of the most commonly used global sensitivity analysis methods, since it studies the influence of each parameter and the interaction of parameters on system output when multiple parameters change simultaneously [39]. Please see Xu et al., for detailed information about EFAST [40].

Reliability Analysis
Reliability analysis is usually performed when a structure is subject to uncertain influences [41]. The reliability, defined as the probability of the structure in a safe state, is where S and F represent the state domain and the failure domain, respectively. The surface separating S and F' is called the failure surface or limit state surface. x represents the possible value of the uncertain component. At present, the commonly used reliability calculation methods are the first-order reliability method (FORM), second-order reliability method (SORM), Monte Carlo method (MCS) and so on [42]. FORM is a popular one due to its simplicity, since only second moment information and the probability distribution type of the random variables are required to estimate the probability [43]. Based on the reliability indicators and checking points of the FORM, the SORM uses the second-order Taylor series expansion at the checking point to replace the original functional function, so as to improve the calculation accuracy of the FORM.

First-Order Reliability Method
The matrix formulation for a correlated normal of the Hasofer-Lind index (β) (also called first-order reliability index) can be expressed as [44] where X stands for the vector representing the set of random variables x i . µ, C and F are the vector of mean values, the covariance matrix and the failure domain, respectively. Low and Tang [45,46] presented an alternative interpretation of β based on the perspective of an expanding ellipsoid in the original space of the basic random variables, expressed as follows: where [R] stands for the correlation matrix. µ N i denotes the equivalent normal mean. σ N i is the equivalent normal standard deviation of random variable x i .
Based on the reliability index, the probability of failure can be evaluated as follows where ϕ(β) refers to the cumulative distribution function of the standard normal variable. P f is the probability of failure.

Second-Order Reliability Method
The SORM of the response surface g(u) = 0 is given by second-order Taylor series expansion at the design point u* in a standard normal U-space as [43] α is the directional vector at the design point in U-space. B is the scaled second-order derivatives of g(u) at u*, known as the scaled Hessian matrix.
The symbols used in this paper are shown in Table 3.

Dynamic Characteristics of WSH System in Steady and Fault States
In the hybrid power system model, two fault points, i.e., H and S, are set up at the output terminal of PSPP and infinite power supply to simulate the situation of the three−phase short circuit fault (TPSCF). The fault occurs at 1 s, is resected at 1.04 s, then the system gradually returns to normal operation. Meanwhile, the three−phase voltage and current of the output terminal of Wind Power Generation (WPG) (labeled as W) and grid-side (labeled as S) are measured, respectively. The corresponding numerical experiments are shown in Figure 8.   Figure 8a,b, it can be seen that the three−phase voltage and current are in a stable periodic motion when t < 1 s. With TPSCF occurring, the three−phase voltage of points W and S decreases, while the three-phase current of points W and S increases. The three−phase voltage and current of each point return to normal values and finally reach a stable state after the fault is removed at t = 1.04 s. Figure 8c,d illustrate the dynamic characteristics of three−phase current and voltage of points W and S with a three−phase short circuit fault of point S occurring at 1.0 s and being cleared at 1.04 s, respectively. From Figure 8c,d, the three-phase voltage and current of each point spread in a stable periodic motion when t < 1 s. The three-phase voltage of points W and S decreases significantly, while the three-phase current of points W and S increases between t = 1 s and t = 1.04 s. Specially, the three−phase voltage of points W and S returns to its normal value and finally reaches a stable state after the TPSCF is removed. The three−phase current of points W and S transforms from the periodic motion to an undulant state during the TPSCF. It worth noting that one of the three−phase currents of point S is separated from the others during the TPSCF. After the TPSCF is removed, the three−phase current of point W rises to its normal value and finally reaches a stable state.
However, one of the three−phase currents of point S is still separated from the others, while the others finally return to a stable state.
The above simulation results have proved that the established model is feasible which can be used to study the reliability analysis of the wind-solar-hydro hybrid power generation system in the following subsections.

Dynamic Performance Indexes (DPIs)
The response quality of the system under disturbance is usually measured by a dynamic performance index of the output variable. It is of practical significance to study the dynamic characteristics of the system by discussing the relationship between the DPIs and the system parameters. The commonly used DPIs are rise time (t r ), settling time (t s ), peak value (p), peak time (t p ) and overshoot (Os), which are used to characterize the response rapidity and stability of the system. Please refer to Appendix A for more details about DPIs. The corresponding results of the DPIs with Ke and Ki changing are shown in Tables 4 and 5. The bar in each cell indicates the relative magnitude of the values with the same color. Ke and Ki are the exciting gain and the integral adjustment coefficient, respectively. t r , t s , p, t p and Os are the rise time, settling time, peak value, peak time and overshoot, respectively. Table 4 shows the statistics of the DPIs of the reactive power and generator terminal voltage with Ke and Ki changing. As for the DPIs of the reactive power, it can be seen that the difference in rise time between different Ke and Ki settings is relatively small compared with that of settling time. The peak value and peak time remain almost unchanged with Ke or Ki changing. Settling time increases with the increases of Ke with a fixed Ki, while rise time shows the opposite trend. The maximum and minimum of settling time are 1.39728 and 0.70833, respectively. It worth noting that Ki has little effect on the DPIs of the reactive power when Ke remains unchanged. With regard to the DPIs of the generator terminal voltage, Ke has influences on the rise time, settling time, peak value and overshoot with different degrees. Specifically, with the increases of Ke, the rise time and peak value decrease, the settling time and overshoot values increase, while the peak time remains unchanged. When Ke = 8, the settling time and overshoot reach the maximum 1.9153 and 2.46662, meaning that a larger Ke value results in a poor rapidity and stability of the system response. The best overall quality of regulation occurs in Ke = 6 which has a relatively smaller settling time and overshoot. In addition, there is little change in the rise time, settling time, peak value, peak time and overshoot with Ki changing, which means that Ki has almost no effect on the DPIs. The bar in each cell indicates the relative magnitude of the values with the same color. Ke and Ki are the exciting gain and the integral adjustment coefficient, respectively. t r , t s , p, t p and Os are the rise time, settling time, peak value, peak time and overshoot, respectively. Table 5 displays the statistics of the DPIs of the guide vane opening and angular velocity with Ke and Ki changing. From Table 5, Ke and K i have little influence on the rise time, settling time, peak value and peak time of the guide vane opening since the difference in each DPI result between different Ke and Ki settings is relatively small. However, the overshoot values of the guide vane opening decrease with the increases of Ke and Ki. The maximum and minimum of the overshoot of the guide vane opening are 189.362 and 181.57, occurring in simulation No. 16 and No. 12, respectively. That is to say that a smaller setting of Ke and Ki causes a slower governor movement, leading to a larger overshoot of the guide vane opening. The above results show that the different Ke and Ki settings have a significant influence on the regulation quality of the guide vane opening. As for the DPIs of the angular velocity, it can be seen that both Ke and Ki have almost no effect on rise time, peak value and peak time. The maximum of settling time occurs in simulation No. 16 where Ke = 6 and Ki = 0.1. In addition, the values of overshoot at Ke = 6 are larger than those at Ke = 7 and Ke = 8. This means that the greater Ke value is, the better the dynamic performance of the system is.

Uncertainty Analysis
In this section, the Monte Carlo method is used to analyze the influence of WSH parameters on the output of PSPP. The iteration step is 1000, and the initial values of Ke and Ki are 7 and 0.25, respectively. Other parameters are set to conform to the normal distribution as shown in Appendix B. For the sake of brevity, only the parameters with definite influence rules on the output of PSPP are given. Graphics on diagonal lines indicate that the values of system parameters are in accordance with the normal distribution. The corresponding numerical results are shown in Figures 9-12. Figure 9 demonstrates the influence of WSH system parameters on the DPIs of the reactive power of PSPP. More specifically, Figure 9a is the influence of Ka on the rise time. Figure 9b is the influence of Ka on the overshoot. Figure 9c is the influence of Td0 on the settling time. Figure 9d is the influence of Tq00 and L1s on the p and pt. From Figure 9a, the rise time decreases with the increases of Ka. The overshoot shows an opposite trend compared with that of the rise time; that is to say that the overshoot increases with Ka increasing as shown in Figure 9b. Figure 9c shows that the settling time increases with Td0 increasing. From Figure 9d, with the increases of Tq00 and Lls, the peak value increases.
However, the value of the peak time is discontinuous, that is, it always equals 0.00666 or 0.00675 regardless of what the parameters are. The above results show that both parameters of PSPP and WPGS have a deterministic effect on the DPIs of reactive power, while PPGS parameters have no regular influence. In other words, WPGS has the ability to regularly affect the DPIs through the coupling effect of subsystems.     Figure 10 shows the influence of WSH system parameters on the DPIs of the generator terminal voltage of PSPP. Figure 10a is the influence of Td0 on the settling time. The changing trend of settling time is consistent with that of Td0, that is, the value of settling time increases with Td0 increasing. The changing trend of the peak value is decreasing as Ka increases in Figure 10b. Figure 10c shows the influence of Lm and H1 on the peak value and peak time. From Figure 10c, the peak value of the generator terminal voltage increases with Lm increasing. No matter how Lm and H1 change, the maximum peak time does not exceed 10. From Figure 10d, the value of overshoot rises due to the increases of Ncellm12 and Ir, indicating that the larger the capacity of photovoltaic power generation is, the greater the adverse effect on the voltage stability. From Figure 10a-d, parameters of WPGS, PSPP, and PPGS have regularity impacts on DPIs of the voltage. That is to say that WPGS and PPGS have a regular impact on the terminal voltage due to the coupling of subsystems. Figure 11 presents the influence of parameters of WSH on the DPIs of the guide vane opening of PSPP. Figure 11a Figure 11c shows that the value of overshoot decreases when the values of F increase. From the above analysis, it can be seen that the influence rule of different parameters is different due to parameter uncertainty. Meanwhile, only the uncertain parameters of PSPP show certainty influence on the DPIs of guide vane opening, which means that parameters of wind and solar subsystems have no regular influence on guide vane opening. Figure 12 shows the influence of WSH parameters on the DPIs of the angular velocity of PSPP. Specifically, Figure 12a From the analysis of Figures 9-12, different subsystems show different influence on the DPIs of PSPP output. It also can be obtained that the effects of PSPP parameters on each DPI have a certain regularity, while the influence degree of each PSPP parameter is different. This is due to the fact that system parameters are uncertain, and DPIs of this paper mainly depend on PSPP. The above results mean that to better maintain the regulation of PSPP, it is important to consider parameters' uncertainty and the coupling effect of subsystems.
To reveal the response speed of the WSH system, the cumulative probability distribution of rise time and settling time are studied in the following contents. The sampling times are 1000 times. The corresponding numerical results are shown in Figures 13 and 14. Figure 13a-d shows the cumulative probability distribution for the rise time of reactive power, generator terminal voltage, guide vane opening and angular velocity, respectively. From Figure 13a, it can be seen that the cumulative probability of rise time of reactive power changes as an "S" curve. No matter what the system parameter values are, the cumulative probabilities are 1.273% and 98.61% with the value of the rise time less than 1.101 × 10 −4 and 1.599 × 10 −4 . These phenomena mean that most of the rise time values are less than 1.599 × 10 −4 . From Figure 13b, the cumulative probability curve is relatively steep, indicating that the rise time value of the generator terminal voltage is comparatively centralized. The cumulative probability of rise time values less than 0.1406 is 99.34%, while the cumulative probability is almost equal to 0 with rise time value less than 0.08581. That is to say that the rise time value is in the range of 0.08581 and 0.1406. From Figure 13c, the cumulative probability is 99.62% when t r < 0.05674, and the cumulative probability is 0.4066% when t r < 0.01501. From Figure 13d, most of the values of rise time are less than 1.468 × 10 −4 , where the cumulative probability is 99.41%. Meanwhile, the slope of rise time cumulative probability of angular velocity curve changes smoothly compared with that of reactive power, generator terminal voltage and guide vane opening. The above phenomena show that there are great differences in the rise time of different output variables, especially reactive power and angular velocity. The cumulative probability curve in Figure 13b changes faster than those in other subgraphs.  Figure 14 shows the cumulative probability distribution for the settling time of the reactive power, generator terminal voltage, guide vane opening and angular velocity. It can be seen that all the cumulative probability curves are similar to "S". The slope of the cumulative probability curve of the angular velocity is larger than that of the reactive power, generator terminal voltage, guide vane opening; that is, the settling time distribution of angular velocity is relatively concentrated. Specifically, when t s < 1, the cumulative probability of angular velocity is larger than that of reactive power, generator terminal voltage and guide vane opening. When t s < 1, the cumulative probabilities of reactive power, generator terminal voltage, guide vane opening and angular velocity are 0, 0, 0.3797% and 66.58%, respectively. When t s < 2, the cumulative probabilities of reactive power, generator terminal voltage, guide vane opening and angular velocity are 0, 0, 99.83% and 100%, respectively. When t s < 4, the cumulative probabilities of reactive power, generator terminal voltage, guide vane opening and angular velocity are 54.91%, 82.48%, 100% and 100%, respectively. From the comparative results, the possible value of the settling time of the guide vane opening and angular velocity are larger than that of reactive power and generator terminal voltage in the case of large probability. The cumulative probability distribution of reactive power, generator terminal voltage, guide vane opening and angular velocity is significantly different from each other. From Figures 13 and 14, for the same DPI, the cumulative probability distributions of different output variables are significantly different from each other. This is because system parameters have different influences on system output as discussed in Section 3.3. Regarding different DPIs, the cumulative probability distributions of the same output variable are also different. In general, the settling time is larger than the rising time. This is due to the fact that the rise time is the time required for the response curve to reach the steady value for the first time, while the settling time is the time when the error between the unit step response and the steady value reaches the accepted value. In other words, the combination of Figures 13 and 14 can more clearly reflect the response speed of the system after disturbance.

Sensitivity Analysis
The sensitivity of model parameters often leads to the uncertainty of the model simulation. Therefore, finding out the sensitive parameters is significant to study the influence of effective parameters on the model output variables. Here, the Extended Fourier Amplitude Sensitivity Test (EFAST) method is used to study the sensitivity. The change law of system parameters is set to normal distribution as shown in Appendix B. The Monte Carlo method is selected for parameter random sampling, and the sampling times are 1000 times. The corresponding results are shown in Figures 15-19 and Table 6. Figure 15 shows the sensitivity index of parameters of the rise time of angular velocity. Specifically, Figure 15a Table 6. From Figure 15b and Table 6, T0 has the strongest impacts on the rise time of the angular velocity with sensitivity index 56.99%, followed by H 1 (9.298%), I r (2.859%), D t (2.719%), T q00 (2.144%), Ncellm12 (1.878%), b p (1.51%), K d (1.507%), T (1.478%) and q nl (1.388%). The total contribution rate of the top 10 sensitive parameters is 81.77%, meaning that these parameters have a direct effect on the rise time of angular velocity and the most significant factors affecting the rise time are identified through sensitivity analysis. The contribution rate of other parameters is less than 1.3%, indicating that the sensitivity of interaction among these parameters is small and the parameters are independent. In addition, it is worth noting that the second and the third sensitivity parameters are H 1 and I r coming from WPGS and PPGS, respectively. These phenomena mean that these parameters have the ability to indirectly influence the angular velocity of PSPP by interacting with other parameters.       Figure 16 demonstrates the sensitivity index of parameters of the settling time of angular velocity. From Figure 16a, it can be seen that different system parameters affect the settling time to a different degree. There are 10 sensitive parameters that have a relatively greater influence on settling time, as shown in Figure 16b and Table 6. Specifically, H has the greatest influence on settling time with sensitivity index 33.7%. The others are T q00 (22.29%), A t (4.851%), F (3.619%), K p (1.829%), H 1 (1.674%), R r (1.671%), F 1 (1.645%), K a (1.596%) and R s (1.57%). The total contribution rate of the top 10 sensitive parameters is 74.45%, meaning that the most significant factors affecting the output are studied and identified through sensitivity. Therefore, the influence of the top 10 sensitive parameters on the settling time should be fully considered in the numerical simulation of the WSH hybrid system. The sensitivity index of other parameters is less than 1.6%, indicating that the sensitivity of interaction among these parameters is small and the parameters are independent. Figure 17 shows the influence degree of parameters on the peak value of the angular velocity. Figure 17a shows the impact of 25 parameters on the peak value of the angular velocity. The top 10 sensitive parameters that have a relatively larger influence on the peak value are shown in Figure 17b. The corresponding results including the sensitivity index and the sensitivity ranking are shown in Table 6. From Figure 17b and Table 6, the most sensitive parameter is H with a contribution rate 72.55%, indicating that the H value directly determines the peak value of the angular velocity. Therefore, more attention should be paid to H to main the system stability. The second one is A t with a sensitivity index 5.44%, followed by F, T q00 , q nl , K a , R s , T d00 , H 1 and K p with sensitivity index 2.693%, 2.663%, 2.073%, 2.053%, 1.779%, 1.581%, 1.333% and 1.286%, respectively. The total contribution rate of the top 10 sensitive parameters to the peak value is 93.45%, which indicates that these parameters have a significant influence on the peak value of the angular velocity. Figure 18 displays the contribution rate of different parameters on the uncertainty of peak time of the angular velocity. Figure 18a is the main effect of the 25 parameters on the peak time of the angular velocity. Figure 18b shows the top 10 sensitive parameters that have a greater impact on the peak time. Rankings of the top 10 sensitive parameters are shown in Table 6. From Figure 18b and Table 6, the sensitivity ranking is H, T 0 , Ncellm12, T q00 , L ls , K a , T, R s , I r and A t . The corresponding sensitivity indexes are 71.59%, 2.355%, 2.231%, 1.988%, 1.844%, 1.681%, 1.444%, 1.422%, 1.318% and 1.28%, respectively. The total contribution rate of the top 10 sensitive parameters is 87.15%, that is, the top 10 sensitive parameters have a significant influence on the peak time of the angular velocity. In other words, the influence of different parameters on peak time varies greatly. The Ncellm12 is the third sensitivity parameter coming from PPGS, indicating that the parameter of PPGS has the ability to indirectly affect the angular velocity by interacting with other parameters. In addition, H has the greatest impact on the peak time consistent with that of peak value and settling time, indicating that the most sensitive parameters of these DPIs are consistent. Figure 19a,b describe the main effects of 25 parameters and the top 10 sensitive parameters on the overshoot of the angular velocity, respectively. The top 10 sensitive parameters are T 0 , q nl , R s , L m , Ncellm12, H 1 , K d , f p , I r and H, respectively. The sensitivity results are shown in Table 6. From Table 6, the maximum sensitivity index is 3.177% coming from T 0 , and the minimum sensitivity index is 1.328% coming from H. It also can be seen that the total contribution rate of the top 10 sensitive parameters is 17.764%. It is worth noting that the sensitivity index value is relatively small compared with that of rise time, settling time, peak value and peak time. This phenomenon means that although many factors affect the overshoot of angular velocity, the difference of influence degree is small.
From Figures 15-19 and Table 6, it is concluded that the sensitivity degree of different DPIs to system parameters is obviously different. This phenomenon means that even the same parameter has a different effect on the response speed and response stability of the system. Moreover, parameters of WPGS and PPGS have a significant influence on DPIs, indicating that these parameters have the ability to indirectly affect the angular velocity of PSPP by interacting with other parameters.

Reliability Analysis
The regulation reliability of WSH hybrid power generation system is directly related to the balance between the power supply and demand. Therefore, it is of great significance to study regulation reliability to maintain the safe and economic operation of the power system. In this section, we study the influence of the system parameters on the output of the WSH system. Here, the peak value of angular velocity is selected as an example, and the number of simulations is 2000. Figure 20 is the distribution of the peak value of angular velocity. From Figure 20, the abscissa represents the peak value of angular velocity, which is distributed between 0.017 and 0.034. The ordinate stands for the number of corresponding peaks. It is clear that the peak value of angular velocity approaches the normal distribution well. Most of the peak values are in the range of 0.022 and 0.024, and the values on both sides are relatively small. Here, the peak value greater than 0.028 is defined as failure range. According to this definition, some values in Figure 20 are in the failure range, and the probability of this part is called the failure probability. To avoid the peak value falling in the failure range, the system parameters should be adjusted based on the results of uncertainty analysis and sensitivity analysis to meet the steady state operation. In addition, the cumulative probability diagram is plotted in Figure 21 to obtain the failure probability.   Figure 21 demonstrates the cumulative probability of the peak value of the angular velocity. The blue solid line is the cumulative probability curve for the peak value, and the red dashed line is the 95% confidence bounds. From Figure 21, it can be seen that the blue solid line is within the 95% confidence bounds, meaning that the simulation results are proved to be reliable. In addition, the cumulative probability of the peak value less than 0.028 is 97.5%, that is, the reliability probability of the power supply is 97.5%. In other words, the probability that the primary objective ensuring the balance between supply and demand cannot be satisfied is 2.5%, meaning that consumers cannot receive the electricity they need. This case may overburden the power system, cause widespread power blackout and finally pose a threat to the safe and economic operation of the power system.

Conclusions
To investigate the regulation reliability of PSPP in a multi-energy power system, a wind-solar-hydro hybrid power system model is established. Based on the established model, uncertainty and sensitivity analysis of system parameters were carried out using the Monte Carlo method and EFAST. Finally, the regulation reliability of the pumped storage power plant in a WSH hybrid power generation system was also discussed. The main results are as follows: (1) The influence rules of the model parameters on the WSH hybrid system are obtained from the uncertainty analysis. Parameters of the wind, solar and hydro subsystem show the different influence on DPIs of the PSPP output due to parameters uncertainty. Both PSPP and WPGS parameters have a deterministic effect on the DPIs of reactive power, while the influence of PPGS has no regularity. The uncertain parameters of WPGS, PSPP and PPGS have regularity influence on the DPIs of the generator terminal voltage. Only PSPP parameters show certainty influence on the DPIs of the guide vane opening and angular velocity. The results also mean that the coupling effect of subsystems has the ability to affect the DPIs of PSPP in a certain case. (2) For the same DPI, the cumulative probability distributions of different output variables are significantly different from each other. Regarding different DPIs, the cumulative probability distributions of the same output variable are also different. In general, the settling time is larger than rising time. (3) The sensitivity degree of different DPIs to system parameters is obviously different, and even the same parameter has a different effect on the response speed and response stability of the angular velocity. The total contribution rate of the top 10 sensitive parameters on the rise time, settling time, peak value, peak time and overshoot of the angular velocity is 81.77%, 74.45%, 72.55%, 87.15% and 17.764%, respectively. Meanwhile, parameters of WPGS and PPGS have the ability to indirectly affect the angular velocity of PSPP by interacting with other parameters. (4) The peak value of angular velocity is distributed between 0.017 and 0.034. Most of the peak value of the angular velocity is in the range of 0.022 to 0.024, and the values on both sides are relatively small. There is a 2.5% probability that the system cannot meet the requirements of operation reliability, which may have a bad impact on the corresponding equipment or even threaten the normal operation of the system. This paper takes parameters uncertainty into account to investigate the regulation characteristics of PSPP in a WSH hybrid power system. Only reactive power, generator terminal voltage, guide vane opening and angular velocity are considered as the research objects. In the future work, similar studies can be conducted to investigate the influence of parameters on other output variables associated with the wind, solar or hydro subsystem.