A Short-Term Optimal Scheduling Model for Wind-Solar-Hydro-Thermal Complementary Generation System Considering Dynamic Frequency Response

This paper proposes a model to realize the coordinated optimal dispatch of wind-solar-hydro-thermal hybrid power generation system, aiming at minimizing the power generation cost of thermal generators and maximizing the water storage value of hydropower stations at the end of the scheduling periods, while considering the dynamic frequency response of wind/solar/hydro/thermal generators. Considering the virtual inertia and droop control of wind farms and PV stations, the dynamic frequency response model of wind-solar-hydro-thermal multi-energy complementary system is derived and the metrics that evaluate the frequency dynamic characteristics of the generation system are presented. Then the dynamic frequency response constraints are incorporated into the traditional optimal scheduling model and the Mixed Integer Linear Programming (MILP) method is used to solve it. Finally, the validity and applicability of the proposed model are verified by simulation examples.


Sc, c
Set and index of scenarios. i Index of reserve types and i ∈ {1, 2, . . . , 6} represents 30s real-time response reserve, AGC reserve, 10min spinning reserve, 30min operating reserve, 60min operating reserve and cold reserve, respectively.
B. PARAMETERS π c The probability of scenario c. a g , b g , c g Operating cost coefficient of thermal generator g. PL t System load at period t. L l, max Maximum transmission power of line l. G l,m Generation shift distribution factor of node m to line l. P h Maximum output of hydro generator h. P h Minimum output of hydro generator h. P g Maximum output of thermal generator g. P g Minimum output of thermal generator g.

RU g
Ramp-up limit of thermal generator g. RD g Ramp-down limit of thermal generator g. UT g Minimum online time for thermal generator g. DT g Minimum offline time for thermal generator g.

I. INTRODUCTION
The development and utilization of clean energy such as wind, solar, and hydropower is the main direction for the transformation of the energy structure of the power system. With increasing penetration of renewable energy sources, the planning and operation of the generation system has become a research hotspot [1]- [4]. Wind and solar power systems are connected to the grid through inverters, which are decoupled from the system frequency, making it difficult to provide rotational inertia. As a result, the large-scale integration of wind and solar power will reduce the system inertia, weaken the frequency modulation capability, and threaten the stability of the system [5]. System inertia plays a vital role in suppressing frequency fluctuations. The decrease of the inertia will increase the system frequency deviation and even induce vicious events such as frequency collapse. Many scholars have carried out research to solve the above-mentioned problem. The consideration of dynamic frequency response in the traditional scheduling model has received extensive attention. Ahmadi and Ghasemi [6] added the nadir frequency constraint to the unit commitment model. By optimizing the number of online synchronous units, the inertial response and Primary Frequency Response (PFR) of the system can be improved. Muzhikyan et al. [7] further considered the constraint of Rate of Change of Frequency (RoCoF). Teng et al. [8] added the quasi steady-state frequency constraints, and the influence of the pumped-hydro storage generators on the PFR was considered. Wen et al. [9], [10] regarded the battery energy storage as an effective regulation source that quickly responds to system frequency deviations. Trovato et al. [11] and Badesa et al. [12] introduced enhanced frequency response constraint, which responds to frequency changes faster than PFR. The RoCoF, nadir frequency and steady-state frequency constraints were established to keep the frequency within an acceptable range after a fault. Paturet et al. [13] studied the stochastic unit commitment in low-inertia grids, and the corresponding frequency constraints were derived from the uniform system frequency dynamic model. However, the above studies only considered the frequency response characteristics of thermal generators.
Akbari and Madani [14]- [16] pointed out that wind turbines can simulate the inertia and droop control of conventional units, which will change the frequency response characteristics of the original system. Considering the virtual inertia and PFR, Yan and Saha [17] proposed an improved Average System Frequency (ASF) model. However, the primary frequency regulation response was described by a static model, which only reflects the magnitude of the steadystate frequency deviation. Wogrin et al. [18] combined the synchronous inertia and virtual inertia of wind farms to study the influence of inertia on power generation expansion planning, but only the minimum inertia demand constraint of the system was considered. Zhang et al. [19] incorporated the frequency regulation characteristics of thermal generators and the frequency support of renewable energy sources in the model. Finally, the analytical formulation for system frequency nadir was derived. Ge et al. [20] considered the virtual inertia and droop control of wind farms in the optimization model. Through time-frequency domain analysis, the whole process constraint of ''inertia support-nadir frequency point-quasi steady-state-secondary frequency modulation'' was achieved.
To cope with the low inertia of solar power system, most researches improve the system frequency response ability by changing the control strategy of the virtual synchronous machine. Xin et al. [21] studied the system control method to make PVs can operate in the droop mode and provide primary frequency support for the power system. Hou et al. [22] studied the influence of virtual inertia and droop control on system frequency response performance, and proposed an adaptive virtual inertia control strategy that can quickly suppress system frequency fluctuations.
However, the above studies only focus on wind-thermal complementary generation systems or solar generation systems, with a single application scenario. Hydropower has many advantages in system operation, such as short startup time, strong flexibility and energy storage. But the research on the dynamic frequency response of system including hydropower generators has not attracted much attention.
In this context, this paper proposes a short-term optimal scheduling model for wind-solar-hydro-thermal complementary generation system considering the dynamic frequency response of wind/solar/hydro/thermal generators. The main contributions of this paper are: 1) Taking into account the synchronous inertia of hydro and thermal generators, the virtual inertia and droop control of wind farms and PV stations, the frequency response equivalent model of the complementary generation system is proposed. The constraints of RoCoF, nadir frequency deviation and quasi steady-state frequency deviation are derived.
2) A mixed integer linear programming model for windsolar-hydro-thermal complementary optimal scheduling is established. The simulation results show the validity of the consideration of dynamic frequency response. In addition, the applicability and accuracy of the proposed model are further verified.
The remainder of the paper is organized as follows. Section II models the dynamic frequency response of the system as well as the corresponding dynamic frequency metrics. Section III presents the short-term optimal scheduling model in detail, including objective function, complex constraints and solution approach. Section IV verifies the validity and applicability of the proposed model through case studies. Section V presents the conclusions.

