A Real-Time Accurate Model and Its Predictive Fuzzy PID Controller for Pumped Storage Unit via Error Compensation

: Model simulation and control of pumped storage unit (PSU) are essential to improve the dynamic quality of power station. Only under the premise of the PSU models reﬂecting the actual transient process, the novel control method can be properly applied in the engineering. The contributions of this paper are that (1) a real-time accurate equivalent circuit model (RAECM) of PSU via error compensation is proposed to reconcile the conﬂict between real-time online simulation and accuracy under various operating conditions, and (2) an adaptive predicted fuzzy PID controller (APFPID) based on RAECM is put forward to overcome the instability of conventional control under no-load conditions with low water head. Respectively, all hydraulic factors in pipeline system are fully considered based on equivalent lumped-circuits theorem. The pretreatment, which consists of improved Suter-transformation and BP neural network, and online simulation method featured by two iterative loops are synthetically proposed to improve the solving accuracy of the pump-turbine. Moreover, the modiﬁed formulas for compensating error are derived with variable-spatial discretization to improve the accuracy of the real-time simulation further. The implicit RadauIIA method is veriﬁed to be more suitable for PSUGS owing to wider stable domain. Then, APFPID controller is constructed based on the integration of fuzzy PID and the model predictive control. Rolling prediction by RAECM is proposed to replace rolling optimization with its computational speed guaranteed. Finally, the simulation and on-site measurements are compared to prove trustworthy of RAECM under various running conditions. Comparative experiments also indicate that APFPID controller outperforms other controllers in most cases, especially low water head conditions. Satisfying results of RAECM have been achieved in engineering and it provides a novel model reference for PSUGS.


Introduction
With the coming of the low-carbon era, the proportion of electricity generated by intermittent renewable energy sources, including photovoltaic power, wind power, and biomass energy, has been growing [1][2][3]. Consequently, there is a greater need for the stability and safety of PSUs, which play the vital roles of generating electricity, peak load shifting, and frequency modulation in the power system [4,5]. In recent years, with the increase in the installed capacity of PSUs and the structural complexity of the hydraulic-mechanical-electrical system, research into transient process and control approaches of PSUs is gaining greater theoretical significance and engineering application value.

Modeling of Full Pipeline System
The dynamic behavior of the transient flow can be described by the mathematical model based on mass and momentum conservation [24]. The solution of hyperbolic partial Equation (1) is inspired by the methods developed in the resolution of the propagation of electrical waves in conductors based on an equivalent scheme representation referred to as the telegraphist's equations [25,26].
The piezometric head H and the discharge Q on the hydraulic side are replaced on the electric side by voltages U and currents I. The analogy between Equation (1) modeling the propagation of pressure waves in hydraulic systems and the telegraphist's equations allows for identifying hydraulic resistance R, hydraulic inductance L, and hydraulic capacitance C. The lineic hydro acoustic parameters are defined as follows: where D is the inner diameter of the pipe, A is the pipe cross-section, c is the wave speed, g is the gravitational acceleration, and λ is the friction coefficient. The capacitive, inductive, and resistive RLC terms correspond to the storage, inertial, and water losses effects, respectively.

Modeling of Full Pipeline System
The dynamic behavior of the transient flow can be described by the mathematical model based on mass and momentum conservation [24]. The solution of hyperbolic partial Equation (1) is inspired by the methods developed in the resolution of the propagation of electrical waves in conductors based on an equivalent scheme representation referred to as the telegraphist's equations [25,26].
The piezometric head H and the discharge Q on the hydraulic side are replaced on the electric side by voltages U and currents I. The analogy between Equation (1) modeling the propagation of pressure waves in hydraulic systems and the telegraphist's equations allows for identifying hydraulic resistance R, hydraulic inductance L, and hydraulic capacitance C. The lineic hydro acoustic parameters are defined as follows: where D is the inner diameter of the pipe, A is the pipe cross-section, c is the wave speed, g is the gravitational acceleration, and λ is the friction coefficient. The capacitive, inductive, and resistive RLC terms correspond to the storage, inertial, and water losses effects, respectively.
Energies 2018, 11, 35 4 of 24 Equation (3) describes the dynamic behavior of the ith element pipe of length dx. Q i /Q i+1 is the input/output discharge, H i /H i+1 is the input/output piezometric heads, and H i+1/2 is the piezometric heads in the middle of the segment. The generalized representation of the length of l pipeline is modeled by n element pipe and dx = l/n. The corresponding spatial discretization is present in Figure 2. Penstock can be derived directly from the equivalent scheme of Figure 2 using Kirchhoff's law as follows: is the piezometric heads in the middle of the segment. The generalized representation of the length of l pipeline is modeled by n element pipe and dx = l/n. The corresponding spatial discretization is present in Figure 2. Penstock can be derived directly from the equivalent scheme of Figure 2 using Kirchhoff's law as follows: For the entire pipeline system, the boundary conditions are the piezometric heads at both ends of the pipe, whereas momentum equations can be merged two by two, introducing the piezometric heads at each node i + 1/2. In addition, the boundary condition of the hydraulic-mechanical junction will be solved in Section 2.2. Using Equations (4) and (5), the implicit matrix form of the system is given by: LsQ n + RQ n + H n+1/2 = H (n−1)+1/2 L 2 sQ n+1 + R 2 Q n+1 + H n+1 = H n+1/2 (4) CsH (n−1)+1/2 + Q n + Q n−1 = 0 For the entire pipeline system, the boundary conditions are the piezometric heads at both ends of the pipe, whereas momentum equations can be merged two by two, introducing the piezometric heads at each node i + 1/2. In addition, the boundary condition of the hydraulic-mechanical junction will be solved in Section 2.2. Using Equations (4) and (5), the implicit matrix form of the system is given by: where the state vectors including the piezometric heads at node i + 1/2; the import and export flow of each element pipe are given by Equation (7). The boundary conditions vector is given by Equation (8).
[A] and [B( → X)] are the system parameter matrices presented in Equation (9). [A 1 ] is the inductance matrix to reflect water inertia; [A 2 ] is the capacitance matrix reflecting the elasticity of the water and wall; [B 1 ] is the impedance matrix, which represents the attrition loss; [B 2 ] and [B 3 ] are the coefficient matrices reflecting the effect on the dynamic characteristics of system with the boundary conditions. It is worth noting that matrix [B( → X)] is a function of the discharge and introduces nonlinear behavior into the system. The discharge Q i in [B( → X)] is retrieved from the Q i−1 [27]: According to the actual piping layout of the power plant, the implicit nonlinear system of PSU is interpreted as a topological network comprising T-type equivalent circuits, as shown in Figure 3. Each hydraulic component, such as the surge tank, valve, and reservoir, can be defined by the equivalent circuit in [27]; the pump-turbine model is described below.
where the state vectors including the piezometric heads at node i + 1/2; the import and export flow of each element pipe are given by Equation (7). The boundary conditions vector is given by Equation (8) [ ] According to the actual piping layout of the power plant, the implicit nonlinear system of PSU is interpreted as a topological network comprising T-type equivalent circuits, as shown in Figure 3. Each hydraulic component, such as the surge tank, valve, and reservoir, can be defined by the equivalent circuit in [27]; the pump-turbine model is described below.

