# Characteristic Analysis and Optimal Regulation of Primary Frequency Regulation Condition in Low Water Head Area Based on Hydraulic-Mechanical-Electrical Coupling Model of Pumped Storage Unit

Authors: Cong Lv; Yanhe Xu; Xin Wu; Qing Zhang
Journal: Complexity
 (2020)

---

## Abstract

Pumped storage power station is an important regulating tool for peak load regulation and frequency regulation of the power grid, especially its primary frequency regulation function, which is of profound significance for the safety and stability of the power grid. As the core equipment of the pumped storage power station, the reversible design of the pump turbine makes it easy to have hydraulic fluctuation and mechanical instability when the unit runs to the “S” characteristic area, which will cause the frequency oscillation of the generator under the condition of primary frequency regulation. Therefore, some innovative work is studied in this paper: (1) an accurate hydraulic-mechanical-electrical coupling system mathematical model of pumped storage unit regulation system (PSURS) is established based on full characteristic curve of the pump turbine and seventh-order synchronous generator and excitation system; (2) the dynamic response characteristics of primary frequency regulation of pumped storage unit (PSU) under different water heads and different frequency disturbances are analyzed by numerical simulation; (3) in view of the unstable condition of the unit under the large frequency disturbance when it operates in the low head and small load area, the objective optimization function considering the ITAE index of hydraulic, mechanical, and electrical factors is proposed; and (4) fractional-order PID controller and the bacterial-foraging chemotaxis gravitational search algorithm (BCGSA) combined optimization strategy is used for PSURS optimization regulation and parameter optimization. The results show that the joint optimization strategy proposed in this paper has smaller objective function value, and makes the PSURS pass through the unbalanced area quickly, with better primary frequency regulation speed and smaller regulation depth.


---

## Body

## 1. Introduction

In recent years, with the rapid growth of China’s power load, the demand for the stability of the power grid is rapidly increasing, and the shortage of peak load regulation capacity has become an important problem restricting the development of the power system. Therefore, the combined operation of pumped storage technology and new energy is widely used in power system [1, 2]. Pumped storage power station is mainly responsible for peak and frequency regulation and peak and valley cutting, which can improve the power supply quality, flexibility, and reliability of the power grid and enable various power resources to be reasonably applied, which is an important part of the power system [3]. The pumped storage power station is dedicated to power management and stability regulation of grid and isolated power system [4].The frequency regulation function of PSU is very important to maintain the stability of power grid frequency, especially for those power grids with small proportion of prompt adjustment units [5]. An et al. proposed an improved frequency dead zone with feed-forward control to optimize the frequency regulation function of PSU [6]. Yang and Yang and Zhao and Ren also studied the role of pumped storage units in enhancing power system stability and frequency control [7, 8]. It can be seen from some accidents in the power system that the units with the input of primary frequency regulation function and the correct setting of primary frequency regulation parameters can greatly restrain the expansion of the accident and ensure the stability of the power system [9–11]. The pumped storage unit regulation system (PSURS) has strong nonlinear and time-varying characteristics [12], and PSU has “S” characteristic area, which leads to the frequency oscillation of the unit under low water head [13, 14]. Therefore, it is of profound significance to optimize the dynamic performance and control quality of PSURS under the condition of low water head [15].At present, PID controller is widely used as controller in hydraulic turbine regulating system due to its good control effect, simple and easy to implement structure, strong robustness, and other advantages [16, 17]. However, when the PSURS is running at a low water head, it is easy to enter the “S” characteristic area, resulting in widespread range fluctuation and unstable operation of the unit’s frequency. Especially when the unit participates in primary frequency regulation [18], the frequency disturbance will greatly increase the probability of wide range fluctuation of the PSU’s frequency. The traditional PID controller cannot adjust the frequency fluctuation of PSU well. In recent years, FOPID controller has been used in a lot of research studies [19, 20]. Xu et al. designed an adaptive fast fuzzy fractional-order PID controller to improve the dynamic performance of pumped storage units under low water head no-load condition [21]. Chen et al. designed a fractional-order PID controller for hydraulic turbine regulating system [22]. Tang et al. designed a fractional-order controller for automatic voltage regulator [23]. Zhang et al. optimized the FOPID controller of PSURS under multiple working conditions [24]; the common first-order generator model has been applied, and electrical factors such as active power, voltage, current, and excitation regulation system have not been considered deeply. Compared with traditional PID controller, the FOPID controller has two more adjustable parameters, which improve its control performance. However, at the same time, additional two adjustable parameters also make the optimal setting of fractional PID controller parameters more difficult [25].Many excellent intelligent optimization algorithms are widely used in parameter optimization of PSU at present [26]. Malik and Zeng proposed a method for parameter optimization of PID governor of hydraulic turbine generator set based on the algorithm of bacteria-particle swarm optimization [27]. Lai et al. used the artificial flock algorithm to optimize the control of PSU [28]. Zhang et al. proposed a novel chaotic grey wolf optimization algorithm to select the optimal control parameters of the pump turbine governing system [29]. Gravitational search algorithm (GSA) is a new intelligent optimization algorithm proposed by Rashedi et al. [30]. It has strong global searching ability and convergence speed and is more effective in PID parameter setting. Premature and falling into local optimum are two common shortcomings of GSA and other swarm intelligence algorithms, which will reduce the accuracy and efficiency of algorithm optimization [31]. Xu et al. proposed a bacterial-foraging chemotaxis gravitational search algorithm (BCGSA) which can effectively improve these problems [21].Accurate model is essential to optimal regulation and control research. Nonlinear structural model of pump turbine based on the characteristic curve can better describe nonlinear hydraulic and mechanical characteristics. The model of seventh-order synchronous generator fully considers electromagnetic transient and subtransient processes and can better respond to frequency changes [32]. Therefore, this paper puts forward an accurate hydraulic-mechanical-electrical fine coupling system mathematical model of PSURS based on and high-order synchronous generator and full characteristic curve of pump turbine. Moreover, the primary frequency regulation characteristics of the PSURS are tested by numerical simulation. The results show that it is easy to enter the unstable area and cause a violent oscillation of the PSURS under the condition of low water head, small load, and large frequency disturbance. Choosing a suitable optimal control method is still an important and challenging issue for this unstable condition. So, the FOPID controller with BCGSA is adopted during frequent disturbance process of +0.5 Hz and −0.5 Hz under low water head and small load condition.The remainder of the paper is structured as follows. Section2 establishes the hydraulic‐mechanical‐electrical coupling model of pumped storage unit regulation system. Section 3 proposes frequency regulation performance analysis of pumped storage units with different water heads and loads. Section 4 describes the optimization method of control parameters including FOPID and BCGSA. Section 5 discusses the experimental results, and Section 6 summarizes the conclusion of the research.


## 2. Hydraulic‐Mechanical‐Electrical Coupling Model

PSURS is a complex closed-loop control system which includes the material and energy exchange of hydraulic, mechanical, and electrical factors which is composed of a water diversion system [33], pump turbine, synchronous generator, excitation system, and governor, as shown in Figure 1.Figure 1
Structure of PSURS.


### 2.1. Model of Water Diversion System