II. SYSTEM DYNAMIC FREQUENCY RESPONSE MODEL
Frequency response describes the behavior of the power system after a sudden loss of a large synchronous unit or power plants. In order to ensure the frequency stability and avoid system collapse, the system must meet the corresponding dynamic frequency metrics [13]: RoCoF, nadir frequency and quasi steady-state frequency.

A. EQUIVALENT MODEL OF SYSTEM FREQUENCY RESPONSE
In this section, taking into account the virtual inertia and droop control of wind farms and PV stations, the frequency response equivalent model of wind-solar-hydro-thermal coordinated generation system is established.
When wind turbines with inertial response and droop control capabilities are connected to the grid, the equivalent model of wind farms participating in frequency regulation is shown in FIGURE 1 [20]. The virtual inertia control of wind farms is a transient process, mainly used to prevent rapid changes of frequency. Droop control is a steady-state process, mainly used to reduce system frequency deviation. The superposition of the above two control methods can achieve fast frequency response and active power modulation.
Due to the fact that both PV stations and wind farms are connected to the grid through inverters, their equivalent model structures are the same.
The simplified equivalent model of ideal hydraulic turbines and governors is shown in FIGURE 2 [23]. The simplified equivalent model of ideal steam turbines and governors is shown in FIGURE 3 [13]. FIGURE 4 shows the equivalent model of system frequency response. Aggregate mechanical power gains K G , K H , K W and K V represent the frequency regulation capabilities of thermal generators, hydro generators, wind farms and PV stations, respectively.
According to FIGURE 4, the transfer function of the system can be obtained as shown in (1). Detailed results are  shown in appendix A.
When there is a power outage in the system, such as the loss of a generator, this process can be regarded as a step response, that is P e = − P s , and we can get: Perform the inverse Laplace transform of (2) and Appendix A shows the specific process. The time domain expression of the system frequency response is derived as:

B. ANALYTIC DERIVATION OF FREQUENCY METRICS
On the basis of (3), the dynamic frequency metrics are derived: RoCoFḟ max , the nadir frequency deviation f max and quasi steady-state frequency deviation f ss .

1) ROCOF
When the system is disturbed, df dt is defined as the initial frequency drop rate, which is mainly used to reflect the inertial response capability of the system.
At t = 0, the frequency drops the fastest. Take derivative of (3) and set t = 0 to obtain the RoCoF: The nadir frequency of the system reflects the degree of frequency drop after the system is disturbed. It is an important indicator for formulating the under-frequency load shedding rules of the system. At the frequency nadir, d f (t) dt = 0 is satisfied. Then, the time to reach the nadir is: Therefore, the nadir frequency deviation is:

3) QUASI STEADY-STATE FREQUENCY DEVIATION
When the frequency reaches the steady state, d f (t)

dt
= 0 and f (t) = a (a is a constant frequency) are satisfied. According to the final value theorem, we can get f ss = lim s→0 [sF (s)]: When a power disturbance occurs, the above frequency metrics should satisfy the following constraints [13]: The first inequality in (7) indicates that the system frequency change rate cannot exceed the maximum permissible RoCoF; the second inequality indicates that the nadir frequency deviation must be within the allowable range, otherwise it will trigger low-frequency load shedding; the last inequality indicates that the quasi steady-state frequency deviation cannot exceed the deviation limit.

III. SHORT-TERM OPTIMAL SCHEDULING MODEL FOR WIND-SOLAR-HYDRO-THERMAL COMPLEMENTARY GENERATION SYSTEM
The short-term optimal scheduling model for wind-solarhydro-thermal complementary generation system considering dynamic frequency constraints is as follows.

A. OBJECTIVE FUNCTION
In order to promote energy conservation and avoid resource waste, the model aims to minimize power generation costs and maximize water storage value at the end of the scheduling periods. The cost of power generation includes the operating cost, startup cost and shutdown cost of thermal generators [24], [25].
The first term in (8) represents the power generation cost of thermal generators. The second item is the water storage value at the end of the scheduling periods. It represents the value of the electric energy that the reservoir water volume can convert in the future, which is determined by the unit value of the water storage, the current water volume of reservoirs and the hydroelectric conversion rate of its own and downstream reservoirs.
c: WATER BALANCE CONSTRAINT The hydraulic connection between the cascade hydropower stations causes the time delay of the flow of water, which lead to an imbalance in the water volume during the scheduling periods and need to be considered in the model. Therefore, the equation of water balance is expressed as:

e: RESERVE CAPACITY CONSTRAINTS
The reserve capacity provided by the hydro generators needs to satisfy the available capacity constraints.
2) WIND AND SOLAR POWER CONSTRAINTS Due to existence of prediction errors, it is not practical to arrange wind and solar power output plans based on forecast values. Equations (18) and (19) indicate that the sum of forecast values and forecast errors is regarded as the maximum power that can be dispatched during the period. Equations (20) and (21) express the constraints on the planned output of wind and solar power.
3) THERMAL POWER CONSTRAINTS a: THERMAL GENERATOR OUTPUT CONSTRAINT U g,t P g ≤ P c g,t ≤ U g,t P g (22) b: RAMP CONSTRAINTS c: MINIMUM ONLINE/OFFLINE TIME CONSTRAINTS

4) SYSTEM CONSTRAINTS a: POWER BALANCE CONSTRAINT
c: RESERVE DEMAND CONSTRAINTS Sufficient reserve capacity should be provided by hydro and thermal generators to remunerate the prediction error of wind and solar power, so as to ensure the stability of the system frequency.
According to the response time and frequency band, the reserve can be divided into 30s real-time response reserve, AGC reserve, 10min spinning reserve, 30 min operating reserve, 60 min operating reserve and cold reserve. In order to realize the reasonable allocation of resources, the discrete Fourier transform method in [26] is employed to quantify the demand of various types of reserves.

d: DYNAMIC FREQUENCY RESPONSE CONSTRAINTS
The constraint (7) derived from the equivalent model of system dynamic frequency response contains a large number of parameters. In order to reduce the calculation burden, the parameters can be aggregated.  (49) Equations (32)

C. SOLUTION APPROACH
The proposed model is a complex nonlinear optimization model, which can be transformed into a mixed integer linear programming model and solved by the commercial optimization software CPLEX.
The operating cost of the thermal generators, minimum online/offline time constraints and the hydropower output function in the hydropower constraints are linearized using the methods in [24] and [27], respectively. The following will focus on the linear conversion method for the dynamic frequency response constraints.
Equation (35) can be transformed into: The right side of the (50) is linear and the left side is the continuous variable K W ,t multiplied by the binary variables I h,t and U g,t respectively. Auxiliary variables KWI h,t and KWU g,t are introduced to represent the product of these two variables respectively. Thus, (50) can be expressed as: Equations (36)-(44) and (48) are all nonlinear constraints, but their linearization methods are the same with (35). Therefore, to save space, the linearization process of (36)-(44) and (48) are omitted here. Furthermore, the nadir frequency deviation constraint expressed by (49) is linearized by the method proposed in [7].