Modeling of Pump-Turbine
The "S" area of the characteristic curves presents uneven distribution, crossing, and aggregation, which causes difficulties when modeling the pump-turbine [28,29]. The simplified transmission coefficient model of the pump-turbine [30] is widely applied in the specified conditions but its accuracy cannot be guaranteed in different operating conditions. The pump-turbine model based on the characteristic curves is hence being extensively used to acquire accurate simulations in all operating conditions. Therefore, pretreatment of the characteristic curves is of vital importance for the modeling of the pump-turbine.
The interpolation model of the pump-turbine based on improved Suter transformation and the BP neural network has been put forward to solve the above shortcomings. Firstly, the improved Suter transformation method is applied to eliminate the crossing, aggregation, and single input-multiple output problems during the process of interpolation. Then, a BP neural network is constructed to enrich the lines of gate opening, data extension, and bad point correction for the treated curve. Finally, the Lagrange interpolation algorithm of double variables-three points form is used to solve the pump-turbine model. Improved Suter transformation adopts dimensionless similarity parameters to describe the characteristic curves as follows: where, a, q, h, m, y, x indicate the relative values of unit speed, unit discharge, unit water head, unit torque, unit guide vane opening and flow angle, respectively. M 11max is the maximum unit torque, M 11r is the rated unit torque, k 1 > |M 11max |/M 11r , k 2 = 0.5-1.2, C y = 0.1-0.3, C h = 0.4-0.6. As shown in Figure 4 (black curves), it is obvious that the improved Suter transformation conquers the drawbacks of crossing, aggregation, and multi-value at "S" area. Nevertheless, the transformed curves still have some shortcomings, as follows: (1) the abscissa of the curves is not aligned and the distribution of the data points in each curve is uneven; (2) unqualified data points exist in the hump area and "S" area because of measurement errors; (3) the measured curves from the manufacturer are too few to meet the requirements of interpolation.

Modeling of Pump-Turbine
The "S" area of the characteristic curves presents uneven distribution, crossing, and aggregation, which causes difficulties when modeling the pump-turbine [28,29]. The simplified transmission coefficient model of the pump-turbine [30] is widely applied in the specified conditions but its accuracy cannot be guaranteed in different operating conditions. The pump-turbine model based on the characteristic curves is hence being extensively used to acquire accurate simulations in all operating conditions. Therefore, pretreatment of the characteristic curves is of vital importance for the modeling of the pump-turbine.
The interpolation model of the pump-turbine based on improved Suter transformation and the BP neural network has been put forward to solve the above shortcomings. Firstly, the improved Suter transformation method is applied to eliminate the crossing, aggregation, and single input-multiple output problems during the process of interpolation. Then, a BP neural network is constructed to enrich the lines of gate opening, data extension, and bad point correction for the treated curve. Finally, the Lagrange interpolation algorithm of double variables-three points form is used to solve the pump-turbine model. Improved Suter transformation adopts dimensionless similarity parameters to describe the characteristic curves as follows: where, a, q, h, m, y, x indicate the relative values of unit speed, unit discharge, unit water head, unit torque, unit guide vane opening and flow angle, respectively. 11max M is the maximum unit torque,  Figure 4 (black curves), it is obvious that the improved Suter transformation conquers the drawbacks of crossing, aggregation, and multi-value at "S" area. Nevertheless, the transformed curves still have some shortcomings, as follows: (1) the abscissa of the curves is not aligned and the distribution of the data points in each curve is uneven; (2) unqualified data points exist in the hump area and "S" area because of measurement errors; (3) the measured curves from the manufacturer are too few to meet the requirements of interpolation.   the coupling relation between W H and W M. In the training process of the BP neural network, the maximum number of iterations is 500, the number of each hidden layer nodes is 15, the learning rate is 0.1, and the training objective error is 2 × 10 −6 . As shown in Figure 4, the inherent law of characteristic can be generalized from the existing data sequence thanks to the strong self-learning and predictive ability of the BP neural network. It is necessary to enrich the lines of guide vane opening that is hard to get from the manufacturer. The maximum relative errors of the prediction curves are only 0.038 and 0.04. As the magnified image in Figure 4 shows, the relative errors have mainly been concentrated in anti-pump, pump braking, and turbine braking conditions. This indicates that the BP neural network does not repeat the curves indiscriminately, but corrects unqualified data points appropriately by understanding the inherent law of characteristic curves. It is helpful to decrease the measuring error and improve the accuracy of interpolation.
In order to make full use of the information from the curves, the Lagrange interpolation algorithm of double variables-three points form is adopted to solve the treated curves, as shown in Equation (11). z(x, y) is unknown points, i.e., W H and W M, x p -x p+2 and y q -y q+2 are relative flow angles and guide vane openings in the interpolation area.
The conventional method for solving water head and torque is given by Equation (12) [31]. One can notice that the relative working head h n+1 at the n + 1 moment is calculated to be zero by Equation (12) when the relative speed and relative flow at the n moment are zero. This is not consistent with the actual water head of PSU, which may lead to unstable and discontinuous simulation in the start-up condition and running condition conversion. Therefore, the improved method is as in Equation (13), which takes into account when the guide vane opening is very small or zero. Hence, it guarantees the consistency between the calculation and the actual operation of the pump-turbine.
The flowchart of the pump-turbine model is illustrated in Figure 5. δ 1 and δ 2 are the values of iterative precision, and N 1 and N 2 are the maximum iterations. The model of the motor-generator is introduced in [9]. Two iterative loops of speed and discharge are applied to ensure the accuracy of the interpolation. Simulation computing is conducted in the Matlab R2013a on a computer with 2.40 GHz Inter(R) Core(TM) i3-2370M CPU, and 4 GB of RAM. It is worth noting that the BP method is only involved in constructing accurate characteristic curves and did not participate in online simulating. The pump-turbine model can not only benefit from strong self-learning and prediction performance of the BP neural network in pretreatment of characteristic curves, but also saves simulation time by avoiding the computational burden of BP functions.