In the modeling of PSUSRS, the water diversion system is usually reduced to a section of equivalent pipe and modeled by water hammer theory. Equivalent simplified models of water diversion system can be divided into rigid water hammer models, approximate elastic water hammer models, and elastic water hammer models. In this paper, the primary frequency regulation condition is small fluctuation condition, and the hydraulic fluctuation is minor variation. Considering the simulation accuracy and hydraulic calculation efficiency, an approximate elastic water hammer model is adopted, and its transfer function is shown in the following equation [34]:(1)Ghs=−hwTrs1+1/2fTrs+1/8Tr2s2.


### 2.2. Model of Pump Turbine

Pump turbine is the key component of PSURS, and there is no accurate analytic expression to describe the dynamic characteristics of this system until now. The flow and torque characteristics of the pump turbine are complex. It is even impossible to obtain the accurately analytic expression of the pump turbine. Through the model experiment, the comprehensive characteristic curve of the pump turbine and the steady-state characteristics of the pump turbine like the pump turbine flow and torque characteristics can be achieved. The instantaneous parameters of the unit are calculated by the characteristic curve, and the dynamic characteristics are obtained.At present, mathematical modeling of the pump turbine mainly includes linear analytical model, internal characteristic model, and interpolation model based on full characteristic curve. The hydraulic-mechanical-electrical coupling model should be able to describe the dynamic characteristics of the process of converting water energy into rotating mechanical energy, so interpolation model based on full characteristic curve is adopted in this paper [21]. The characteristic curve of the pump turbine takes the true value of the guide vane opening y as parameter to establish the flow characteristic curve Q11∼N11 and torque characteristic curve M11∼N11 of the pump turbine, as shown in the following equation [35]:(2)Q11=Q11y,N11,M11=M11y,N11.Due to the “S” area of the characteristic curve of the pump turbine, the problem of “one input and many outputs” exists in the interpolation calculation. The logarithmic curve projection (LCP) method not only solves the multivalue problem of zero opening degree, small opening degree, and “S” characteristic region of guide vane opening in the description of runner characteristics but also facilitates the solution and has high precision. The LCP is described as follows [36]:(3)n1=N11N11r,(4)v1=Q11Q11r,(5)m1=M11M11r,(6)x=n1ev1,where Q11r, N11r, and M11r are rated unit flow, rated unit rotational speed, and rated unit torque.Characteristic curved surfaces obtained by LCP are shown in Figure2. The flow Q11t and torque M11t calculation process of the pump turbine mathematical model is summarized as the following:Step (1): Initialization. Initialize the guide opening yt and unit rotational speed n11t at t time and unit flow Q11t−Δt at t−Δt time.Step (2): Transformation. Let Q11t0=Q11t−Δt and calculate the abscissa xtn using transform formulas (3)–(6) of LCP.Step (3): Interpolation. For curves x−Q11 and x−M11 of Figure 2, calculate M11tn and Q11tn using yt and xtn as inputs by the interpolation method.Step (4). Calculate absolute value error of Q11tn+1 and Q11tn, while n changes from 1 to Niter. If Q11tn+1−Q11tn<σ, directly output Q11t and M11t; else update Q11tn+1 by Q11tn+1=Q11tn−τ⋅Q11tn−Q11tn−1, where τ is the scale factor and σ is the calculation accuracy coefficient.Step (5): Iterative Calculation. Repeat Step (2) to Step (4) until the stop criteria reached. When n=Niter and the stop criteria are not reached, if Q11tn+1−Q11tn≤10σ, directly output Q11t and M11t; else output values according to the given formulas Q11t=Q11tn+5σ and M11t=M11tn+5σ.Figure 2
Logarithmic curve projection of the characteristic curves.
(a)
(b)


### 2.3. Model of Synchronous Generator

Synchronous generator is the core of the power system. Its dynamic performance is very complex. In the mathematical model of synchronous motors, the zero shaft winding is not considered. The electromagnetic transition process of the five windingsd, q, f, D, and Q (winding flux or current as the state variable) and the rotor mechanical transition process (δ, ω as the state variable) constitute the synchronous motor as a seventh-order model [32, 37]. In this paper, the combined transition process of water and electricity is studied. Considering the electromagnetic transient and subtransient processes, the seventh-order generator model is selected.The voltage equation of each winding of the synchronous motor is shown in the following equation (Park equation) [38]:(7)vd=−raid+ψ•d−ωψq,vd=−raiq+ψ•q−ωψd,vf=rfif+ψ•f,vD=0=rDiD+ψ•D,vQ=0=rQiQ+ψ•Q,where the subscript d represents the d axis winding, the subscript q represents the q axis winding, the subscript f represents the rotor excitation winding, the subscript D represents the d axis equivalent damping winding, and the subscript Q represents the q axis equivalent damping winding.The flux linkage equation of axisd and axis q reflects the magnetic coupling relationship between stator winding and rotor winding. The specific equation is shown in the following equation:(8)ψd=−xdid+xadif+iD,ψq=−xqiq+xaqiQ,ψf=xfif+xadiD−id,ψD=xDiD+xadif−id,ψQ=xQiQ−xaqiq.The rotor motion equation is shown in the following equation [39, 40]:(9)dδdt=ω−1ω0,dωdt=T−mψdiq−ψqid2H.The second-order model of the generator rotor and the relationship between generator inertia and output power and power angle are considered, and the relationship between power angle and rotation speed is established. The second-order model of the rotor and the voltage equation of equation (7) constitute the seventh-order model of the generator.


### 2.4. Model of Excitation System

The excitation system provides the excitation power to the generator, acts as regulating the voltage, keeping the voltage at the end of the generator or the pivot point constant, and can control the reactive power distribution of parallel running generators. It can improve the stability limit of the power system and greatly improve the transient stability of the power system.There are various kinds of excitation systems for synchronous generators. Different excitation models can be chosen according to the excitation mode of the hydropower station. The transient limiter is used in IEEEST1A model, and the intrinsic excitation time constant is very small, so the stability of the excitation system can be ignored. Because this paper focuses on the dynamic stability of generator output frequency under the condition of primary frequency regulation, the stability of the excitation system can be ignored. So, this paper adopts the international standard IEEE ST1A excitation system model [41].


### 2.5. Model of Speed Governor

Microcomputer speed governor is commonly used in the speed control of the pump turbine in modern pumped storage power station. The speed governor is comprised of a microcomputer regulator and electro-hydraulic follow-up device. Currently, the widely used microcomputer regulator adopts the parallel PID controller shown in Figure3, and its transfer function is shown in the following equation:(10)GPIDs=Kp+KdsTds+1+Kisnref−n+bpyref−y.Figure 3
The parallel PID controller.Electro-hydraulic actuator is the actuator of the governor. The servo mechanism and the main valve of actuator can be described by two first-order inertia links in series. The transfer function of the electro-hydraulic actuator is shown in the following equation:(11)Gs=k0TyBs+1×1Tys.The purpose of PSURS optimization is to suppress the instability of various operating conditions and make the system safe and stable. To accurately reflect the strong nonlinearity of the regulation system of PSU under the coupling of hydraulic, mechanical, and electrical factors, the nonlinear model structure of the regulation system was adopted in this paper, as shown in Figure4.Figure 4
Hydraulic-mechanical-electrical coupling model.


## 2.1. Model of Water Diversion System