IV. SIMULATION RESULTS
In order to verify the validity of the proposed model, analysis and calculations are carried out with the combination of ten thermal generators, two cascade hydropower stations consisting of seven hydropower generators, a wind farm with an installed capacity of 1000MW and a PV station with an installed capacity of 900MW.
The parameters of thermal generators are revised from the IEEE 10-unit 39-bus test system given in [24]. Refer to Table 5-8 in the appendix for the data of thermal generators and hydro stations. It is assumed that the unit value of water storage in the future is 1$/MWh, the hydroelectric conversion rate of hydro stations is 0.2 MWh/10 4 m 3 , and the water flow delay time is 1 hour. Forecasted value of load, wind power and solar power are shown in FIGURE 5. The base frequency is 50Hz, the maximum permissible RoCoF is 0.2Hz/s, the under-frequency load shedding trigger is 49.6Hz and the maximum permissible steady-state frequency deviation is 0.3Hz. The system power outage in each period is shown in appendix TABLE 9. The scheduling model is solved by CPLEX12.5 on a personal computer with Intel Core i7 3.0GHz CPU and 4GB RAM.

A. VALIDITY OF CONSIDERING DYNAMIC FREQUENCY RESPONSE CONSTRAINTS
To verify the validity of the proposed model, the comparison of Method 1 and Method 2 has been analyzed.
Method 1: Model solution without considering dynamic frequency response constraints.
Method 2: Based on Method 1, with the dynamic frequency response constraints considered.    6 and FIGURE 7 show the units state of the two methods. Comparing Method 2 with Method 1, we can find that in response to the frequency change caused by the power outage, additional thermal generators will be online to improve the system's ability of resisting frequency stability risks. The dynamic frequency metrics of the two methods are analyzed below.

1) ROCOF
The initial frequency change rate of system is mainly affected by the inertia. FIGURE 8 compares the RoCoF in Method 1 and Method 2.
From FIGURE 8, we find that when considering frequency response constraints, each period does not exceed the maximum permissible RoCoF (0.2Hz/s). When the frequency constraints are not considered, except for periods 19 and 21, VOLUME 9, 2021 other periods cannot meet the frequency change rate requirement. Taking period 11 as an example, in Method 1, thermal generator G1 and G2 are online, and the total inertia of system is 7s, which cannot meet the minimum inertia requirement of 15.25s. Therefore, the system frequency drops rapidly and the frequency stability is threatened. When Method 2 is employed, the thermal generators G3, G7 and hydro generator H4 are added to be online during this period. The total inertia reaches 21s, which improved the inertia response capability of the system and effectively prevented the frequency drop rate from being too fast.

Comparison of nadir frequency deviation in Method 1and
Method 2 is given in FIGURE 9. From FIGURE 9, when the frequency response constraints are considered, the nadir frequency deviation of each period is less than 0.4Hz, which does not reach the under-frequency load shedding trigger frequency of 49.6Hz. While in Method 1, the frequency deviation exceeds the maximum allowable deviation value of 0.4Hz during almost half of the periods. This shows that the model proposed in this paper reduces the risk of system load shedding and effectively improves the frequency regulation ability of the system.

3) QUASI STEADY-STATE FREQUENCY DEVIATION
Comparison of quasi steady-state frequency deviation in Method 1 and Method 2 is given in FIGURE 10. As can be seen from FIGURE 10, in Method 2, the quasi steady-state frequency deviation can be controlled within the allowable deviation range. However, in Method 1, the steady-state frequency deviation of all periods exceeds the maximum permissible steady-state frequency deviation (0.3Hz), which shows the importance of considering the dynamic frequency constraints.  We can see from FIGURE 11 and TABLE 1 that: i) When there is a power outage in the system, the number of online generators in Method 2 increases. Therefore, the overall output of thermal generators has increased, but the increase is small. The reason is that in order to reduce the cost of the system, the newly added units only operate at the minimum output level, thereby providing inertia support. As a conclusion, when considering the dynamic frequency response constraints, the value of the objective function increased by 20.34%. In addition, due to the balance of system supply and demand, hydropower output and the acceptance level of wind and solar power have decreased. It can be seen that the abandoned wind and solar power in Method 2 is approximately 1.38 times that of Method 1.
ii) When in dry season, hydropower generators are mainly online during periods of high load demand and low wind and solar power output. Affected by the small amount of natural inflow water, the initial and final reservoir volume was the same and has no value of water storage.