Nonlinear Servo-Mechanism
All non-linear factors are taken into account for the servo-mechanism, as shown in Figure 6, including main distribution valve, main servomotor, the dead zone of the servo-device, the saturation of the pressure regulating valve, the rate limit of the main servomotor, and the output limit of guide vane opening. 0 k is the gain coefficient of servo-mechanism, 1 y T is the main distribution valve response time, and y T is the main servomotor response time.

Space-Time Discrete Analysis of RAECM
The Courant's stability criteria linking the space and time discretization, pipe space step dx and simulation time step dt , respectively, through the wave speed c, is given by

Nonlinear Servo-Mechanism
All non-linear factors are taken into account for the servo-mechanism, as shown in Figure 6, including main distribution valve, main servomotor, the dead zone of the servo-device, the saturation of the pressure regulating valve, the rate limit of the main servomotor, and the output limit of guide vane opening. k 0 is the gain coefficient of servo-mechanism, T y1 is the main distribution valve response time, and T y is the main servomotor response time.

Nonlinear Servo-Mechanism
All non-linear factors are taken into account for the servo-mechanism, as shown in Figure 6, including main distribution valve, main servomotor, the dead zone of the servo-device, the saturation of the pressure regulating valve, the rate limit of the main servomotor, and the output limit of guide vane opening. 0 k is the gain coefficient of servo-mechanism, 1 y T is the main distribution valve response time, and y T is the main servomotor response time.

Space-Time Discrete Analysis of RAECM
The Courant's stability criteria linking the space and time discretization, pipe space step dx and simulation time step dt , respectively, through the wave speed c, is given by

Space-Time Discrete Analysis of RAECM
The Courant's stability criteria linking the space and time discretization, pipe space step dx and simulation time step dt, respectively, through the wave speed c, is given by dt < dx/c = dT, where dT is the time step of space discretization. The MOC with the limitation of Courant's stability condition is widely used in simulations of unsteady flow in a pipeline with small space steps [32]. However, it is Energies 2018, 11, 35 9 of 24 unsuitable for calculations with large space steps because of the limitations of strict adherence to the time-step-space relationship.

Discrete Analysis of Variable Space-Step
Different from MOC, the time step of space discretization is not restricted by the simulation time step in the RAECM. The dimensions of the discrete system and time consumption decrease significantly with the increase in dT. Consequently, RAECM can guarantee real-time simulation of PSUGS by reasonably increasing dT under the premise of ensuring accuracy.
In this section, we set the time step of space discretization dT to be 0.02 s, 0.05 s, 0.1 s, and 0.5 s independently. The equivalent circuit model of variable space-steps is compared with MOC under 100% load rejection operation, as shown in Table 1 and Figure 7. With the increase of dT, the maximum value of volute pressure decreases and the minimum value of pressure at the draft tube increases gradually. Meanwhile, the friction and elasticity of the unit pipe section from the center to the end must be considered. Accordingly, the modified formulas of water head for compensating error are proposed in RAECM. If the correction node is on the right side of the unit pipe, the correcting formula is given by Equation (14); otherwise, it is as given by Equation (15). Here, H R and H L are the corrected water head of the monitoring node in jth section unit pipe. H j,C is the water head of the center node in jth section unit pipe. Q t j,L and Q t−1 j,L are the discharge in the left of jth section unit pipe at t and t − 1. Q t j,R and Q t−1 j,R are the discharge in the right of jth section unit pipe at t and t − 1.
As demonstrated in Table 1 and Figure 7, the calculated water pressure in volute and draft tube before and after the correction are similar in the equivalent circuit model with small-scale spatial discretization (dT < 0.05 s), which means very small truncation error. That is because the number of unit pipes is very large, and hence accurate calculation of the water pressure can be achieved. However, it is worth noting that the truncation error in the equivalent circuit model with large-scale spatial discretization (dT ≥ 0.05 s) is too large to be neglected. Therefore, truncation error compensation is essential to RAECM for the purpose of precisely reflecting the actual transient process. Table 1. Comparison between RAECM with variable space-step discretization and MOC.

Model
Step

Discrete Analysis of Variable Time Step
An existing algorithm of fourth order explicit R-K has been applied to solve the global set of differential equations [33]. In this paper, the implicit RadauIIA algorithm is developed to solve nonlinear ordinary differential equations. The implicit RadauIIA formula is depicted as: Here, n y and 1 + n y are the outputs of system at n and n+1 moment, t Δ is time step of simulation, i c , i b , ij a are got from Butcher matrix that 1 12 1 / 12 a = − , 21 3 / 4 a = , 22 Equation (6) can be reorganized as follows: where [E] is the characteristic matrix of RAECM. The eigenvalues of [E] are given by formula (18), k s reflects the system stability. The real part k σ is damping coefficient, the imaginary part k w represents the oscillation frequency of the system. The real eigenvalues represent the damping mode corresponding to the rigid water hammer of PSU and the complex eigenvalues represent the vibration mode corresponding to the elastic water hammer. Taking load rejection condition as an example, all the eigenvalues of RAECM with dT = 0.02 s, dT = 0.05 s and dT = 0.5 s are represented in the complex plan in Figure 8. With the increase of dT, the number of characteristic points decreases and the computational efficiency of RAECM is improved. All the eigenvalues are on the left hand side of the complex plan, corresponding to the negative damping of system [33], i.e., the model of PSU is stable. All the complex eigenvalues are symmetrically distributed in the complex plan, which indicated that the reciprocating motion of the water hammer wave in the elastic pipe.