In the modeling of PSUSRS, the water diversion system is usually reduced to a section of equivalent pipe and modeled by water hammer theory. Equivalent simplified models of water diversion system can be divided into rigid water hammer models, approximate elastic water hammer models, and elastic water hammer models. In this paper, the primary frequency regulation condition is small fluctuation condition, and the hydraulic fluctuation is minor variation. Considering the simulation accuracy and hydraulic calculation efficiency, an approximate elastic water hammer model is adopted, and its transfer function is shown in the following equation [34]:(1)Ghs=−hwTrs1+1/2fTrs+1/8Tr2s2.


## 2.2. Model of Pump Turbine

Pump turbine is the key component of PSURS, and there is no accurate analytic expression to describe the dynamic characteristics of this system until now. The flow and torque characteristics of the pump turbine are complex. It is even impossible to obtain the accurately analytic expression of the pump turbine. Through the model experiment, the comprehensive characteristic curve of the pump turbine and the steady-state characteristics of the pump turbine like the pump turbine flow and torque characteristics can be achieved. The instantaneous parameters of the unit are calculated by the characteristic curve, and the dynamic characteristics are obtained.At present, mathematical modeling of the pump turbine mainly includes linear analytical model, internal characteristic model, and interpolation model based on full characteristic curve. The hydraulic-mechanical-electrical coupling model should be able to describe the dynamic characteristics of the process of converting water energy into rotating mechanical energy, so interpolation model based on full characteristic curve is adopted in this paper [21]. The characteristic curve of the pump turbine takes the true value of the guide vane opening y as parameter to establish the flow characteristic curve Q11∼N11 and torque characteristic curve M11∼N11 of the pump turbine, as shown in the following equation [35]:(2)Q11=Q11y,N11,M11=M11y,N11.Due to the “S” area of the characteristic curve of the pump turbine, the problem of “one input and many outputs” exists in the interpolation calculation. The logarithmic curve projection (LCP) method not only solves the multivalue problem of zero opening degree, small opening degree, and “S” characteristic region of guide vane opening in the description of runner characteristics but also facilitates the solution and has high precision. The LCP is described as follows [36]:(3)n1=N11N11r,(4)v1=Q11Q11r,(5)m1=M11M11r,(6)x=n1ev1,where Q11r, N11r, and M11r are rated unit flow, rated unit rotational speed, and rated unit torque.Characteristic curved surfaces obtained by LCP are shown in Figure2. The flow Q11t and torque M11t calculation process of the pump turbine mathematical model is summarized as the following:Step (1): Initialization. Initialize the guide opening yt and unit rotational speed n11t at t time and unit flow Q11t−Δt at t−Δt time.Step (2): Transformation. Let Q11t0=Q11t−Δt and calculate the abscissa xtn using transform formulas (3)–(6) of LCP.Step (3): Interpolation. For curves x−Q11 and x−M11 of Figure 2, calculate M11tn and Q11tn using yt and xtn as inputs by the interpolation method.Step (4). Calculate absolute value error of Q11tn+1 and Q11tn, while n changes from 1 to Niter. If Q11tn+1−Q11tn<σ, directly output Q11t and M11t; else update Q11tn+1 by Q11tn+1=Q11tn−τ⋅Q11tn−Q11tn−1, where τ is the scale factor and σ is the calculation accuracy coefficient.Step (5): Iterative Calculation. Repeat Step (2) to Step (4) until the stop criteria reached. When n=Niter and the stop criteria are not reached, if Q11tn+1−Q11tn≤10σ, directly output Q11t and M11t; else output values according to the given formulas Q11t=Q11tn+5σ and M11t=M11tn+5σ.Figure 2
Logarithmic curve projection of the characteristic curves.
(a)
(b)


## 2.3. Model of Synchronous Generator

Synchronous generator is the core of the power system. Its dynamic performance is very complex. In the mathematical model of synchronous motors, the zero shaft winding is not considered. The electromagnetic transition process of the five windingsd, q, f, D, and Q (winding flux or current as the state variable) and the rotor mechanical transition process (δ, ω as the state variable) constitute the synchronous motor as a seventh-order model [32, 37]. In this paper, the combined transition process of water and electricity is studied. Considering the electromagnetic transient and subtransient processes, the seventh-order generator model is selected.The voltage equation of each winding of the synchronous motor is shown in the following equation (Park equation) [38]:(7)vd=−raid+ψ•d−ωψq,vd=−raiq+ψ•q−ωψd,vf=rfif+ψ•f,vD=0=rDiD+ψ•D,vQ=0=rQiQ+ψ•Q,where the subscript d represents the d axis winding, the subscript q represents the q axis winding, the subscript f represents the rotor excitation winding, the subscript D represents the d axis equivalent damping winding, and the subscript Q represents the q axis equivalent damping winding.The flux linkage equation of axisd and axis q reflects the magnetic coupling relationship between stator winding and rotor winding. The specific equation is shown in the following equation:(8)ψd=−xdid+xadif+iD,ψq=−xqiq+xaqiQ,ψf=xfif+xadiD−id,ψD=xDiD+xadif−id,ψQ=xQiQ−xaqiq.The rotor motion equation is shown in the following equation [39, 40]:(9)dδdt=ω−1ω0,dωdt=T−mψdiq−ψqid2H.The second-order model of the generator rotor and the relationship between generator inertia and output power and power angle are considered, and the relationship between power angle and rotation speed is established. The second-order model of the rotor and the voltage equation of equation (7) constitute the seventh-order model of the generator.


## 2.4. Model of Excitation System

The excitation system provides the excitation power to the generator, acts as regulating the voltage, keeping the voltage at the end of the generator or the pivot point constant, and can control the reactive power distribution of parallel running generators. It can improve the stability limit of the power system and greatly improve the transient stability of the power system.There are various kinds of excitation systems for synchronous generators. Different excitation models can be chosen according to the excitation mode of the hydropower station. The transient limiter is used in IEEEST1A model, and the intrinsic excitation time constant is very small, so the stability of the excitation system can be ignored. Because this paper focuses on the dynamic stability of generator output frequency under the condition of primary frequency regulation, the stability of the excitation system can be ignored. So, this paper adopts the international standard IEEE ST1A excitation system model [41].


## 2.5. Model of Speed Governor

Microcomputer speed governor is commonly used in the speed control of the pump turbine in modern pumped storage power station. The speed governor is comprised of a microcomputer regulator and electro-hydraulic follow-up device. Currently, the widely used microcomputer regulator adopts the parallel PID controller shown in Figure3, and its transfer function is shown in the following equation:(10)GPIDs=Kp+KdsTds+1+Kisnref−n+bpyref−y.Figure 3
The parallel PID controller.Electro-hydraulic actuator is the actuator of the governor. The servo mechanism and the main valve of actuator can be described by two first-order inertia links in series. The transfer function of the electro-hydraulic actuator is shown in the following equation:(11)Gs=k0TyBs+1×1Tys.The purpose of PSURS optimization is to suppress the instability of various operating conditions and make the system safe and stable. To accurately reflect the strong nonlinearity of the regulation system of PSU under the coupling of hydraulic, mechanical, and electrical factors, the nonlinear model structure of the regulation system was adopted in this paper, as shown in Figure4.Figure 4
Hydraulic-mechanical-electrical coupling model.