B. THE IMPACT OF WET/DRY SEASON ON SYSTEM SCHEDULING
The above section is the analysis of the system scheduling results of the dry season. This section will explore the influence of the wet season on the scheduling model.
Method 3: Based on Method 2, with the natural inflow of hydro stations meets the characteristics of wet season. The unit state is shown in FIGURE 12 and TABLE 2 shows the optimization results of this method. Comparing the optimization results of Method 3 with Method 2, we can get the following conclusions:  i) In wet season, the proportion of thermal generators in online operation units is relatively high. However, due to the natural inflow water is abundant in wet season, the dynamic frequency response requirements of the system can be met by the online hydropower generators. Moreover, the maximum volume water of reservoirs is reached at the end of scheduling period, so that the value of water storage in the wet season is greatly improved.
ii) The advantages of hydropower are more obvious in wet season. The hydropower units not only increased its output, but also assumed the main task of providing reserve, which can reduced the output of thermal generators and bring desirable economic benefits. Moreover, hydropower is clean, and it also brings better environmental benefits.

C. ANALYSIS OF MODEL APPLICABILITY
Generally, the minimum startup and shutdown time of thermal generators are relatively long. To verify the applicability of the model in this paper, on the basis of Method 2, the minimum startup and shutdown parameters of the thermal generators are changed. The specific parameters settings are shown in Table 3.
Method 4: Based on Method 2, with the minimum startup and shutdown time of thermal generators increased.
For the comparison of units state, RoCoF, nadir frequency deviation and quasi-steady state frequency deviation in case of whether frequency response constraints are considered, see FIGURE 13-17 in the appendix.
According to the analysis of the scheduling results, we can find that after changing the parameters of the thermal generators, regardless of whether the dynamic frequency response constraints are considered or not, the nadir frequency deviation do not exceed 0.4Hz. But when considering the frequency response, the frequency deviation value is smaller. Regarding the RoCOF and quasi-steady state frequency deviation, when the dynamic frequency constraints are not considered, there will always be several periods exceeding the maximum limit, which will threaten the safe and stable operation of the system.
It can be seen that the model proposed in this paper can reasonably arrange the unit commitment. For the generators with a long continuous startup and shutdown time in practice, it still has the advantages of improving the inertial response and frequency regulation of the system.

D. ANALYSIS OF MODEL ACCURACY
When the dynamic frequency response is not considered in Method 1, the model has 30932 constraints, 79691 variables. For the remaining three methods with dynamic frequency constraints, the model has 106580 constraints and 81851 variables. In practice, the operators should make a quick decision, so we set a time limit of 10min as the stopping criterion. The optimization results obtained from the four methods are shown in Table 4, where the final gap is expressed as: The final gaps of Method 1-4 are 0.17%, 0.58%, 1.21% and 1.45%, respectively. The results are quite acceptable because of the large number of variables in the models. The final gap means the objective solution is at least within the percentage of the global optimal solution. The longer run time may only tighten the bounds but not improve the objective solution [28], [29].

V. CONCLUSION
Aiming at the power system with high penetration of wind and solar power, this paper fully considers the impact of its low inertia and proposes a short-term optimal dispatch model for wind-solar-hydro-thermal coordinated generation system considering dynamic frequency response.
Taking into account the virtual inertia and droop control of wind farms and PV stations, the dynamic frequency response model is established, which can describe the dynamic change process of frequency under disturbance in detail. The derived frequency constraints are added to the scheduling model. The simulation results show that by rationally arranging the unit commitment, the rate of change of frequency after disturbance can be effectively reduced and the frequency response capability of the system can be enhanced.
The comparison of optimization results of wet and dry season shows that the advantages of hydropower are more obvious in wet season. In wet season, only the hydropower generators are online to meet the frequency response requirements of the system, which effectively reduces the cost of system power generation. Meanwhile, the objective function takes into account the future value of water storage at the end of scheduling period, which is of great significance for improving the utilization of clean energy.
The model proposed in this paper is still applicable when the thermal generators with a long continuous startup and shutdown time. Moreover, the accuracy of the model is also guaranteed while ensuring the solution speed.
With the continuous development of energy storage technology, how to coordinate the energy storage and the output of conventional units to effectively improve the inertial response characteristics of the system is the next issue to be studied.