Discrete Analysis of Variable Time Step
An existing algorithm of fourth order explicit R-K has been applied to solve the global set of differential equations [33]. In this paper, the implicit RadauIIA algorithm is developed to solve nonlinear ordinary differential equations. The implicit RadauIIA formula is depicted as: Here, y n and y n+1 are the outputs of system at n and n+1 moment, ∆t is time step of simulation, Equation (6) can be reorganized as follows: → .
where [E] is the characteristic matrix of RAECM. The eigenvalues of [E] are given by formula (18), s k reflects the system stability. The real part σ k is damping coefficient, the imaginary part w k represents the oscillation frequency of the system. The real eigenvalues represent the damping mode corresponding to the rigid water hammer of PSU and the complex eigenvalues represent the vibration mode corresponding to the elastic water hammer. Taking load rejection condition as an example, all the eigenvalues of RAECM with dT = 0.02 s, dT = 0.05 s and dT = 0.5 s are represented in the complex plan in Figure 8. With the increase of dT, the number of characteristic points decreases and the computational efficiency of RAECM is improved. All the eigenvalues are on the left hand side of the complex plan, corresponding to the negative damping of system [33], i.e., the model of PSU is stable. All the complex eigenvalues are symmetrically distributed in the complex plan, which indicated that the reciprocating motion of the water hammer wave in the elastic pipe.
Energies 2017, 10, x FOR PEER REVIEW 11 of 24 , 1 ,2, In order to ensure the accuracy of the RAECM in all operating conditions, it is necessary to analyze the stability of different solution algorithms under variable time step discretizations. The stability functions of fourth-order explicit R-K and implicit RadauIIA algorithm, respectively, are defined as follows: is the stability criteria of Equation (17) with numerically solving [33].
Following variable time step dt = 0.02 s, 0.05 s, 0.1 s, the characteristic points of RAECM with dT = 0.05 s are illustrated in Figure 9. At ≤ dt dT , i.e., Courant's criteria, the fourth-order R-K and the implicit RadauIIA method can guarantee the stability of the system, as shown in Figure 9a,b. At ≥ dt dT , the value of z is enlarged with the increase of dt, which makes the distribution of characteristic points in the complex plan exceed the stability region of the fourth-order R-K algorithm and causes solution instability. Nevertheless, the implicit RadauIIA algorithm can still ensure the absolute convergence of RAECM with large time step discretization. Img ω In order to ensure the accuracy of the RAECM in all operating conditions, it is necessary to analyze the stability of different solution algorithms under variable time step discretizations. The stability functions of fourth-order explicit R-K and implicit RadauIIA algorithm, respectively, are defined as follows: where, z = dt × s k . R(|z|) ≤ 1 is the stability criteria of Equation (17) with numerically solving [33]. Following variable time step dt = 0.02 s, 0.05 s, 0.1 s, the characteristic points of RAECM with dT = 0.05 s are illustrated in Figure 9. At dt ≤ dT, i.e., Courant's criteria, the fourth-order R-K and the implicit RadauIIA method can guarantee the stability of the system, as shown in Figure 9a,b. At dt ≥ dT, the value of z is enlarged with the increase of dt, which makes the distribution of characteristic points in the complex plan exceed the stability region of the fourth-order R-K algorithm and causes solution instability. Nevertheless, the implicit RadauIIA algorithm can still ensure the absolute convergence of RAECM with large time step discretization.
Following variable time step dt = 0.02 s, 0.05 s, 0.1 s, the characteristic points of RAECM with dT = 0.05 s are illustrated in Figure 9. At ≤ dt dT , i.e., Courant's criteria, the fourth-order R-K and the implicit RadauIIA method can guarantee the stability of the system, as shown in Figure 9a,b. At ≥ dt dT , the value of z is enlarged with the increase of dt, which makes the distribution of characteristic points in the complex plan exceed the stability region of the fourth-order R-K algorithm and causes solution instability. Nevertheless, the implicit RadauIIA algorithm can still ensure the absolute convergence of RAECM with large time step discretization.

Implicit RadauIIA vs. Fourth-Order Explicit R-K
The black and grey areas shown in Figure 10 indicate the unconditionally stable region of RadauIIA and fourth-order explicit R-K; the darker the margin, the more stable the region. It can be easily observed that the unconditional stable region of fourth-order explicit R-K has a very limited distribution on the left-hand side of the complex plan. On the contrary, the unconditional stable region of implicit RadauIIA covers the whole complex plain except for a small area on the right-hand side of the complex axis. Especially when the guide vane closed from full opening to zero with fast speed under extreme conditions such as load rejection and pump outage, the imaginary part of s k is very large because of the high-frequency pressure fluctuation, which is hazardous to the stability of the fourth-order explicit R-K for RAECM. However, the implicit RadauIIA method can overcome this obstacle and reflect the actual transition process of all operating conditions of PSU.

Implicit RadauIIA vs. Fourth-Order Explicit R-K
The black and grey areas shown in Figure 10 indicate the unconditionally stable region of RadauIIA and fourth-order explicit R-K; the darker the margin, the more stable the region. It can be easily observed that the unconditional stable region of fourth-order explicit R-K has a very limited distribution on the left-hand side of the complex plan. On the contrary, the unconditional stable region of implicit RadauIIA covers the whole complex plain except for a small area on the right-hand side of the complex axis. Especially when the guide vane closed from full opening to zero with fast speed under extreme conditions such as load rejection and pump outage, the imaginary part of k s is very large because of the high-frequency pressure fluctuation, which is hazardous to the stability of the fourth-order explicit R-K for RAECM. However, the implicit RadauIIA method can overcome this obstacle and reflect the actual transition process of all operating conditions of PSU. Theoretically, PSUGS is a typical implicit system because of its time-varying feature, nonlinearity, and non-minimum phase characteristic. Therefore, the implicit RadauIIA method is more suitable for PSUGS, while the fourth-order explicit R-K is unstable for the implicit system of PSU in most conditions, especially when the step of space-time discretization is large and extreme.