## 3. Primary Frequency Regulation Performance Analysis

The primary frequency regulation function of PSU is very important to maintain the stability of power grid frequency, which can provide power support rapidly when the power grid suddenly changes in large load, improve the reliability of power system, optimize system dispatching, and stabilize the power grid frequency [42]. In order to study the dynamic regulation process of PSURS, the simulation experiment is performed in this section. The refined mathematical model of the regulation system of the pumped storage unit is shown in Section 2. The model of the pump turbine uses the data of the actual unit of a pumped storage power station in China. The water head and model parameters of the power station are shown in Table 1.Table 1
Water head and parameters of simulation model of a power station.
Water head (m)
Transfer parameters
H
max
H
r
H
min
h
w
T
yB
T
y
T
r
k
0
f
217.0
195.0
187.0
1.7
0.02
0.2
1.0
1.0
0.75The control parametersKp = 3, Ki = 2, and Kd = 4 are used to analyze the dynamic characteristics of the system when the PSURS receives frequency disturbances of 0.2 Hz and 0.5 Hz at different water heads with different loads. In this paper, the maximum opening of the pump turbine ymax is 43.1°, and the setting opening under different guide vane opening and different loads during the experimental analysis is shown in Table 2. The time domain performance analysis indicator is integral of time multiplied absolute error (ITAE). ITAE values of frequency f and guide vane opening y under various working conditions are shown in Tables 3–6. For convenience of comparison, data in the table are normalized. By these Tables, it is shown that under the same water head, the lower the unit with the load, the greater the ITAE_f ( the ITAE value of the frequency ) and the ITAE_y (the ITAE value of the opening),whether it has big frequency disturbance (Δf=−0.5Hz) or small frequency disturbance (Δf=−0.2Hz). If the ITAE value is relatively large, then the system is easier to instability and dynamic characteristics is worse, especially as Δf=−0.5Hz. The ITAE_f and ITAE_y with 25% load are much higher than the value with greater load. At the same time, under the same load, when Δf=−0.2Hz, ITAE value also increases slowly with the increase of water head. But when Δf=−0.5Hz, the ITAE values of low water head (Hw = 190 m) with 25% load and 50% load are larger than other water head, and the instability degree is higher at 25% load.Table 2
The set value of guide vane opening at different water heads and loads.
y
max = 43.1°
H
w
m
190
195
200
210
Opening of 25% load (%)
31
30
27
24
Opening of 50% load (%)
45
42
40
38
Opening of 75% load (%)
64
60
56
52
Opening of 100% load (%)
100
95
90
78Table 3
The ITAE value of the output frequency ITAE_f when Δf=−0.5Hz.
Δ
f
=
−
0.5
Hz
H
w
m
190
195
200
210
25% load
0.9065
0.5899
0.6225
0.6025
50% load
0.1527
0.1055
0.1119
0.1217
75% load
0.0649
0.0877
0.0960
0.1031
100% load
0.0032
0.0129
0.0359
0.0429Table 4
The ITAE value of the output frequency ITAE_f when Δf=−0.2Hz.
Δ
f
=
−
0.2
Hz
H
w
m
190
195
200
210
25% load
0.5666
0.7216
0.8265
0.8966
50% load
0.4498
0.4569
0.4652
0.5117
75% load
0.0751
0.3650
0.4291
0.4525
100% load
0.0209
0.0185
0.0153
0.3211Table 5
The ITAE value of the output opening ITAE_y when Δf=−0.5Hz.
Δ
f
=
−
0.5
Hz
H
w
m
190
195
200
210
25% load
0.9041
0.5703
0.6221
0.6276
50% load
0.1764
0.0098
0.1135
0.0099
75% load
0.0084
0.0086
0.0089
0.0092
100% load
0.0061
0.0065
0.0069
0.0075Table 6
The ITAE value of the output opening ITAE_y when Δf=−0.2Hz.
Δ
f
=
−
0.2
Hz
H
w
m
190
195
200
210
25% load
0.6971
0.7117
0.7254
0.7436
50% load
0.5992
0.6141
0.6302
0.6453
75% load
0.4548
0.4909
0.5205
0.5497
100% load
0.2023
0.2362
0.2710
0.3643After analysis of water head and frequency, the comprehensive ITAE index was set up as shown in equation (12), and its specific value and change trend are shown in Tables 7 and 8 and Figure 5. According to Figures 5(a) and 5(b), when the unit is carrying 25% load, the PSU is more likely to enter into an unstable state after frequency disturbance, and its performance is worse, especially under low water head (Hw = 190 m) and large Δf=−0.5Hz.Table 7
The comprehensive ITAE Index whenΔf=−0.5Hz.
Δ
f
=
−
0.5
Hz
H
w
m
190
195
200
210
25% load
1.8106
1.1603
1.2446
1.2301
50% load
0.3290
0.1153
0.2254
0.1316
75% load
0.0733
0.0963
0.1049
0.1123
100% load
0.0093
0.0194
0.0428
0.0504Table 8
The comprehensive ITAE Index whenΔf=−0.2Hz.
Δ
f
=
−
0.2
Hz
H
w
m
190
195
200
210
25% load
1.2637
1.4333
1.5519
1.6401
50% load
1.0489
1.0710
1.0954
1.1631
75% load
0.5299
0.8559
0.9495
1.0022
100% load
0.2232
0.2547
0.2863
0.6854Figure 5
The comprehensive ITAE index of different water heads with different loads. (a)Δf=−0.5Hz. (b) Δf=−0.2Hz.
(a)
(b)Figure6 shows the comparison between low water head Hw = 190 m and high water head (Hw = 210 m) of guide vane opening y, water hammer pressure h, and dynamic PE under different loads. Comparative analysis shows that the lower the load, the worse the dynamic characteristics of the system. At 25% load, the PSURS experienced violent oscillation, and when Hw = 190 m, the amplitude of oscillation was larger, especially violent.Figure 6
The comparison diagram ofy, h, and pe. (a) 25% load. (b) 50% load. (c) 75% load. (d) 100% load.
(a)
(b)
(c)
(d)By analyzing the simulation experiment of frequency disturbance of the regulating system with different water heads and different load conditions, it can be known that PSURS is easy to enter the unstable area and cause violent oscillation of the system under the condition of low water head, small load, and large frequency disturbance. At this time, the control capability of the PSURS is required to increase, and appropriate control strategies and parameters are particularly important to maintain the stability of the frequency of the units and the power grid.


## 4. Control Parameter Optimization Method

### 4.1. Fractional-Order PID Controller

The fractional-order PID controller is denoted as PIλDμ [22]. It is a generalization of classical integer-order controller, which is more concise and accurate for complex practical systems [43]. Due to the introduction of differential order λ and integral order μ, the whole controller has two more adjustable parameters, so the setting range of controller parameters becomes larger, the controller can control the controlled object more flexibly, and better control results can be expected. Its transfer function is shown in equation (12) [44]. In order to achieve a trade-off between the complexity and accuracy of the achievable FOPID controller, 5-th order Oustaloup’s recursive approximation is done for the integral-differential operators within a chosen frequency band of ω∈10−2,102rad/s.(12)Gcs=Kp+Kisλ+Kdsμ,where λ and μ are, respectively, coefficients of integral and differential operators, which are any real numbers, and λ, μ > 0.The fractional-order PID is applied to the parallel PID controller model mentioned above, and the improved parallel FOPID governor model is obtained, as shown in Figure7 [45].Figure 7
The parallel FOPID controller.


