Low Frequency Oscillations in a Hydroelectric Generating System to the Variability of Wind and Solar Power

Abstract

The penetration of multiple integrated renewable energies to the power grid are relevant for decision making in energy policy, environment and business. Such an electricity penetration is affected by the intermittent and volatile characteristics of integrated energies, mostly significantly related to the safe and stable electricity production and supply in real world. Here, this paper focuses on the low frequency oscillation analysis of the hydropower generation response to the wind and solar variability. To enable this analysis, a hybrid model of hydropower system integrating with the wind and solar power system is presented. The Nyquist and root-locus stability methods are used to investigate the sensitivity performance of the hydropower governor to the fluctuation of the integrated renewable energies. Additionally, to quantify the risk of the hybrid system, the low frequency oscillation response of hydropower system to wind/solar/hydropower quota and transmission line distance ratio is extensively investigated in this study. The results show that under the case of the wind, solar and hydropower ratio is 40:1:150, the optimal values for maximally reducing hydropower low frequency oscillation are finally determined as kp = 0.8, ki = 0.25 and kd = 0.5. Regarding a certain wind/solar/hydropower quota, it is a promising strategy to increase the solar-load transmission line in order to achieve the safe and stable operation of the hybrid system and a relatively excellent dynamic regulation capacity of the hydropower governor. The model, methods and results implemented in this study are exploited to markedly improve new knowledge applications, policy management, low carbon emissions and investment competitiveness of future energy systems.

1. Introduction

Renewable energy-based electricity generation is on the rise to supplement and replace conventional thermal-based power production. The worldwide newly added renewable electricity capacity accounted for about 50% started in 2015—of which wind and solar power shared almost 70% [1]. This stimulates the integration of multiple renewable energy generation, for example, Denmark, Texas, China, and the United States and California are the typical leaders to push the development of this matter [2]. In conclusion, the penetration of renewable energy integration to the power grid, especially for the wind and solar energies, until recently thought to be a long-distant future scenario, is taking place right now [3,4,5].

However, the volatility and fluctuation of wind speed and solar irradiation in different hours and days present some challenges to the electricity operation and management, which raises a new problem on how to balance the production and consumption in the electricity industry [6,7,8]. In other words, the unpredictable wind and solar power injection increases the uncertainty and unbalance of the power grid; meanwhile, this uncertainty significantly leads to the increase of the reserve and imbalance operation costs and the fast growth of global power price fluctuation [9,10,11,12]. One of the most promising solutions to this challenge is the application of energy storage systems (ESSs) [13,14]. The ESS plays a damper and regulator role in grid-connection of renewable energies, and its storage status is generally summarized in three aspects [15]: (i) decoupling the timing of generation and consumption, (ii) ensuring the peak shaving and switching balance from one power energy to another one, and (iii) regulating electricity end-use quality. To underline these goals, so far the hydropower system, which stores the surplus water energy in peak hours and consumes this energy at peak-load hours, is recognized as an ideal ESS for enhancing the economic and environmental profits, the transmission congestion mitigation, and the reduction of wind/solar curtailment in renewable energy integration [16,17].

The participation of a hydropower system is a proper supplement to tackle the uncertainty problems of wind and solar generation. The explanation for this is the capabilities of a hydropower system in the excess energy storage in noticeable amounts and the fast response to the use of operating reserves through the regulation between upstream and downstream reservoirs [18,19]. Hydropower plays an important role in improving energy utilization efficiency and reducing output fluctuation in a hydro–wind–photovoltaic system. In [20], the proposed three indices have the ability to show the stability and complementarity characteristics of the hybrid system with the time scale varying from second to hour. Wang et al. [21] refined the model of a hydro–wind–photovoltaic system and investigates the complementary mechanism underlying multi-energy systems. Hence, the issue of understanding on how wind and solar power injection affects the operational stability of hydropower system, as one of the practical challenges of the electricity network, is the key focus of this work [22,23,24]. Other merits of the hydropower system in the hybrid renewable energy generation are concluded as follows:

• Fast responses to frequency and voltage variation under the appropriate controller. Energy planers and regulation entities cheer [25,26].
• Balancing power price volatility and reducing capital and operating expenses. Storage companies and energy investors cheer [27,28].
• Hydropower greatly promotes the large-scale applications of renewable energies, considerably reducing greenhouse gas (GHG) emissions and improving public health benefits. Governments and environmental organizations cheer [29,30].

Considering the significant contributions of hydropower system to the electrical generation of renewable energy integration, some of the research papers in this regard are from an optimized scheduling viewpoint which involves in both wind and solar joint energies. In [31,32,33], the optimized algorithms aimed at providing the certain operational observations, hydropower generation schedule and scheduling cost assessment for the wind/solar/hydropower hybrid system. However, while accounting for hydropower electricity, these studies do not focus on the impact of the frequent optimal dispatch on the fatigue property of hydropower components. In [34,35,36,37], the synergetic scheduling model and control strategy of wind/solar/hydropower units are the main focuses, but not involving in the detailed hydraulic–electric–mechanical internal characteristics. Moreover, some papers are from the economic, policy and environment viewpoint. In [38,39,40,41,42], the market value and environmental benefits of the hybrid system is successfully quantified from the perspective of data statistics and the corresponding policy support is also recommended. Nevertheless, a quantitative analysis of the coordinating influences between the hydropower system and wind/solar systems is lacking. This paper aims to fill these gaps.

Based on the above considerations, this paper is devoted to the stability interaction between the hydropower system and the intermittent wind/solar power. The purpose of this study is to analyze the internal stability mechanism in the hybrid system and to quantify the influence of wind/solar power uncertainty on the hydropower low frequency oscillation. Innovation aspects of this paper are:

• A promising hybrid system of the hydropower generation integrating with wind farm and solar photovoltaic system is established using MATLAB/Simulink, in order to enable the stability analysis. This contributes to the current international pool of the integration modelling knowledge.
• The sensitivity of hydropower low frequency oscillations to its governor regulation capacity is quantified under the volatility influence of wind and solar energies. The main adopted methods include Nyquist response and root-locus analysis.
• To understand the stability conditions of the hybrid system, the influence of different wind/solar/hydropower quotas (i.e., W: S: H) and the various transmission line distance ratios on the low frequency oscillation mode of hydropower system are also quantified. The assessment indicators in this part include the three-phase parallel RLC load, the grounding transformer, the three-phase PI section line and the wind-farm transmission line, and the assessment criteria are oscillation frequency and damping ratio.

The remainder of this paper is organized as follows. Section 2 presents a detailed model of the integrated wind/solar/hydropower power system. The low frequency oscillation sensitivity of hydropower system to PID governor regulation capacity is then investigated in Section 3. In Section 4, results for the low frequency oscillation response to renewable-quota and transmission line distance ratio are provided and discussed. Finally, conclusions are drawn in Section 5.

2. Mathematical Model of the Hybrid Power System

In this section, a hybrid system comprising the hydropower system, wind farm and solar photovoltaic system is established using MATLAB/Simulink [43] (see the schematic diagram in Figure 1). The hydropower system is composed of the penstock, governor, hydro-turbine, synchronous generator and excitation sector. The wind farm mainly employs a doubly fed induction generator (DFIG) and wind turbine. Additionally, the solar photovoltaic system focuses on the blocks of the photovoltaic array, photovoltaic inverter and maximum power point tracing (MPPT) controller. The detailed modeling process is illustrated as follows.

2.1. Hydropower System

2.1.1. Penstock

The penstock includes two sectors, i.e., the water diversion pipe and steel pressure pipe, and its dynamic behavior is expressed by elastic water hammer equations, as follows:

$$\left\{\begin{array}{l}\frac{\partial Q}{\partial t}-gA\frac{\partial H}{\partial L}+\frac{f{Q}^2}{2AD}=0\\ a^2\frac{\partial Q}{\partial L}-gA\frac{\partial H}{\partial t}=0\end{array}\right.$$

(1)

where Q, H, L, f, D, A and a denote the flow of a pipe section, head of a pipe section, distance between the studied pipe section and calculated origin, water flow friction coefficient, pipe diameter, pipe area and water-hammer velocity.

If $q=∆Q/Q_r$, $h=∆H/H_r$ and $l=∆L/L_r$, Equation (1) can be rewritten as:

$$\left\{\begin{array}{l}{T}_W\frac{\partial q}{\partial t}-\frac{\partial h}{\partial l}+f_1(q_0^2+2q_{0}q)=0\\ \frac{4T_W}{T_r^2}\frac{\partial q}{\partial l}-\frac{\partial h}{\partial t}=0\end{array}\right.$$

(2)

where T_W is the water flow inertia time constant, $T_W=\frac{L_r{Q}_r}{{gAH}_r}$. T_r is the water hammer pressure time constant, $T_r=\frac{2L_r}{a}$. The variable f₁ is represented by $f_1=\frac{f{Q}_r^2{L}_r}{2D{A}^{2}g{H}_r}$, wherein parameters L_r, Q_r and H_r are the total pipe length, design flow and design head.

Equation (2) experiences the Laplace transform and then is expressed as:

$$\left\{\begin{array}{l}h(1,S)=ch(\frac{T_r}{2}s+f)h(0,S)+2h_{W}sh(\frac{T_r}{2}S+f)q(0,S)\\ q(1,S)=\frac{1}{2h_W}sh(\frac{T_r}{2}s+f)h(0,S)+ch(\frac{T_r}{2}S+f)q(0,S)\end{array}\right.$$

(3)

where the hyperbolic function is $sh\left(x\right)=\frac{e^x-e^{-x}}{2},ch\left(x\right)=\frac{e^x+e^{-x}}{2}$.

If the upstream pipe head $h(1,S)$ and the downstream pipe flow $q(0,S)$ are measured, thus the upstream pipe flow and the downstream pipe head for the water diversion pipe can be obtained as:

$$\left\{\begin{array}{l}h(0,S)=\frac{1}{ch(\frac{T_r}{2}S+f)}h(1,S)-2h_{W}th(\frac{T_r}{2}S+f)q(0,S)\\ q(1,S)=\frac{1}{ch(\frac{T_r}{2}S+f)}q(0,S)+\frac{1}{2h_W}th(\frac{T_r}{2}S+f)h(1,S)\end{array}\right.$$

(4)

Similarly, the Laplace transform expression of the steel pressure pipe is obtained as:

$${{\displaystyle h}}_q={{\displaystyle Z}}_{01}q\tanh \left({{\displaystyle T}}_{01}s\right)$$

(5)

where h_q, Z₀₁ and T₀₁ are transient head of hydraulic turbine, standardized value of hydraulic surge impedance of pipeline and elastic time, respectively.

Ignoring the higher order term $\tanh \left({{\displaystyle T}}_{01}s\right)$ in Equation (5), thus the transfer function from the hydro-turbine flow to hydro-turbine head can be rewritten as:

$${{\displaystyle h}}_q\left(s\right)={{\displaystyle Z}}_{01}\frac{{{\displaystyle \pi }}^2{{\displaystyle T}}_{01}s+{{\displaystyle T}}_{01}^3{{\displaystyle s}}^3}{{{\displaystyle \pi }}^2+4{{\displaystyle T}}_{01}^2{{\displaystyle s}}^2}q\left(s\right)$$

(6)

2.1.2. Governor

In this work, the governor adopts the proportional–integral–differential (PID) controller, where the proportional loop aims to shorten response time, the integral loop decreases the steady error, and the differential loop improves the overshoot response. Generally, the output signal of a universal PID controller is expressed as:

$$Y_{PID}=Y_P+Y_I+Y_D$$

(7)

and

$$\left\{\begin{array}{l}{Y}_P=K_P\Delta F\\ Y_I=\frac{K_I}{S}\Delta I=\frac{K_I}{S}\left[\Delta F+b_p(Y_c-Y_{PID})\right]\\ Y_D=\frac{K_{D}S}{1+T_{1v}S}\Delta F\end{array}\right.$$

(8)

where K_P, K_I and K_D are the proportional adjustment coefficient, integral adjustment coefficient and differential adjustment coefficient, respectively. △F, Y_c, b_p and T_1v are the relative deviation of the input frequency, relative predefined guide vane opening, permanent interpolation coefficient and differential loop time constant, respectively. If the permanent interpolation coefficient (b_p) equals to zero, thus Equation (8) is rewritten as:

$$G(S)=\frac{Y_{PID}(s)}{\Delta F(s)}=\left(K_P+K_I\frac{1}{S}+K_{D}S\right)$$

(9)

Based on the above considerations, the block of the PID controller is shown in Figure 2.

One of the important sector in hydro-electric subsystem is the hydraulic cylinder. The hydraulic cylinder sector that transforms the output signal of PID controller to the hydraulic signal can be modeled as Figure 3.

2.1.3. Hydro-Turbine

The hydro-turbine model uses the recommendation expression from IEEE Group, which is expressed as:

$${{\displaystyle p}}_m={{\displaystyle A}}_{t}h\left(q-{{\displaystyle q}}_{nl}\right)-{{\displaystyle D}}_{t}y\Delta \omega$$

(10)

where A_t, q, h, f_p, p_m, q_nl, △ω, D_t, and y denote the hydro-turbine gain, relative deviation of hydro-turbine flow, relative deviation of hydro-turbine head, relative head loss coefficient in penstock, mechanical output power, relative no-load flow, relative deviation of hydro-turbine rotational speed, damping factor, and relative deviation of guide vane opening, respectively. Thus, the block diagram of the hydro-turbine is performed in Figure 4.

2.1.4. Synchronous Generator

The three-phase synchronous generator is achieved in the DQ rotor framework, and its stator winding uses Y type connection form to link with the internal neutral point. The three-phase synchronous generator works at the generator/electric motor mode, and the electrical part of the machine is represented by a sixth-order state-space model and the mechanical part is the same as in the simplified synchronous machine block [44]. The model takes into account the dynamics of the stator, field, and damper windings [45]. The voltage equation and flux equation of the six-order synchronous generator are:

$$\left\{\begin{array}{l}{{\displaystyle V}}_d={{\displaystyle R}}_s{{\displaystyle i}}_d+\frac{d}{dt}{{\displaystyle \phi }}_d-{{\displaystyle \omega }}_R{{\displaystyle \phi }}_q\hfill \\ {{\displaystyle V}}_q={{\displaystyle R}}_s{{\displaystyle i}}_q+\frac{d}{dt}{{\displaystyle \phi }}_q-{{\displaystyle \omega }}_R{{\displaystyle \phi }}_d\hfill \\ {{\displaystyle V}}_{fd}^{\prime }={{\displaystyle R}}_{fd}^{\prime }{{\displaystyle i}}_{fd}^{\prime }+\frac{d}{dt}{{\displaystyle \phi }}_{fd}^{\prime }\hfill \\ {{\displaystyle V}}_{kd}^{\prime }={{\displaystyle R}}_{kd}^{\prime }{{\displaystyle i}}_{kd}^{\prime }+\frac{d}{dt}{{\displaystyle \phi }}_{kd}^{\prime }\hfill \\ {{\displaystyle V}}_{kq1}^{\prime }={{\displaystyle R}}_{kq1}^{\prime }{{\displaystyle i}}_{kq1}^{\prime }+\frac{d}{dt}{{\displaystyle \phi }}_{kq1}^{\prime }\hfill \\ {{\displaystyle V}}_{kq2}^{\prime }={{\displaystyle R}}_{kq2}^{\prime }{{\displaystyle i}}_{kq2}^{\prime }+\frac{d}{dt}{{\displaystyle \phi }}_{kq2}^{\prime }\hfill \end{array}\right.$$

(11)

and

$$\left\{\begin{array}{l}{{\displaystyle \phi }}_d={{\displaystyle L}}_d{{\displaystyle i}}_d+{{\displaystyle L}}_{md}\left({{\displaystyle i}}_{fd}^{\prime }+{{\displaystyle i}}_{kd}^{\prime }\right)\hfill \\ {{\displaystyle \phi }}_q={{\displaystyle L}}_q{{\displaystyle i}}_q+{{\displaystyle L}}_{mq}{{\displaystyle i}}_{kq}^{\prime }\hfill \\ {{\displaystyle \phi }}_{fd}^{\prime }={{\displaystyle L}}_{fd}^{\prime }{{\displaystyle i}}_{fd}^{\prime }+{{\displaystyle L}}_{md}\left({{\displaystyle i}}_d+{{\displaystyle i}}_{kd}^{\prime }\right)\hfill \\ {{\displaystyle \phi }}_{kd}^{\prime }={{\displaystyle L}}_{kd}^{\prime }{{\displaystyle i}}_{kd}^{\prime }+{{\displaystyle L}}_{md}\left({{\displaystyle i}}_d+{{\displaystyle i}}_{fd}^{\prime }\right)\hfill \\ {{\displaystyle \phi }}_{kq1}^{\prime }={{\displaystyle L}}_{kq1}^{\prime }{{\displaystyle i}}_{kq1}^{\prime }+{{\displaystyle L}}_{mq}{{\displaystyle i}}_q\hfill \\ {{\displaystyle \phi }}_{kq2}^{\prime }={{\displaystyle L}}_{kq2}^{\prime }{{\displaystyle i}}_{kq2}^{\prime }+{{\displaystyle L}}_{mq}{{\displaystyle i}}_q\hfill \end{array}\right.$$

(12)

where V, I, R and $\phi$ are the voltage, current, resistance and flux linkage of winding, respectively; ${\omega }_R$ is rotor rotation angular velocity; L_d, L_q, L^’_fd, L^’_kd, L^’_kq1 and L^’_kq2 are self-inductance of winding; and L_md and L_mq are mutual inductance of winding.

2.1.5. Excitation Sector

The excitation system provides the excitation power to regulate the generator voltage and control the reactive power distribution of the generator in parallel operation states. It has a great influence on the dynamic behaviors of the generator and has the ability to improve the stability of the power system. The model of the excitation system is shown in Equation (13). The block diagram of the excitation system is presented in Figure 5.

$$\mathrm{The} \mathrm{excitation} \mathrm{system} \mathrm{model}:\left\{\begin{array}{l}{\dot{U}}_R=\frac{1}{T_A}\left[K_A\left(U_{ref}+U_{PSS}-U_t-U_F\right)-U_R\right]\\ {\dot{U}}_F=\frac{1}{T_F}\left(K_F{\dot{E}}_f-U_F\right)\\ {\dot{E}}_f=\frac{1}{T_L}\left[U_R-\left(S_E+K_L\right)E_f\right]\end{array}\right.$$

(13)

Thus, in light of the above subsystem block modeling, the global hydropower system (run-of-river hydropower plant) is obtained in Figure 6.

2.2. Wind Farm

2.2.1. Wind Turbine Model

The blades in the wind turbine are used to capture the wind, and the output power of the wind turbine is closely related to the wind energy efficiency (C_p) and the swept area (A_r) of the blade. The operational phenomenon of wind turbines can be described by the aerodynamics theory, and based on this theory, the power captured by blades of the wind turbine can be expressed as:

$$P_W=A_r{\rho }_{air}{C}_p{v}_{\omega }^3$$

(14)

Where P_W, ρ, C_p, A_r, and v_ω are the mechanical power output of the wind turbine, air density, wind energy efficiency, swept area of the blade, and wind velocity, respectively. Here, the wind energy efficiency C_p can expressed as the ratio between the power captured by the wind turbine (P_w) and the total power through the wind turbine (P_a), i.e.,

$$C_p=\frac{P_w}{P_a}$$

(15)

As discussed in [36], Equation (15) can be expressed in the form of Equation (16), which is easier to perform the calculation and is written as follows:

$$\left\{\begin{array}{l}{C}_p=0.22(\frac{116}{{\lambda }_a}-0.4\beta -5)e^{-12.5/{\lambda }_a}\\ {\lambda }_a={(\frac{1}{\lambda +0.08\beta }-\frac{0.035}{{\beta }^3+1})}^{-1}\end{array}\right.$$

(16)

where λ_a is the tip velocity ratio and β is the pitch angle of the blade. Using all these equations, the wind characteristic curves can be obtained and described as shown in Figure 7.

2.2.2. Mechanical Drive Shaft Model

The mechanical drive shaft is an important component of the wind farm (WF), which coupled the wind turbine with the DFIGURE. In this work, a two-mass drive shaft model, as presented in Figure 8, is used, as it describes the energy transfer between medium/high and low speed shafts while comparing with the one-mass model. The dynamic characteristics of a two-mass mechanical drive shaft model can be expressed in terms of the following equations [37]:

$$\left[\begin{array}{l}{\dot{\omega }}_t\\ {\dot{\omega }}_g\\ {\dot{T}}_{ls}\end{array}\right]=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right]\left[\begin{array}{c}{\omega }_t\\ {\omega }_g\\ T_{ls}\end{array}\right]+\left[\begin{array}{c}{b}_{11}\\ b_{21}\\ b_{31}\end{array}\right]T_a+\left[\begin{array}{c}{b}_{12}\\ b_{22}\\ b_{32}\end{array}\right]T_m$$

(17)

with

$$\left\{\begin{array}{l}{b}_{11}=\frac{1}{J_r},b_{12}=0\\ b_{21}=0,b_{22}=-\frac{1}{J_g}\\ b_{31}=\frac{B_{ls}}{J_r},b_{32}=\frac{B_{ls}}{n_g{J}_g}\end{array}\right.$$

(18)

and

$$\left\{\begin{array}{l}{a}_{11}=-\frac{B_r}{J_r},a_{12}=0,a_{13}=-\frac{1}{J_r}\\ a_{21}=0,a_{22}=-\frac{B_g}{J_g},a_{23}=\frac{1}{n_g{J}_g}\\ a_{31}=(K_{ls}-\frac{B_{ls}{B}_r}{J_r}),a_{32}=\frac{1}{n_g}(\frac{B_{ls}{B}_r}{J_g}-K_{ls}),a_{33}=-B_{ls}(\frac{J_r+n_g^2{J}_g}{n_g^2{J}_g{J}_r})\end{array}\right.$$

(19)

where T_ls, T_t, T_m, ω_t, ω_g, J_r, J_g, B_r, B_g, B_ls and n_g are the low speed shaft torque, aerodynamic torque, generator electromagnetic torque, rotor speed, generator speed, rotor inertia, generator inertia, rotor external damping, generator external damping, low speed shaft damping, and ratio constant, respectively.

The combined block diagram representation of the wind turbine and mechanical drive shaft model is shown in Figure 9. As shown in Figure 9, the input of the wind turbine are the blade angle and wind speed, and the output is the aerodynamic torque. The input of the drive train are the aerodynamic torque and generator speed, and the output are the wind turbine speed and torque transmitted through the shaft. Then, the torque is transformed to the torque at the input of the generator.

2.2.3. DFIG Model in dq Frame

The DFIG exhibits a flexible working mode due to its AC voltage excitation system which easily achieves the variable-speed constant frequency (VSCF) generation. The Park transformation can be used to represent the dynamical model of the DFIG in a synchronously rotating dq frame and in this frame, the DFIG model can be expressed as [38]:

$$\left\{\begin{array}{l}{v}_{sd}=R_s\cdot i_{sd}+\frac{d{\Phi }_{sd}}{dt}-{\omega }_s\cdot {\Phi }_{sq}\\ v_{sq}=R_s\cdot i_{sq}+\frac{d{\Phi }_{sq}}{dt}-{\omega }_s\cdot {\Phi }_{sd}\\ v_{rd}=R_r\cdot i_{rd}+\frac{d{\Phi }_{rd}}{dt}-{\omega }_r\cdot {\Phi }_{rq}\\ v_{rq}=R_r\cdot i_{rq}+\frac{d{\Phi }_{rq}}{dt}-{\omega }_r\cdot {\Phi }_{rd}\end{array}\right.$$

(20)

where v_sd,q, v_rd,q, i_sd,q, i_rd,q, Φ_sd,q, Φ_rd,q, R_s, R_r, ω_s, and ω_r are the stator voltage in dq frame, rotor voltage in dq frame, stator current in dq frame, rotor current in dq frame, stator flux in dq frame, rotor flux in dq frame, stator resistance, rotor resistance, stator speed, and rotor speed, respectively.

The d- and q-axis are magnetically decoupled, and thus the flux-linkage equations can be written as:

$$\left\{\begin{array}{l}{\Phi }_{sd}=L_s\cdot i_{sd}+M\cdot i_{rd}\\ {\Phi }_{sq}=L_s\cdot i_{sq}+M\cdot i_{rq}\\ {\Phi }_{rd}=L_r\cdot i_{rd}+M\cdot i_{sd}\\ {\Phi }_{rq}=L_r\cdot i_{rq}+M\cdot i_{sq}\end{array}\right.$$

(21)

where L_s, L_r, and M are the stator inductance, rotor inductance, and mutual inductance, respectively.

Additionally, the electromagnetic torque in dq frame can be expressed as:

$$T_g=\frac{3}{2}{n}_p[{\Phi }_{sd}{i}_{sq}-{\Phi }_{sq}{i}_{sd}]$$

(22)

where n_p is the pole pairs of the DFIG.

2.2.4. PWM Converter Model

The PWM converter can be described through a state-average model, which represents the product of the voltage and duty ratio. Based on the principle of the state-space average model, the three-phase PWM converter in a DFIG-based WF is expressed as:

$$\left\{\begin{array}{l}{L}_f\frac{d{i}_a}{dt}=-R_f{i}_a+e_a-(d_a{v}_{dc}+v_{N0})\\ L_f\frac{d{i}_b}{dt}=-R_f{i}_b+e_b-(d_b{v}_{dc}+v_{N0})\\ L_f\frac{d{i}_c}{dt}=-R_f{i}_c+e_c-(d_c{v}_{dc}+v_{N0})\end{array}\right.$$

(23)

where d_aT and d^’_aT are the conduction times of switch S1 and switch S2, respectively; with T and d_a as the switching cycle of the PWM switches and the duty ratio of phase d_a; i_abc represents the currents flowing through three phases; e_abc represents the voltages at the grid connection point for three phases; R_f and L_f are the filter resistance and inductance, respectively; v_dc is the DC-link voltage at the input of the PWM converter; and v_N0 is the neutral voltage at the grid connection point. The subscripts a, b, and c with different variables represent the phase a, phase b, and phase c, respectively.

By applying Kirchhoff’s current law at the input side of the PWM converter where the DC-link capacitor is connected, it can be written as:

$$i_{dc}=d_a{i}_a+d_b{i}_b+d_c{i}_c$$

(24)

where i_dc is the Dc current at the input of the PWM converter.

The block diagram representation of the state-average PWM converter model is shown in Figure 10, whose outputs are actually average voltage control signals used for controlling the PWM switches in the converter.

2.3. Solar Photovoltaic System

The solar photovoltaic system involves in three components, i.e., the photovoltaic array, photovoltaic inverter and maximum power point tracing (MPPT) controller. Herein, the photovoltaic array contains several parallel modules, and the current-voltage property of each module is expressed by five operating parameters, i.e., current source (IL, which can utilize solar energy to generate current), diode parameters (I0 and NI), series resistance (rs) and shunt resistance (rsh). These parameters are closely to the variation of temperature and irradiation. Besides, the photovoltaic inverter is an important element to ensure system stability during the grid-connected transient process, which can convert the DC to sinusoidal AC in order to successfully link with power grid. Moreover, maximum power point tracing (MPPT) controller is applied in this study to maximize the solar electricity production. The block diagram if the proposed solar photovoltaic system is performed in Figure 11.

Based on the above proposed wind, solar and hydropower subsystems, a hybrid power system is finally obtained by MATLAB/Simulink [43], as shown in Figure 12.

3. Low Frequency Oscillation Response to Hybrid Regulation

The purpose of this section is to investigate the sensitivity of hydropower low frequency oscillations to regulation capacity of the grid-connected hybrid system. The hydropower governor, as an effective coordinator, has a powerful ability to control and adjust the power fluctuation as the intermittent wind and solar powers inject to the power grid. Based on its merit, the PID governor parameters, involving in the proportional adjustment coefficient (k_p), the integral adjustment coefficient (k_i) and the differential adjustment coefficient (k_d), are selected to analyze their regulation performance on low frequency oscillation of hydropower generation based on the Nyquist response and root-locus analysis.

3.1. Nyquist Response to PID Regulation

3.1.1. Nyquist Profile

The stability characteristic of Nyquist response is described as

R(s) = 1/K+ G(s)H(s)

(25)

where R(s), K, G(s) and H(s) represent equivalent function, adjustable gains, feedback path system function and forward path system function, respectively. According to the Nyquist encirclement property, the clockwise number around the point of −1/K equals to the difference between the closed-loop poles in the right half plane and the poles from G(s)H(s) plane. In other words, the stable condition for the closed-loop system is that there is no closed-loop pole in the right half plane, while in this situation the open-loop system can maintain stability. Thus, the Nyquist stability criterion is concluded as: To keep the system stability, the clockwise number around the point of −1/K equals to the negative pole numbers in the G(s)H(s) right half plane. If there is no poles in the G(s)H(s) right half plane, the clockwise number around the point of −1/K should be zero in order to keep the system stability.

3.1.2. Influences of Governor Parameters on Nyquist and Step Responses

The studied domains of the proportional adjustment coefficient (k_p), the integral adjustment coefficient (k_i) and the differential adjustment coefficient (k_d) are determined. In this case study, the wind, solar and hydropower ratio keeps 40 vs. 1 vs. 150 to benefits the identification of leading variate. The Nyquist and step response results are shown in Figure 13, Figure 14 and Figure 15.

From Figure 13a, there are similar behaviors for Nyquist results as the proportional adjustment coefficient (k_p) increases from 0.8 to 2.4. The Nyquist trajectory is a transversely zygomorphic closed curve starting at the origin, and it surrounds twice in the clockwise direction and then returns back to the origin. The curvature of Nyquist trajectory increases after it first go through the real-axis, but the curvature gradually decreases after the second real-axis traverse. All Nyquist trajectories intersect at the origin or at a point near the origin, whereas the peak response is continuously enhanced with the decrease of the proportional adjustment coefficient (k_p). Moreover, it is clearly observed in Figure 13b, there are some similarities for the step response of the hydropower system under a different proportional adjustment coefficient (k_p). For example, the step response within k_p = 2.4 first increases to the maximum overshoot (i.e., 0.014) and then gradually stabilized at the value of 0.002. However, notr that there is the smaller proportional adjustment coefficient, the smaller maximum overshoot. This reveals that the hydropower block has an excellent regulation capacity to the low frequency oscillation of the hybrid system within k_p = 0.8 compared with the situation within k_p = 2.4

The Nyquist response in Figure 14a experiences two clockwise rotations to return to the origin, meaning that the Nyquist trajectory is a transversely zygomorphic closed curve for the variation range of the integral adjustment coefficient k_i = [0.25, 1.25]. These curves only intersect at the origin, although the curvature of the first circle is greater than that of the second circle for a specific value of k_i. Additionally, regarding all different values of k_i, their second circles move right along with the real-axis and the corresponding peak response improved from k_i = 0.25 to k_i = 1.25. As performed in Figure 14b, all step responses have a decaying trend under different k_i conditions, differently, the peak value and final value of the step response increase gradually with the increase of integral adjustment coefficient (k_i). For instance, the peak value and final value for k_i = 1.25 are 0.0115 and 0.0081 in comparison with the relatively smaller values of 0.0098 and 0.001 for k_i = 0.25. However, fortunately, all overshoots for various parameter settings of k_i are less than the value of 10⁻³, thus the hydropower system can keep its stability during the selected variation domain of k_i.

Figure 15 shows that there is a similar change rule for both Nyquist and step responses for the change of the differential adjustment coefficient between k_d = 0.5 and k_d = 1.5. For a certain value of k_d, its Nyquist trajectory starts and also finally converges at the origin by turning clockwise twice; meanwhile, the curvature of the first circle is greater than that of the second circle. The peak response of various Nyquist trajectories increase as the differential adjustment coefficient rises from k_d = 0.5 to k_d = 1.5. Additionally, the step responses of the hydropower system perform decaying trends under various k_d conditions, and these response curves finally converge to a certain value, around 0.0017. The peak values for different values of k_d are extremely close to the mean value of 0.0098, although the peak value increases slightly with the change of the differential adjustment coefficient from k_d = 0.5 to k_d = 1.5. This implies that the studied value domain of the differential adjustment coefficient between (k_d) has a lower sensitivity to the low frequency oscillation of hydropower system.

3.1.3. Root-Locus Profile

To further study the pole and zero point characteristics of the low frequency oscillation stability of hydropower system for the solar/wind/hydropower hybrid system, the root-locus analysis is extensively conducted in this section. Herein, it is emphasized on that the zero point in the root-locus trajectory is the pole from the open-loop system while the pole in the root-locus trajectory equals to the closed-loop pole in this work. The detailed root-locus mechanism can be described as follows:

If the step response of the system is expressed as:

$$C\left(t\right)={\ell }^{-1}\left[C\left(s\right)\right]=A_0+{\displaystyle \sum _{i=1}^n{A}_i{e}^{s_{i}t}}$$

(26)

where C(s), A₀ and A_i denote the Laplace transformation of the step response, and residues corresponding to zero and pole, respectively. Hence, the general stability rule influenced by the pole and zero point in the root-locus trajectory is summarized two aspects: (1) The stability condition needs all closed-loop poles be located in the left half plane, meaning that all root-locus trajectories are in the left half plane, and (2) the rapidness of system response depends on the distance between the closed-loop pole (s_i = σ + jω) and the imaginary-axis. If such a distance is shorter, there exists a rapider decaying speed of the e^sit and thus results in an excellent rapidness of system response.

3.1.4. Influences of Governor Parameters on Root-Locus Response

To clearly perform the PID parameters regulation on the root-locus response, the control variate method is adopted. That is, two of the PID parameters remain unchanged, but only one of them changes gradually. Besides, the study is conducted based on the wind, solar and hydropower quota 40:1:150. The root-locus results are shown in Figure 16, Figure 17 and Figure 18.

Both zero points and poles (Figure 16a) lie in the real-axis under the condition of k_p = 1.6, k_i = 0.75 and k_d = 0.5, and the zero point coordinates are (−0.0925, 0) and (−1.69, 0) as well as the pole coordinates are (−1, 0), (−1.05, 0) and (−15, 0). There are three root-locus responses: (i) The first root-locus trajectory begins at the pole (−1, 0) and ends at the zero point (−0.0925, 0). (ii) The second root-locus trajectory starts from the pole (−1.05, 0) and is extended towards the zero point (−1.69, 0), and (iii) the third root-locus trajectory stems from the pole (−15, 0) and goes to the negative infinity. All root-locus trajectories are located in the left half plane, and thus the hydropower system is stable in this situation. With the increase of the differential adjustment coefficient from k_d = 0.8 to k_d = 2.4, the pole coordinate (Figure 16b) remains unchanged, but the zero point is gradually closer to the imaginary-axis. This indicates that there is an excellent rapidness response of the hydropower system to the disturbance of low frequency oscillations within a smaller k_p.

From Figure 17a, it is found that the all poles and zero points are in the real-axis under the condition of k_p = 1.6, k_i = 0.75 and k_d = 1. The zero point coordinates include (−0.255, 0) and (−1.84, 0), and the pole coordinates involve in (−1, 0), (−1.1, 0) and (−15, 0). There are three root-locus trajectories: the first trajectory is between the pole (−1, 0) and zero point (−0.255, 0), the second trajectory is from the pole (−1.1, 0) to the zero point (−1.84, 0), and the last trajectory starts from the pole (−15, 0) and links to the negative infinity. The hydropower system operates safely and stably since these root-locus trajectories are in the left halt plane. Figure 17b illustrates the regulation capacity of the integral adjustment coefficient from k_i = 0.25 to k_i = 1.25. With the decrease of parameter k_i, the pole coordinate keeps constant while the zero point moves closely to the imaginary-axis. This indicates that the hydropower system has a better rapidness as the parameter k_i increases to 1.25.

Similarly, the poles and zero points are in the real-axis for k_p = 1.6, k_i = 0.75 and k_d = 1 in Figure 18a. The zero points coordinates are (−0.255, 0) and (−1.84, 0), while the pole coordinates are (−1, 0), (−1.1, 0) and (−15, 0). The first root-locus trajectory is from the pole (−1, 0) to the zero point (−0.255, 0), the second root-locus trajectory is from the pole (−1.1, 0) to the zero point (−1.84, 0), and the last root-locus trajectory experiences from the pole (−15, 0) to negative infinity. All the three trajectories lie in the left half plane, resulting in a steady operational condition of the hydropower system responding to the low frequency oscillations. Moreover, as shown in Figure 18b, the zero point has an opposite trend with the pole with the increase of differential adjustment coefficient (k_d), where the zero point moves closely to the imaginary-axis, but the pole is far away from the imaginary-axis. This reveals that there is a rapid response of the hydropower system to the low frequency oscillation within k_d = 1 in comparison to the condition of k_d = 1.6.

4. Low Frequency Oscillation Response to Renewable-Quota

This section focuses on investigating the influence of different wind/solar/hydropower quota (i.e., W:S:H) on the low frequency oscillation mode and system stability. Simultaneously, the line distance ratio between the solar-load transmission and the wind-load transmission (i.e., SL:WL) is also considered. Based on this, the wind/solar/hydropower quota include 20:1:150, 30:1:150 and 40:1:150 (the combination is ratio of installed capacity), and the line distance are determined as 1:1, 2:1 and 3:1. These renewable quotas are selected in this paper to study the evolution law between the increase of installed capacity of wind power and the low-frequency oscillation of power grid under different renewable quotas. Similarly, the transmission line length ratios are selected in this paper to study the influence of the solar-load line on the low frequency oscillation in the hydropower-based hybrid renewable energy system. The low frequency oscillation mainly be represented by the three-phase parallel RLC load, the grounding transformer, the three-phase PI section line and the wind-farm transmission line. The results are demonstrated in

Table 1. Low frequency oscillation response of hydropower system to wind/solar/hydropower quota and line distance ratio.

Line Length Ratio (S:W)	W:S:H	Eigenvalue	Frequency/HZ	Damping Ratio/%	Remark
1:1	20:1:150	−1.14 × 10−1 + j0	0.0182	100	Three-Phase Parallel RLC Load
		−1.26 × 101 + j0	2.01	100	Grounding Transformer
		−9.18 × 102 − j4.47 × 103	7.26 × 102	20.1	Three-Phase PI Section Line
		−3.15 × 103 − j9.43 × 103	1.58 × 103	31.7	Windformsystem/30 km line
	30:1:150	−1.13 × 10−1 + j0	0.018	100	Three-Phase Parallel RLC Load
		−1.65 × 101 + j0	2.63	100	Grounding Transformer
		−9.47 × 102 − j4.61× 103	7.48 × 102	20.1	Three-Phase PI Section Line
		−3.15× 103 − j9.43× 103	1.58 × 103	31.7	Windformsystem/30 km line
	40:1:150	−1.13 × 10−1 + j0	0.018	100	Three-Phase Parallel RLC Load
		−1.99 × 101 + j0	3.17	100	Grounding Transformer
		−9.77 × 102 − j4.74 × 103	7.71 × 102	20.2	Three-Phase PI Section Line
		−3.15 × 102 − j9.43 × 103	1.58 × 103	31.7	Windformsystem/30 km line
2:1	20:1:150	−1.14 × 10−1 + j0	0.0182	100	Three-Phase Parallel RLC Load
		−1.26 × 101 − j0	2.01	100	Grounding Transformer
		−4.62 × 103 − j0	7.36 × 102	100	Three-Phase PI Section Line
		−3.23 × 102 − j7.11 × 103	1.13 × 103	4.54	Windformsystem/30 km line
	30:1:150	−1.13 × 10−1 + j0	0.018	100	Three-Phase Parallel RLC Load
		−1.65 × 101 + j2.37 × 10−10	2.63	100	Grounding Transformer
		−9.67 × 102 − j4.58 × 103	7.45 × 102	20.7	Three-Phase PI Section Line
		−3.22 × 102 − j7.11 × 103	1.13 × 103	4.53	Windformsystem/30 km line
	40:1:150	−1.13 × 10−1 + j0	0.018	100	Three-Phase Parallel RLC Load
		−1.99 × 101 + j0	3.17	100	Grounding Transformer
		−9.98 × 102 − j4.71 × 103	7.66 × 102	20.7	Three-Phase PI Section Line
		−3.22 × 102 − j7.11 × 103	1.13 × 103	4.52	Windformsystem/30 km line
3:1	20:1:150	−1.14 × 10−1 + j0	0.0182	100	Three-Phase Parallel RLC Load
		−1.26 × 101 + j0	2.01	100
		−8.82 × 102 − j4.39 × 103	7.13 × 102	19.7	Three-Phase PI Section Line
		−5.85 × 102 + j627 × 103	1.00 × 103	9.28	Windformsystem/30 km line
	30:1:150	−1.13 × 10−1 + j0	0.018	100	Three-Phase Parallel RLC Load
		−1.65 × 101 + j0	2.63	100	Grounding Transformer
		−1.72 × 102 − j4.48 × 103	7.13 × 102	3.85	Three-Phase PI Section Line
		−5.85 × 102 + j6.27 × 103	1.00 × 103	9.28	Windformsystem/30 km line
	40:1:150	−1.13 × 10−1 + j0	0.018	100	Three-Phase Parallel RLC Load
		−1.99 × 101 − j0	3.17	100	Grounding Transformer
		−1.76 × 102 − j4.47 × 103	7.12 × 102	3.95	Three-Phase PI Section Line
		−5.85 × 102 + j6.27 × 103	1.00 × 103	9.28	Windformsystem/30 km line

.

With respect to the line distance ratio between the solar-load transmission and the wind-load transmission (i.e., SL:WL) is 1:1, the evaluation results (Table 1) show that the oscillation frequency of the three-phase parallel RLC load and the wind-farm transmission line roughly keeps unchanged as the wind/solar/hydropower quota (i.e., W:S:H) changes from 20:1:150, 30:1:150 to 40:1:150. Conversely, the oscillation frequency of the grounding transformer and the three-phase PI section line continuously increases in this situation, for instance, the oscillation frequency of grounding transformer is 2.01 vs. 2.63 vs. 3.17, and the oscillation frequency of three-phase PI section line is 7.26 × 10² vs. 7.48 × 10² vs. 7.71 × 10². Besides, the damping ratios for these four stability indicators are constant. This suggests that the appropriate increase of the wind farm installed capacity is beneficial for the electricity production without system stability defects. As regards the condition of SL:WL = 2:1, the oscillation frequency of the three-phase parallel RLC load and the wind-farm transmission line is unchanged while the oscillation frequency of the grounding transformer and the three-phase PI section line obviously increases with the capacity-increasing improvement of wind farm. The oscillation damping ratio of the three-phase parallel RLC load and the grounding transformer has no variation in contrast to the decline damping ratio of the three-phase PI section line and the wind-farm transmission line. It is indicated that the wind farm installed capacity is unable to design too large, otherwise the stability problem is likely to occur under the condition of SL:WL = 2:1. Regarding the condition of SL:WL = 3:1, there is no variation for the oscillation frequency of the three-phase parallel RLC load, the three-phase PI section line and the wind-farm transmission line, but the oscillation frequency of the grounding transformer rises. Differently, the oscillation damping ratio of the three-phase parallel RLC load, the grounding transformer and the wind-farm transmission line totally remains unchanged, whereas the oscillation damping ratio of the grounding transformer first declines and then rises. Thus, it is possible to increase the installed capacity of the wind farm since the hydropower generation has a low sensitivity to the low frequency oscillations in this case.

In addition, as the line distance ratio between the solar-load transmission and the wind-load transmission increases from 1:1 to 3:1, the oscillation frequency and damping ratio of the three-phase parallel RLC load and the grounding transformer are generally unchanged; however, such indicators of the three-phase PI section line has an increase trend and then gradually declines. Additionally, the oscillation frequency of the wind-farm transmission line decreases, although its damping ratio decreases first and then maintains an increase trend. In other words, when the renewable-quota is determined, the system stability and its dynamic behaviors are able to be enhanced by increasing the length of the solar-load transmission line.

5. Conclusions

Renewable energies, especially for wind and solar energies, are major electricity generation strategies to address the constraints of resources and the environment for energy development. However, the intermittency, uncontrollability and randomness demerits of wind/solar generations pose large challenges for the stable integration into the power grid. Based on this, a hybrid system of the hydropower integrating with wind and solar sectors is presented in order to investigate the low frequency oscillation stability of hydropower systems responding to the fluctuation of wind and solar power injections. The main adopted approaches are Nyquist response and root-locus analysis to evaluate the regulation capacity of governor parameters (i.e., the proportional adjustment coefficient k_p, the integral adjustment coefficient k_i, and the differential adjustment coefficient k_d) to the system stability profile. Furthermore, the influence of different renewable-quota and transmission line distance ratios on the low frequency oscillation mode of the hydropower system is quantified. The major conclusions can be drawn as follows.

• Under the case where the wind, solar and hydropower ratio is 40:1:150, it is interesting that the smaller the governor parameters (k_p, k_i, and k_d), the smaller the Nyquist overshoot and step fluctuation. Herein, the studied domains for k_p, k_i, and k_d are [0.8, 2.4], [0.25, 1.25] and [0.5, 1.5], and thus the optimal values for maximally reducing hydropower low frequency oscillation are finally determined as k_p = 0.8, k_i = 0.25 and k_d = 0.5.
• The wind/solar/hydropower hybrid system keeps global stability in the studied governor parameter domains since the Nyquist and root-locus low frequency oscillation responses meet the relevant stability criteria, i.e., the clockwise number around the point of −1/K equals to the negative pole numbers in the G(s)H(s) right half plane, as well as all root-locus trajectories are in the left half plane. Despite this merit, the overshoot problem is expected to arouse great attention and discussion to reduce the fatigue damage of hydropower components.
• Aiming at different wind/solar/hydropower quotas (i.e., 20: 1: 150, 30: 1: 150, and 40: 1: 150), the four quantified indicators (i.e., the three-phase parallel RLC load, the grounding transformer, the three-phase PI section line, and the wind-farm transmission line) show that it is beneficial to increase the wind farm installed capacity to maximize the electricity production without system stability defects under the solar-load and wind-load line ratios for 1:1 and 3:1 excepting for 2:1. This is contributed by the smaller quantified values of oscillation frequency and damping ratio.
• Regarding a certain wind/solar/hydropower quota, it is a promising strategy to increase the solar-load transmission line in order to achieve the safe and stable operation of the hybrid system and a relatively excellent dynamic regulation capacity of the hydropower governor.

Hydropower systems, such as the vital energy storage system and power balance controller, play a key role towards the successful achievement of the world’s long-term clean energy goals. For this, the results provide important guidance for realizing the maximization of economic and environmental benefits in a hybrid power grid. The proposed solutions evaluate the safety situation of the hybrid system in quantitative terms, which is cheered by the energy stakeholders of governments, energy planers, regulation entities, storage companies and environmental organizations.