Adaptive Predicted Fuzzy PID Controller
PID and improved PID controllers are based on current and past system signals, but lack the prediction of the future states of the system. MPC relies heavily on the system model, and a non-analytical nonlinear model such as MOC is inappropriate in the prediction control theory framework [34]. Aiming to make up for the above shortcomings, an adaptively predicted fuzzy PID based on the integration of FPID and MPC is proposed to regulate PSU accurately in diverse operating conditions with its computational speed guaranteed.

Fuzzy PID Controller
The FPID controller consists of a basic PID controller and a fuzzy inference machine, which adaptively adjusts the parameters of the PID controller as shown in Figure 11. The process of fuzzy control mainly comprises fuzzification, fuzzy rule base, membership functions used in the rules, reasoning mechanism by use of fuzzy logic operators, and defuzzification [34]. Theoretically, PSUGS is a typical implicit system because of its time-varying feature, nonlinearity, and non-minimum phase characteristic. Therefore, the implicit RadauIIA method is more suitable for PSUGS, while the fourth-order explicit R-K is unstable for the implicit system of PSU in most conditions, especially when the step of space-time discretization is large and extreme.

Adaptive Predicted Fuzzy PID Controller
PID and improved PID controllers are based on current and past system signals, but lack the prediction of the future states of the system. MPC relies heavily on the system model, and a non-analytical nonlinear model such as MOC is inappropriate in the prediction control theory framework [34]. Aiming to make up for the above shortcomings, an adaptively predicted fuzzy PID based on the integration of FPID and MPC is proposed to regulate PSU accurately in diverse operating conditions with its computational speed guaranteed.

Fuzzy PID Controller
The FPID controller consists of a basic PID controller and a fuzzy inference machine, which adaptively adjusts the parameters of the PID controller as shown in Figure 11. The process of fuzzy control mainly comprises fuzzification, fuzzy rule base, membership functions used in the rules, reasoning mechanism by use of fuzzy logic operators, and defuzzification [34].

Variation Factor of PID Parameters
Based on fuzzy rules, the PID parameters can be tuned dynamically. Nevertheless, the rapid change of PID parameters may lead to vibration in the control outputs. The variation factor ( ) k γ is proposed to restrain the change of PID parameters according to the frequency deviation to ensure the stability of control process, which is defined as follows: where e(k), e(k−1) and e(k−2) are the speed errors of PSUGS at time k, k−1, and k−2. The variation factor follows the linear decreasing rule as in Equation (21). It can gradually weaken the influence of fluctuation of PID parameters for the control quality. The real-time deviations of PID parameters [ p K′ , i K ′ , d K ′ ] are deduced with a fuzzy inference engine. Thus, the PID parameters of the current operating condition at k-step can be computed as follows: Figure 11. Block diagram of the APFPID controller.
A typical two-input, three-output fuzzy controller is applied in this paper. Speed tracking error and change rate of tracking error, i.e., e and e c , are the inputs of the controller; the deviation of PID parameters K p , K i , and K d are the outputs, as shown in Figure 11. Mamdani-type inferencing is used to describe the fuzzy rules. The variation of outputs is deduced by the fuzzy If-Then rules as follows:

Variation Factor of PID Parameters
Based on fuzzy rules, the PID parameters can be tuned dynamically. Nevertheless, the rapid change of PID parameters may lead to vibration in the control outputs. The variation factor γ(k) is proposed to restrain the change of PID parameters according to the frequency deviation to ensure the stability of control process, which is defined as follows: where e(k), e(k−1) and e(k−2) are the speed errors of PSUGS at time k, k−1, and k−2. The variation factor follows the linear decreasing rule as in Equation (21). It can gradually weaken the influence of fluctuation of PID parameters for the control quality. The real-time deviations of PID parameters [K p , K i , K d ] are deduced with a fuzzy inference engine. Thus, the PID parameters of the current operating condition at k-step can be computed as follows:

Rolling Prediction
The predictive functions are always applied in predictive control to obtain the future information [34]. However, an appropriate predictive function is hard to find to match the complex PSU under various conditions. In addition, traditional predictive control only retains the first predictive control output and abandons the predicted value the rest of the time because of the online computation for rolling optimization. In order to overcome the above problems, rolling prediction is proposed to replace rolling optimization because the RAECM can exactly match the actual transient process in various conditions. The incremental sequences of control signals in the future N p -step can be calculated in the sampling period of industrial frequency owing to the real-time performance of RAECM.
The rolling prediction process is described by Figure 11. According to the real-time control signal u(k/k) at moment k, the real-time PID parameters, and the current rotational speed of the system x(k/k), we can get predictive control output u(k + 1/k) of the (k + 1)th step at the k moment. The increment of the (k + 1)th step prediction control output at the k moment, i.e., ∆u(k + 1/k) is defined as follows: Without loss of generality, ∆u(k + i/k) is expressed as follows: The increment sequences of control signals in the future N p -step are obtained by Equation (24) as follows: ∆u(k + 1/k), i = 0, i = 0, 1, · · · , N p − 1 . It should be noted that the rolling sequences not only contain the past and current information but also possess the predictive trends of control process in the future. In the practical control process, the current control increment is more important than the forecasting value, and the effect of predictive control increments gradually abates with the increase of predictive step i. Therefore, a weighted formula of nonlinear attenuation is proposed for control signal increment sequences, as follows: where the attenuation coefficient λ = N p −i N p , N p is the predictive horizon. The real-time control output reasonably contains the possible state deviation of the future system and subsequent control law using Equation (26). In addition, the softening control signal of APFPID is conducive to system stability with suddenly large external disturbance.
Due to limited space, the flowchart of the APFPID controller is depicted in Figure 12.