### 4.2. Proposed Optimization Objective Function

Various time domain integral performance indexes like integral of time multiplied absolute error (ITAE), integral of time multiplied squared error (ITSE), integral of squared time multiplied error whole squared (ISTES), and integral of squared time multiplied squared error (ISTSE) are considered in the problem similarly, and every index has certain advantages in control system design. For ITAE criterion, time multiplied absolute error of the control system is judged to minimize, and the maximum percentage of overshoot is also minimized since the absolute error is included in the ITAE. The ITSE criterion penalizes the error more than the ITAE, and due to the time multiplication term, the oscillation damps out faster. Other integral performance indexes like ISTES and ISTSE both have higher power of time and error terms. When the pumped storage unit is running under turbine conditions, it mainly is dependent on the parameters of the adjustment device to improve the dynamic quality and stability of the regulation system. ITAE can well express the adjustment speed and depth factors concerned in the optimization of primary frequency regulation parameters. Therefore, in order to improve the dynamic performance of the system after the disturbance is generated, the time domain performance analysis indicators of ITAE is selected in this paper. As shown in equation (13), the weighted structure of the three parameters of water hammer pressure, frequency, and guide vane opening output by the PSURS is adopted. In particular, ITAE values of the three parameters are normalized and weight is redistributed due to their different ranges.(13)J=w1×ITAE_h+w2×ITAE_y+w3×ITAE_f,where ITAE_h, ITAE_y, and ITAE_f are the ITAE values of water hammer pressure, guide vane opening, and frequency, respectively, and w1, w2, and w3 are the weights of the control objective function.


### 4.3. BCGSA

The GSA is a kind of global optimization algorithm which makes use of the universal gravitation interaction between particles to search intelligently [46]. Premature and falling into local optimum are two common shortcomings of GSA and other swarm intelligence algorithms, which will reduce the accuracy and efficiency of algorithm optimization. The BCGSA introduces Pbest-Gbest directed particle update strategy of PSO algorithm, adaptive elastic ball valve boundary processing method, and chemotactic operation in bacterium foraging algorithms (BFA). These three improvement strategies can effectively avoid the algorithm premature and fall into the local optimum and they also can improve the global optimization ability. The advantage of BCGSA for multiparameter optimization has been tested [21].The computational process of BCGSA algorithm can be summarized as shown in Figure8. When updating the velocity and position of particles, PSO algorithm’s Pbest-Gbest oriented particle updating strategy is applied, and the mathematical expressions are shown in equations (14) and (15):(14)vidt+1=r1⋅vidt+aidt+c1⋅r2⋅Pibestdt−xidt+c2⋅r3⋅Gibestdt−xidt,(15)xidt+1=xidt+vidt,where r1, r2, and r3 are random numbers between [0, 1]; c1 and c2 are learning factors whose value range is [0, 2]; vidt is the i particle’s velocity in the d dimension; xidt represents the position of the i particle in the d dimension; Pibest(t) is the optimal value of the current algebra of the i particle; Gbest(t) is the global optimal value of the group; Gworst(t) is the global optimal value of the group; Xbest is the global optimal particle of the current algebra; and Xworst is the global worst particle of the current algebra.Figure 8
BCGSA optimization steps.When the velocity and position of each particle are updated, for the particles that cross the threshold boundary, the adaptive pinball method is adopted to process, simulate the bounce characteristics of the elastic ball, and introduce the adaptive elastic coefficientζ, as shown in the following equation:(16)ifxidt>Ubd,up=xidt−Ubdandxidt=Ubd−ζ⋅up,ifxidt<Lbd,down=Lbd−xidtandxidt=Lbd+ζ⋅down,where Ub(d) and Lb(d) are the upper and lower boundaries of the threshold value of the dimension. For the few particles still clustered at the boundary, their positions will be changed according to the following equation:(17)Ifxidt>Ubdorxidt<Lbd,xidt=rand⋅Ubd−Lbd+Lbd.The fitness function value of each particle is calculated, as shown in equations (18) and (19). Chemotaxis local search operation is carried out for global optimal particles and global worst particles of the current algebra according to a certain number of chemotactic steps and moving steps.(18)xi,j+1=xi,j+Ciϕj,(19)ϕj=ΔiΔTiΔi,where x(i, j) is the position of the ith particle at the jth chemotactic, C(i) is the moving step length after the random direction, ϕ(j) is the flipped random direction, and the vector Δ(i) composed of random number elements between [−1, 1] is randomly generated. Global optimal particle and global worst particle can be updated only when the new particle is better than the original value after the completion of the chemotactic step.


## 4.1. Fractional-Order PID Controller

The fractional-order PID controller is denoted as PIλDμ [22]. It is a generalization of classical integer-order controller, which is more concise and accurate for complex practical systems [43]. Due to the introduction of differential order λ and integral order μ, the whole controller has two more adjustable parameters, so the setting range of controller parameters becomes larger, the controller can control the controlled object more flexibly, and better control results can be expected. Its transfer function is shown in equation (12) [44]. In order to achieve a trade-off between the complexity and accuracy of the achievable FOPID controller, 5-th order Oustaloup’s recursive approximation is done for the integral-differential operators within a chosen frequency band of ω∈10−2,102rad/s.(12)Gcs=Kp+Kisλ+Kdsμ,where λ and μ are, respectively, coefficients of integral and differential operators, which are any real numbers, and λ, μ > 0.The fractional-order PID is applied to the parallel PID controller model mentioned above, and the improved parallel FOPID governor model is obtained, as shown in Figure7 [45].Figure 7
The parallel FOPID controller.


## 4.2. Proposed Optimization Objective Function

Various time domain integral performance indexes like integral of time multiplied absolute error (ITAE), integral of time multiplied squared error (ITSE), integral of squared time multiplied error whole squared (ISTES), and integral of squared time multiplied squared error (ISTSE) are considered in the problem similarly, and every index has certain advantages in control system design. For ITAE criterion, time multiplied absolute error of the control system is judged to minimize, and the maximum percentage of overshoot is also minimized since the absolute error is included in the ITAE. The ITSE criterion penalizes the error more than the ITAE, and due to the time multiplication term, the oscillation damps out faster. Other integral performance indexes like ISTES and ISTSE both have higher power of time and error terms. When the pumped storage unit is running under turbine conditions, it mainly is dependent on the parameters of the adjustment device to improve the dynamic quality and stability of the regulation system. ITAE can well express the adjustment speed and depth factors concerned in the optimization of primary frequency regulation parameters. Therefore, in order to improve the dynamic performance of the system after the disturbance is generated, the time domain performance analysis indicators of ITAE is selected in this paper. As shown in equation (13), the weighted structure of the three parameters of water hammer pressure, frequency, and guide vane opening output by the PSURS is adopted. In particular, ITAE values of the three parameters are normalized and weight is redistributed due to their different ranges.(13)J=w1×ITAE_h+w2×ITAE_y+w3×ITAE_f,where ITAE_h, ITAE_y, and ITAE_f are the ITAE values of water hammer pressure, guide vane opening, and frequency, respectively, and w1, w2, and w3 are the weights of the control objective function.