Experiments and Analysis of Results
To verify the effectiveness, accuracy, and stability of the RAECM attached to a APFPID controller, comparative experiments and analyses are conducted in this section. According to the actual layout of the pipeline system and the real physical parameters of the pumped storage plant in China (see Table A4 and A5), the RAECM is designed as shown in Section 2. The RAECM via error compensation is dispersed with 0.5 dT s = , as described in Section 3. The APFPID controller is constructed as introduced in Section 4.

Pump Outage Condition
All the initial simulation settings are the same as the actual settings in this condition: upper reservoir water level 731.9 m As demonstrated in Figure 13, the simulation reflects the real operating conditions well, and the trend of parametric curves with simulation is consistent with the practical process. From further

Experiments and Analysis of Results
To verify the effectiveness, accuracy, and stability of the RAECM attached to a APFPID controller, comparative experiments and analyses are conducted in this section. According to the actual layout of the pipeline system and the real physical parameters of the pumped storage plant in China (see Tables A4 and A5), the RAECM is designed as shown in Section 2. The RAECM via error compensation is dispersed with dT = 0.5 s, as described in Section 3. The APFPID controller is constructed as introduced in Section 4.

Pump Outage Condition
All the initial simulation settings are the same as the actual settings in this condition: upper reservoir water level H u = 731.9 m, lower reservoir water level H d = 167.7 m, discharge of pump-turbine Q = −49.37 m 3 /s, unit power P = 279.93 MW, and the opening is 80% of the rated guide van opening. The pump condition cedes to pump outage condition when the breaker is suddenly tripped at 7.82 s. The guide vane opening values are sampled on-site with the HBM data acquisition system (QuantumXMX840A-P), and then inputted to the RAECM. Then the different physical quantities of the transient process can be obtained.
As demonstrated in Figure 13, the simulation reflects the real operating conditions well, and the trend of parametric curves with simulation is consistent with the practical process. From further analysis, as shown in Table 2, it is evident that the extreme indices of pressure at volute and draft tube simulated with RAECM are closer to the on-site values compared with MOC. The running track of pump outage condition with RAECM passes through the pump area, the pump braking area, the turbine area, and a small part of the turbine braking area, as shown in Figure 14. For the simulated pressure, shown in Figure 13b,c, there is a small static deviation from the measurements, which will be discussed in the next section. As shown in Figure 13b, during the operating period with zero guide vane opening, the pressure oscillations gradually attenuate and are smaller than those of MOC owing to the implicit RadauIIA method. Meanwhile, unattenuated oscillations occur in the volute with MOC. As shown in the dashed frame in Figure 13c, the pulsating pressure in the draft tube might be difficult to compare to modeled values, due to the swirl not being considered. It is worth mentioning that the RAECM will take into account the excitation source of the pump-turbine as well as the effect of the vortex rope compressibility in the future work. That cannot be done by the MOC [35]. analysis, as shown in Table 2, it is evident that the extreme indices of pressure at volute and draft tube simulated with RAECM are closer to the on-site values compared with MOC. The running track of pump outage condition with RAECM passes through the pump area, the pump braking area, the turbine area, and a small part of the turbine braking area, as shown in Figure 14. For the simulated pressure, shown in Figure 13b,c, there is a small static deviation from the measurements, which will be discussed in the next section. As shown in Figure 13b, during the operating period with zero guide vane opening, the pressure oscillations gradually attenuate and are smaller than those of MOC owing to the implicit RadauIIA method. Meanwhile, unattenuated oscillations occur in the volute with MOC. As shown in the dashed frame in Figure 13c, the pulsating pressure in the draft tube might be difficult to compare to modeled values, due to the swirl not being considered. It is worth mentioning that the RAECM will take into account the excitation source of the pump-turbine as well as the effect of the vortex rope compressibility in the future work. That cannot be done by the MOC [35].

Start-Up and No-Load Operations
PID controller has good performance for PSU under a conventional water head. However, PSU is apt to fall into the "S" characteristic area under the middle-low water head, which may lead to intensive oscillations [9]. Start-up to no-load conditions under a low water head are studied not only to verify the effectiveness of the APFPID controller but also to test the stability and accuracy of RAECM in small disturbance conditions.

Study of APFPID Controller Performance
The operating parameters of PSUGS under start-up and no-load operating conditions are given in the Table A5. An APFPID controller together with PID and FPID controller is applied to the PSUGS for better dynamic performance. In order to perform a fair comparison in fitness evaluation, the parameters , , p i d K K K of PID controller are set to 1.5, 0.2, 0.1 via optimization of GSA [17]. Then they are directly set as the initial parameters of the APFPID and FPID controllers. The low working water head is set to 535 m, the no-load guide vane opening is 18.4% of the rated value, and the predictive horizon p N is set to 5. When the governor receives the start-up command, the guide vane opens to no-load opening with the fastest speed and holds this opening angle until the rotational speed reaches 90% of the rated speed. Then, the controller is switched to closed-loop control. The simulation time is set to 100 s. The transient process and performance indices with the considered controllers are illustrated in Table 3 and Figure 15. It is obvious that the APFPID controller has better dynamic and static performance than the other controllers in start-up and no-load conditions under a low water head. From Table 3 and Figure 15, since the operation point is closer to the "S" characteristic area at a low water head (i.e., 535 m), the PID controller still has bad index values with more oscillation times and large overshoot under low water head conditions, although the PID parameters have been optimized by GSA. The adjusting time of the FPID controller is longer than others, which may be due to the fact that the initial PID parameters are not individually optimized. The APFPID controller shows superior control performance over PID and FPID controllers with less overshoot and more rapid adjustment because of rolling prediction. Figure 14b helps to confirm the above analysis; the operation point of PSU is pulled back into the no-load operation area after crossing the runaway speed curve under start-up to no-load conditions with a low water head.