## 4.3. BCGSA

The GSA is a kind of global optimization algorithm which makes use of the universal gravitation interaction between particles to search intelligently [46]. Premature and falling into local optimum are two common shortcomings of GSA and other swarm intelligence algorithms, which will reduce the accuracy and efficiency of algorithm optimization. The BCGSA introduces Pbest-Gbest directed particle update strategy of PSO algorithm, adaptive elastic ball valve boundary processing method, and chemotactic operation in bacterium foraging algorithms (BFA). These three improvement strategies can effectively avoid the algorithm premature and fall into the local optimum and they also can improve the global optimization ability. The advantage of BCGSA for multiparameter optimization has been tested [21].The computational process of BCGSA algorithm can be summarized as shown in Figure8. When updating the velocity and position of particles, PSO algorithm’s Pbest-Gbest oriented particle updating strategy is applied, and the mathematical expressions are shown in equations (14) and (15):(14)vidt+1=r1⋅vidt+aidt+c1⋅r2⋅Pibestdt−xidt+c2⋅r3⋅Gibestdt−xidt,(15)xidt+1=xidt+vidt,where r1, r2, and r3 are random numbers between [0, 1]; c1 and c2 are learning factors whose value range is [0, 2]; vidt is the i particle’s velocity in the d dimension; xidt represents the position of the i particle in the d dimension; Pibest(t) is the optimal value of the current algebra of the i particle; Gbest(t) is the global optimal value of the group; Gworst(t) is the global optimal value of the group; Xbest is the global optimal particle of the current algebra; and Xworst is the global worst particle of the current algebra.Figure 8
BCGSA optimization steps.When the velocity and position of each particle are updated, for the particles that cross the threshold boundary, the adaptive pinball method is adopted to process, simulate the bounce characteristics of the elastic ball, and introduce the adaptive elastic coefficientζ, as shown in the following equation:(16)ifxidt>Ubd,up=xidt−Ubdandxidt=Ubd−ζ⋅up,ifxidt<Lbd,down=Lbd−xidtandxidt=Lbd+ζ⋅down,where Ub(d) and Lb(d) are the upper and lower boundaries of the threshold value of the dimension. For the few particles still clustered at the boundary, their positions will be changed according to the following equation:(17)Ifxidt>Ubdorxidt<Lbd,xidt=rand⋅Ubd−Lbd+Lbd.The fitness function value of each particle is calculated, as shown in equations (18) and (19). Chemotaxis local search operation is carried out for global optimal particles and global worst particles of the current algebra according to a certain number of chemotactic steps and moving steps.(18)xi,j+1=xi,j+Ciϕj,(19)ϕj=ΔiΔTiΔi,where x(i, j) is the position of the ith particle at the jth chemotactic, C(i) is the moving step length after the random direction, ϕ(j) is the flipped random direction, and the vector Δ(i) composed of random number elements between [−1, 1] is randomly generated. Global optimal particle and global worst particle can be updated only when the new particle is better than the original value after the completion of the chemotactic step.


## 5. Case Study and Analysis

In this section, a comprehensive simulation model of PSURS, which can fully reflect the hydraulic, mechanical, electrical, and other state characteristics of the unit described in Section2, is built in MATLAB. This paper adopts the actual parameters of a pumped storage power station with a unit capacity of 300 MW in Hubei province, China, which has been applied in practical studies for many times. The traditional PID controller and FOPID controller are, respectively, used in the speed governor part to simulate the control process of the frequency disturbance process under the condition of low water head and small load, and the objective function mentioned in BCGSA and Section 4.2 is used to optimize the control parameters. Furthermore, through comparison with the traditional PID control method, the superiority of the control performance of the selected scheme is verified.The water head and PSU model parameters of the pumped storage power station are shown in Table1. It is easy to enter the “S” characteristic area when the unit is running with low water head near the rated water head. Especially in low water head and low water load conditions, the PSU is particularly easy to enter unstable areas when subjected to large frequency disturbances. In this paper, 190 m water head and 25% rated load are selected for the simulation experiment. The population size of BCGSA is 30, and the maximum number of iterations of the algorithm is 100, c1 = c2 = 2. The control parameters of PSURS are essentially the same as those of the conventional governor. The time range of model simulation is 55 s. PID controller parameter setting range according to relevant regulations and experience is KP∈0.5,4, Ki∈0.05,4, and Kd∈0.05,4. The adjustment order of the integral order and differential order of FOPID is taken as λ,μ∈0.1,1.1.


### 5.1. Frequency Disturbance of −0.5 Hz under Low Water Head and Small Load Condition

In the control optimization of nonlinear PSURS, when BCGSA is used to set the control parameters of complex time-varying regulation system, the convergence curves of the objective function of the particle swarm optimization (PSO) algorithm, gravitational search algorithm(GSA), and BCGSA are shown in Figure9. Under the same parameter setting, the overall value of the objective function calculated by BCGSA for the unit optimization controlled by FOPID is lower, and the convergence value is also smaller. According to the convergence curves of PID-PSO, PID- GSA, PID-BCGSA, FOPID-PSO, FOPID- GSA, and FOPID-BCGSA, it is clear that the performance of BCGSA is better than PSO and the standard GSA. It also verifies that the FOPID control effect is generally better than PID control when the unit suffers 0.5 HZ frequency disturbance under the condition of low water head and 25% load.Figure 9
Convergence curves of objective function withΔf=−0.5Hz.With reference to the above parameter settings, based on the refined model of PSU, BCGSA is used to optimize the control parameters, and the simulation experiment of frequency disturbance under low water head and small load conditions is carried out. The optimized FOPID and PID control parameters are shown in Table9. The optimal convergence value of the objective function and the ITAE values of the frequency f, water hammer pressure h, power pe, and guide vane opening y are shown in Table 10. Regulating time ts and the time tf at which the output enters violent oscillation of each output are shown in Table 11. tf is the time when the first wave peak or trough occurs when the output oscillates violently. It can be found from the experiment that there is no significant difference in the length of time when the unit is in a violent oscillation. Therefore, the time when the output enters the first wave peak or trough in the continuous oscillation region is taken as the indicator to consider the regulation effect of PSU. The larger the tf value is, the slower the unit’s response speed to frequency disturbance is, and the smaller the tf is, the faster the unit’s response speed is. ts is the regulating time in the regional energy index. The transition process curve for each output is shown in Figure 10.Table 9
The optimized controller parameters withΔf=−0.5Hz.
Case
Controller parameters
K
p
K
i
K
d
λ
μ
PID
3.8507
1.7732
1.2081
—
—
FOPID
2.8946
4.0834
1.7732
0.8049
0.8715Table 10
Performance indicators withΔf=−0.5Hz.
Index
J
ITAE
_
f
ITAE
_
h
ITAE
_
pe
ITAE
_
y
PID
7.9122
4.1133 × 103
2.9083 × 104
7.0998 × 105
2.2788 × 104
FOPID
6.9413
2.9674 × 103
1.7478 × 104
1.0996 × 105
1.0177 × 104Table 11
Comparison of performance indexes withΔf=−0.5Hz.
Indexes
PID
FOPID
f
pe
h
y
f
pe
h
y
t
f
s
15.05
15.12
15.06
15.07
10.42
10.71
10.40
10.53
t
s
s
25.20
43.98
20.03
40.71
17.63
35.56
15.11
30.03Figure 10
Dynamic response of each output withΔf=−0.5Hz. (a) Frequency. (b) Water hammer pressure. (c) Power. (d) Guide vane opening.
(a)
(b)
(c)
(d)It can be seen from Table10 that the ITAE values of the output parameters of the unit under the optimal control of FOPID are relatively small, which effectively improves the regulation quality when the unit suffers from 0.5 Hz frequency disturbance under the condition of low water head and 25% load. The small objective function of FOPID also means that the unit gets good dynamic characteristics. As can be seen from Figure 10 and Table 11, after the disturbance occurs, PSU will still enter the oscillation region, but the frequency of oscillation of the unit under the optimal control of FOPID is reduced, and the tf value of each parameter is far smaller than that of the unit under PID control, making the oscillation time shorter. The ts value of the frequency under the control of FOPID is 17.63 s, far less than the 25.20 s under the control of PID, and the frequency is easier to arrive at the stable state. Similarly, the ts value of the power, water head, and opening degree under the control of FOPID is all less than that of the unit under the control of PID, which indicates that PSU can reach the stable state faster under the control of FOPID. Although the overshoot of FOPID control is slightly larger than that of PID, the dynamic indexes under both controllers meet the requirements of national standards.