Start-Up and No-Load Operations
PID controller has good performance for PSU under a conventional water head. However, PSU is apt to fall into the "S" characteristic area under the middle-low water head, which may lead to intensive oscillations [9]. Start-up to no-load conditions under a low water head are studied not only to verify the effectiveness of the APFPID controller but also to test the stability and accuracy of RAECM in small disturbance conditions.

Study of APFPID Controller Performance
The operating parameters of PSUGS under start-up and no-load operating conditions are given in the Table A5. An APFPID controller together with PID and FPID controller is applied to the PSUGS for better dynamic performance. In order to perform a fair comparison in fitness evaluation, the parameters K p , K i , K d of PID controller are set to 1.5, 0.2, 0.1 via optimization of GSA [17]. Then they are directly set as the initial parameters of the APFPID and FPID controllers. The low working water head is set to 535 m, the no-load guide vane opening is 18.4% of the rated value, and the predictive horizon N p is set to 5. When the governor receives the start-up command, the guide vane opens to no-load opening with the fastest speed and holds this opening angle until the rotational speed reaches 90% of the rated speed. Then, the controller is switched to closed-loop control. The simulation time is set to 100 s.
The transient process and performance indices with the considered controllers are illustrated in Table 3 and Figure 15. It is obvious that the APFPID controller has better dynamic and static performance than the other controllers in start-up and no-load conditions under a low water head. From Table 3 and Figure 15, since the operation point is closer to the "S" characteristic area at a low water head (i.e., 535 m), the PID controller still has bad index values with more oscillation times and large overshoot under low water head conditions, although the PID parameters have been optimized by GSA. The adjusting time of the FPID controller is longer than others, which may be due to the fact that the initial PID parameters are not individually optimized. The APFPID controller shows superior control performance over PID and FPID controllers with less overshoot and more rapid adjustment because of rolling prediction. Figure 14b helps to confirm the above analysis; the operation point of PSU is pulled back into the no-load operation area after crossing the runaway speed curve under start-up to no-load conditions with a low water head.

Robustness Analysis of the APFPID Controller
Hydraulic factor is the primary factor influencing PSUGS control performance under no-load conditions. The low water head can cause poor robustness of the controller, as mentioned above. Hence, various water heads, i.e., 540 m, 530 m, and 527 m, are considered in order to verify the robustness of the APFPID controller. From Figure 16, with the decrease of water head, the APFPID controller can ensure the stability of the system better than the PID controller, which indicates that the proposed APFPID controller shows sufficient robustness to hydraulic parameter variation of the system under no-load conditions.

Robustness Analysis of the APFPID Controller
Hydraulic factor is the primary factor influencing PSUGS control performance under no-load conditions. The low water head can cause poor robustness of the controller, as mentioned above. Hence, various water heads, i.e., 540 m, 530 m, and 527 m, are considered in order to verify the robustness of the APFPID controller. From Figure 16, with the decrease of water head, the APFPID controller can ensure the stability of the system better than the PID controller, which indicates that the proposed APFPID controller shows sufficient robustness to hydraulic parameter variation of the system under no-load conditions.

Robustness Analysis of the APFPID Controller
Hydraulic factor is the primary factor influencing PSUGS control performance under no-load conditions. The low water head can cause poor robustness of the controller, as mentioned above. Hence, various water heads, i.e., 540 m, 530 m, and 527 m, are considered in order to verify the robustness of the APFPID controller. From Figure 16, with the decrease of water head, the APFPID controller can ensure the stability of the system better than the PID controller, which indicates that the proposed APFPID controller shows sufficient robustness to hydraulic parameter variation of the system under no-load conditions.

No-Load Stability Analysis of RAECM
The no-load operation point is far from the designed operation point, which is close to the "S" area. In addition to satisfying the requirements of large disturbance conditions, the stability and accuracy of RAECM must be analyzed in the no-load condition. The dynamic response and performance indicators of the no-load operating condition with diverse model and filtered on-site measurement are obtained in Figure 17 and Table 4. From Figure 17a and Table 4, it is revealed that the real pressure of volute after filtering is in the state of fluctuation. The RAECM considering the elasticity of wall is consistent with the on-site measurement, and is more tolerant with ample fluctuation than MOC. From Figure 17b and Table 4, it is obvious that the rigid water hammer model fails to reflect the small fluctuation characteristic of rotational speed under no-load condition, while RAECM is able to roughly reflect the tendency of the filtered on-site measured speed, which demonstrates that RAECM can simulate the real operating conditions well.

No-Load Stability Analysis of RAECM
The no-load operation point is far from the designed operation point, which is close to the "S" area. In addition to satisfying the requirements of large disturbance conditions, the stability and accuracy of RAECM must be analyzed in the no-load condition. The dynamic response and performance indicators of the no-load operating condition with diverse model and filtered on-site measurement are obtained in Figure 17 and Table 4. From Figure 17a and Table 4, it is revealed that the real pressure of volute after filtering is in the state of fluctuation. The RAECM considering the elasticity of wall is consistent with the on-site measurement, and is more tolerant with ample fluctuation than MOC. From Figure 17b and Table 4, it is obvious that the rigid water hammer model fails to reflect the small fluctuation characteristic of rotational speed under no-load condition, while RAECM is able to roughly reflect the tendency of the filtered on-site measured speed, which demonstrates that RAECM can simulate the real operating conditions well.

Real-Time Analysis and Engineering Application
The RAECM proposed in this paper has been successfully embedded into the PSUGS performance test instrument for a pumped storage plant in Jiangxi Province, China. As shown in Figure 18, the semi-physical simulation equipment is composed of DSP and a fully isolated integrated array device, which has online simulation and on-site recording functions. For the online simulation, the equipment is treated as a digital PSU connected with the governor through digital

Real-Time Analysis and Engineering Application
The RAECM proposed in this paper has been successfully embedded into the PSUGS performance test instrument for a pumped storage plant in Jiangxi Province, China. As shown in Figure 18, the semi-physical simulation equipment is composed of DSP and a fully isolated integrated array device, which has online simulation and on-site recording functions. For the online simulation, the equipment is treated as a digital PSU connected with the governor through digital I/O interfaces and an A/D conversion module. It gets switch inputs and opening signals from the governor and outputs frequency to the governor, which constitutes a closed-loop system.  In each step of the simulation test of PSUGS, the model unit operation time is defined as the total time consumption including the reception of guide vane opening, the calculation of RAECM, and the output of the rotational speed signal from RAECM. Simulation computing is conducted in DSP (TMS320F28335, produced by TI Company located in Dallax, TX, USA). A total of 15,000 steps are taken in the simulation and the average model unit operation time of RAECM is 0.68 ms. In the engineering applications, the sampling interval or the adjustment period is usually determined by a power frequency cycle of 20 ms, so the RAECM can fully meet the real-time simulation requirement of engineering.
By simulating the real PSU, the dynamic real-time test of PSUGS can be carried out under various operating conditions. This is useful in the commissioning of a power plant because it enables quick testing to reduce the time and cost of tuning system parameters and control settings. Satisfying results for RAECM have been achieved in practical application; this lays the foundation for the application of the APFPID controller in engineering.

Discussion
The reliability and accuracy of RAECM are verified under various operating conditions, compared with on-site measurement, MOC, and a simplified model. In the pump outage condition, although the transient process of RAECM goes through extreme operating areas (hump area and "S" zone), the simulation results coincide with the actual measurements. This indicates that RAECM is effective in various operating areas and stable in extreme conditions. In the start-up and no-load conditions, RAECM is able to reflect the undulation of the filtered on-site measured rotational speed and water pressure owing to the influence of hydraulic factors such as the actual pipeline layout, the elasticity of wall, the inertia of water, and the friction loss. Furthermore, there is a small static deviation between the simulation and the on-site measurement. It might be a consequence of the following aspects: (1) the difference between model test data and the real unit characteristic curve; (2) the imprecise parameters of the pipeline system such as pipe roughness or friction coefficient. In Section 5.1, RAECM reflects the real transient process better than MOC under zero guide vane opening conditions owing to the implicit RadauII A method. In each step of the simulation test of PSUGS, the model unit operation time is defined as the total time consumption including the reception of guide vane opening, the calculation of RAECM, and the output of the rotational speed signal from RAECM. Simulation computing is conducted in DSP (TMS320F28335, produced by TI Company located in Dallax, TX, USA). A total of 15,000 steps are taken in the simulation and the average model unit operation time of RAECM is 0.68 ms. In the engineering applications, the sampling interval or the adjustment period is usually determined by a power frequency cycle of 20 ms, so the RAECM can fully meet the real-time simulation requirement of engineering.
By simulating the real PSU, the dynamic real-time test of PSUGS can be carried out under various operating conditions. This is useful in the commissioning of a power plant because it enables quick testing to reduce the time and cost of tuning system parameters and control settings. Satisfying results for RAECM have been achieved in practical application; this lays the foundation for the application of the APFPID controller in engineering.

Discussion
The reliability and accuracy of RAECM are verified under various operating conditions, compared with on-site measurement, MOC, and a simplified model. In the pump outage condition, although the transient process of RAECM goes through extreme operating areas (hump area and "S" zone), the simulation results coincide with the actual measurements. This indicates that RAECM is effective in various operating areas and stable in extreme conditions. In the start-up and no-load conditions, RAECM is able to reflect the undulation of the filtered on-site measured rotational speed and water pressure owing to the influence of hydraulic factors such as the actual pipeline layout, the elasticity of wall, the inertia of water, and the friction loss. Furthermore, there is a small static deviation between the simulation and the on-site measurement. It might be a consequence of the following aspects: (1) the difference between model test data and the real unit characteristic curve; (2) the imprecise parameters of the pipeline system such as pipe roughness or friction coefficient. In Section 5.1, RAECM reflects the real transient process better than MOC under zero guide vane opening conditions owing to the implicit RadauII A method.
Due to the existence of the "S" area, the stability of the controller is vital for the power quality of PSU under a no-load condition with a low water head. The APFPID controller shows superior control performance with less overshoot and more rapid adjustment than the PID and FPID controllers. Moreover, the variation factor and softening control signal, which take future error trends in the APFPID controller into consideration, can effectively eliminate the oscillation of torque and keep the system stable. Consequently, the APFPID enables the PSUGS to connect to the power grid rapidly and smoothly. In Section 5.2.2, the APFPID controller shows better robustness than the PID controller under low water head conditions.
In the future, RAECM can provide a novel model reference for the optimization of control parameters and the guide vane opening closing laws owing to the high computational efficiency and accuracy. Furthermore, since RAECM is easy to understand and implement by non-professionals, it can also be extended in the power grid to research the interaction between PSUGS and the power system.

Conclusions
A real-time accurate model and its APFPID controller for PSU via error compensation are studied in this paper. The modeling of the pipeline system takes full account of the pipe layout, inertia, elasticity, and friction loss. The pretreatment, which consists of improved Suter transformation and a BP neural network, and online simulation featuring two iterative loops, are proposed to improve the solving accuracy of the pump-turbine. The modified formula of the water head for compensating error is derived with variable spatial discretization to improve the precision of the real-time simulation. Next different solving algorithms are analyzed for RAECM. Then, an APFPID controller based on RAECM is proposed to improve the quality of the transient process under a no-load condition with a low water head. Finally, simulation tests under various operating conditions are utilized to demonstrate the validity and reliability of RAECM compared with MOC, onsite measurements, and the linear model. The excellent control performance of APFPID is also verified under a no-load condition as compared with the PID and FPID controllers. The major conclusions are summarized as follows: (1) RAECM can not only meet the real-time requirement of online simulation but also ensure accuracy and stability under various operating conditions. (2) RAECM with the implicit RadauIIA method is more suitable for simulating a transient process than MOC, especially during the operation period under zero guide vane opening of the extreme conditions. (3) APFPID controller based on RAECM can effectively restrain the oscillation of rotational speed of PSU under a no-load condition with a low water head. (4) RAECM has been applied in semi-physical simulation. It is useful in the commissioning of a power plant because it enables quick testing to reduce the time and cost for tuning system parameters and control settings.