### 5.2. Frequency Disturbance of +0.5 Hz under Low Water Head and Small Load Condition

When frequency disturbance of +0.5 Hz under low water head and small load occurs, different algorithms are used to optimize the controller parameters and the iterative curve of the objective function is shown in Figure11. On the whole, lower objective function value fully demonstrates the excellent FOPID control with BCGSA under this condition.Figure 11
Convergence curves of objective function withΔf=+0.5Hz.At the same time, the simulation experiment of frequency disturbance of +0.5 Hz under low water head and 25% load is carried out. The optimized control parameters are shown in Table12, the optimal convergence value of fitness function and the ITAE value of each output are shown in Table 13, the regulating time ts and the time tf are shown in Table 14, and the transition process curve of each output is shown in Figure 12.Table 12
The optimized controller parameters withΔf=+0.5Hz.
Case
Controller parameters
K
p
K
i
K
d
λ
μ
PID
3.9992
2.5840
2.4346
—
—
FOPID
2.8370
3.5545
1.5089
1.0007
0.9924Table 13
Performance indicators withΔf=+0.5Hz.
Index
J
ITAE
_
f
ITAE
_
h
ITAE
_
pe
ITAE
_
y
PID
9.1106
81.7019
4.6269 × 103
1.8867 × 104
7.9778 × 103
FOPID
6.7108
76.1314
1.5369 × 103
0.8624 × 104
1.3310 × 103Table 14
Comparison of performance indexes withΔf=+0.5Hz
Indexes
PID
FOPID
f
pe
h
y
f
pe
h
y
t
f
s
11.22
11.10
11.26
11.06
9.942
9.838
10.00
9.950
t
s
s
15.50
32.87
27.56
42.19
14.12
20.02
22.33
15.56Figure 12
Dynamic response of each output withΔf=+0.5Hz. (a) Frequency. (b) Water hammer pressure. (c) Power. (d) Guide vane opening.
(a)
(b)
(c)
(d)Similarly, it can be obtained from Table13 that the ITAE value under FOPID optimization control is smaller, indicating that the dynamic performance of the unit is more excellent at this time. As can be seen from Figure 12 and Table 14, the opening of the pump turbine and the output power of the synchronous generator are particularly rapid in response speed after frequency disturbance. Under the control of FOPID, the ts value of power is 20.02 seconds, far less than 32.87 s, and the ts value of the opening is less than half of the PID control. It takes less time for each output to reach the stable state, which indicates that FOPID control has better dynamic characteristics and regulation quality.


### 5.3. Robustness Analysis of the Obtained Method

The robustness of PSURS is very important for the stable operation of the unit. The PSURS has the function of water diversion system and water flow inertia constantTw (the value of Tw is hw × Tr) which is an important factor affecting the dynamic performance of PSURS. Increasing the value of Tw will aggravate the fluctuation of water hammer pressure, which can be used to verify the robustness of the proposed method. Therefore, in the test and experiment of the obtained frequency regulation optimization control method, we adjusted the time constant of water flow inertia constant Tw to be 2.1. The simulation experiment of frequency disturbance of −0.5 Hz under low water head and 25% load is carried out. The dynamic response curve of each output is shown in Figure 13. When Tw increases, the controller and control parameters in Table 8 are adopted. Compared with Tw of 1.7, the dynamic response process of frequency, water head, power, and opening of the PSURS increases in a small range of amplitude in the oscillation region, but the overall stability time and steady-state error do not change. The optimal primary frequency control method adopted in this paper is robust to the change of hydraulic parameters.Figure 13
Dynamic response of each output at differentTw with Δf=−0.5Hz. (a) Frequency. (b) Water hammer pressure. (c) Power. (d) Guide vane opening.
(a)
(b)
(c)
(d)


## 5.1. Frequency Disturbance of −0.5 Hz under Low Water Head and Small Load Condition

In the control optimization of nonlinear PSURS, when BCGSA is used to set the control parameters of complex time-varying regulation system, the convergence curves of the objective function of the particle swarm optimization (PSO) algorithm, gravitational search algorithm(GSA), and BCGSA are shown in Figure9. Under the same parameter setting, the overall value of the objective function calculated by BCGSA for the unit optimization controlled by FOPID is lower, and the convergence value is also smaller. According to the convergence curves of PID-PSO, PID- GSA, PID-BCGSA, FOPID-PSO, FOPID- GSA, and FOPID-BCGSA, it is clear that the performance of BCGSA is better than PSO and the standard GSA. It also verifies that the FOPID control effect is generally better than PID control when the unit suffers 0.5 HZ frequency disturbance under the condition of low water head and 25% load.Figure 9
Convergence curves of objective function withΔf=−0.5Hz.With reference to the above parameter settings, based on the refined model of PSU, BCGSA is used to optimize the control parameters, and the simulation experiment of frequency disturbance under low water head and small load conditions is carried out. The optimized FOPID and PID control parameters are shown in Table9. The optimal convergence value of the objective function and the ITAE values of the frequency f, water hammer pressure h, power pe, and guide vane opening y are shown in Table 10. Regulating time ts and the time tf at which the output enters violent oscillation of each output are shown in Table 11. tf is the time when the first wave peak or trough occurs when the output oscillates violently. It can be found from the experiment that there is no significant difference in the length of time when the unit is in a violent oscillation. Therefore, the time when the output enters the first wave peak or trough in the continuous oscillation region is taken as the indicator to consider the regulation effect of PSU. The larger the tf value is, the slower the unit’s response speed to frequency disturbance is, and the smaller the tf is, the faster the unit’s response speed is. ts is the regulating time in the regional energy index. The transition process curve for each output is shown in Figure 10.Table 9
The optimized controller parameters withΔf=−0.5Hz.
Case
Controller parameters
K
p
K
i
K
d
λ
μ
PID
3.8507
1.7732
1.2081
—
—
FOPID
2.8946
4.0834
1.7732
0.8049
0.8715Table 10
Performance indicators withΔf=−0.5Hz.
Index
J
ITAE
_
f
ITAE
_
h
ITAE
_
pe
ITAE
_
y
PID
7.9122
4.1133 × 103
2.9083 × 104
7.0998 × 105
2.2788 × 104
FOPID
6.9413
2.9674 × 103
1.7478 × 104
1.0996 × 105
1.0177 × 104Table 11
Comparison of performance indexes withΔf=−0.5Hz.
Indexes
PID
FOPID
f
pe
h
y
f
pe
h
y
t
f
s
15.05
15.12
15.06
15.07
10.42
10.71
10.40
10.53
t
s
s
25.20
43.98
20.03
40.71
17.63
35.56
15.11
30.03Figure 10
Dynamic response of each output withΔf=−0.5Hz. (a) Frequency. (b) Water hammer pressure. (c) Power. (d) Guide vane opening.
(a)
(b)
(c)
(d)It can be seen from Table10 that the ITAE values of the output parameters of the unit under the optimal control of FOPID are relatively small, which effectively improves the regulation quality when the unit suffers from 0.5 Hz frequency disturbance under the condition of low water head and 25% load. The small objective function of FOPID also means that the unit gets good dynamic characteristics. As can be seen from Figure 10 and Table 11, after the disturbance occurs, PSU will still enter the oscillation region, but the frequency of oscillation of the unit under the optimal control of FOPID is reduced, and the tf value of each parameter is far smaller than that of the unit under PID control, making the oscillation time shorter. The ts value of the frequency under the control of FOPID is 17.63 s, far less than the 25.20 s under the control of PID, and the frequency is easier to arrive at the stable state. Similarly, the ts value of the power, water head, and opening degree under the control of FOPID is all less than that of the unit under the control of PID, which indicates that PSU can reach the stable state faster under the control of FOPID. Although the overshoot of FOPID control is slightly larger than that of PID, the dynamic indexes under both controllers meet the requirements of national standards.


## 5.2. Frequency Disturbance of +0.5 Hz under Low Water Head and Small Load Condition

When frequency disturbance of +0.5 Hz under low water head and small load occurs, different algorithms are used to optimize the controller parameters and the iterative curve of the objective function is shown in Figure11. On the whole, lower objective function value fully demonstrates the excellent FOPID control with BCGSA under this condition.Figure 11
Convergence curves of objective function withΔf=+0.5Hz.At the same time, the simulation experiment of frequency disturbance of +0.5 Hz under low water head and 25% load is carried out. The optimized control parameters are shown in Table12, the optimal convergence value of fitness function and the ITAE value of each output are shown in Table 13, the regulating time ts and the time tf are shown in Table 14, and the transition process curve of each output is shown in Figure 12.Table 12
The optimized controller parameters withΔf=+0.5Hz.
Case
Controller parameters
K
p
K
i
K
d
λ
μ
PID
3.9992
2.5840
2.4346
—
—
FOPID
2.8370
3.5545
1.5089
1.0007
0.9924Table 13
Performance indicators withΔf=+0.5Hz.
Index
J
ITAE
_
f
ITAE
_
h
ITAE
_
pe
ITAE
_
y
PID
9.1106
81.7019
4.6269 × 103
1.8867 × 104
7.9778 × 103
FOPID
6.7108
76.1314
1.5369 × 103
0.8624 × 104
1.3310 × 103Table 14
Comparison of performance indexes withΔf=+0.5Hz
Indexes
PID
FOPID
f
pe
h
y
f
pe
h
y
t
f
s
11.22
11.10
11.26
11.06
9.942
9.838
10.00
9.950
t
s
s
15.50
32.87
27.56
42.19
14.12
20.02
22.33
15.56Figure 12
Dynamic response of each output withΔf=+0.5Hz. (a) Frequency. (b) Water hammer pressure. (c) Power. (d) Guide vane opening.
(a)
(b)
(c)
(d)Similarly, it can be obtained from Table13 that the ITAE value under FOPID optimization control is smaller, indicating that the dynamic performance of the unit is more excellent at this time. As can be seen from Figure 12 and Table 14, the opening of the pump turbine and the output power of the synchronous generator are particularly rapid in response speed after frequency disturbance. Under the control of FOPID, the ts value of power is 20.02 seconds, far less than 32.87 s, and the ts value of the opening is less than half of the PID control. It takes less time for each output to reach the stable state, which indicates that FOPID control has better dynamic characteristics and regulation quality.


## 5.3. Robustness Analysis of the Obtained Method

The robustness of PSURS is very important for the stable operation of the unit. The PSURS has the function of water diversion system and water flow inertia constantTw (the value of Tw is hw × Tr) which is an important factor affecting the dynamic performance of PSURS. Increasing the value of Tw will aggravate the fluctuation of water hammer pressure, which can be used to verify the robustness of the proposed method. Therefore, in the test and experiment of the obtained frequency regulation optimization control method, we adjusted the time constant of water flow inertia constant Tw to be 2.1. The simulation experiment of frequency disturbance of −0.5 Hz under low water head and 25% load is carried out. The dynamic response curve of each output is shown in Figure 13. When Tw increases, the controller and control parameters in Table 8 are adopted. Compared with Tw of 1.7, the dynamic response process of frequency, water head, power, and opening of the PSURS increases in a small range of amplitude in the oscillation region, but the overall stability time and steady-state error do not change. The optimal primary frequency control method adopted in this paper is robust to the change of hydraulic parameters.Figure 13
Dynamic response of each output at differentTw with Δf=−0.5Hz. (a) Frequency. (b) Water hammer pressure. (c) Power. (d) Guide vane opening.
(a)
(b)
(c)
(d)


## 6. Conclusion

For the PSU, when it encounters a large frequency disturbance in the low water head and small load condition, it is easy to enter the unstable region and the oscillation occurs. The characteristic curve of the pump turbine is transformed by LCP, and the PSURS refinement model is established based on the seventh-order synchronous generator with the excitation system. This model can describe the nonlinear relationship between the speed-guide vane opening and the flow and torque, the characteristic of the speed governor's limiting and dead zone, and the complex coupling dynamic response law of the hydraulic-mechanical-electrical system of the regulating system is established. Then, based on the refinement model, the dynamic response characteristics of the PSURS under the condition of primary frequency modulation are analyzed by numerical simulation. The results show that when the unit operates at low water head and 25% load, it suffers from frequency disturbance, and the system is easy to enter the “S” zone and then fall into the state of hydraulic, mechanical, and electrical coupling instability. Therefore, the objective function of control parameter optimization considering hydraulic, mechanical, and electrical factors is established. The FOPID controller is used in the optimization control, and the control parameters are optimized by BCGSA. The simulation results show that the FOPID controller shows strong control performance compared with the traditional PID control when PSURS encounters large frequency disturbance in the low water head and small load condition. It effectively improves the dynamic performance of the unit’s transition process. For low water head and small load conditions, the unit has a good inhibitory effect when it encounters a large frequency disturbance and enters an unstable area and violently oscillates.



---

*Source: 6161784-2020-01-11.xml*