Inverse Optimal Control Using Metaheuristics of Hydropower Plant Model via Forecasting Based on the Feature Engineering

Optimal operation of hydropower plants (HP) is a crucial task for the control of several variables involved in the power generation process, including hydraulic level and power generation rate. In general, there are three main problems that an optimal operation approach must address: (i) maintaining a hydraulic head level which satisfies the energy demand at a given time, (ii) regulating operation to match with certain established conditions, even in the presence of system’s parametric variations, and (iii) managing external disturbances at the system’s input. To address these problems, in this paper we propose an approach for optimal hydraulic level tracking based on an Inverse Optimal Controller (IOC), devised with the purpose of regulating power generation rates on a specific HP infrastructure. The Closed–Loop System (CLS) has been simulated using data collected from the HP through a whole year of operation as a tracking reference. Furthermore, to combat parametric variations, an accumulative action is incorporated into the control scheme. In addition, a Recurrent Neural Network (RNN) based on Feature Engineering (FE) techniques has been implemented to aid the system in the prediction and management of external perturbations. Besides, a landslide is simulated, causing the system’s response to show a deviation in reference tracking, which is corrected through the control action. Afterward, the RNN is including of the aforementioned system, where the trajectories tracking deviation is not perceptible, at the hand of, a better response with respect to use a single scheme. The results show the robustness of the proposed control scheme despite climatic variations and landslides in the reservoir operation process. This proposed combined scheme shows good performance in presence of parametric variations and external perturbations.


Introduction
Hydropower plants management is a problem which complexity is related to the increase of the population on the vicinities of these power generation systems. Taking decisions related to dam volumes and their reservoir of destination poses a real challenge, especially when there are inhabited areas that could be at risk of flooding [1]. Some of the factors involved in dam volumes management include some climatic phenomena, particularly those which may cause landslides and, as a result, changes in storage capacity [2]. Such perturbations in storage capacity are problematic because these cause the system to be unable to satisfy the required energy demand [3], thus, it is necessary to implement techniques to handle said external perturbations and to manage the hydraulic level without hindering power generation or the network's energy supply.
To tackle this problem, it is necessary to model the dynamics hydropower plant [4][5][6][7], in order to calculate the power generated by the system [8]. Recently, several works related to the power regulation problem in hydropower plants have been published; in [9], for example, an application of a robust controller to regulate the speed of the governor was proposed. Following this trend, in [10] a PID controller that allows one to maximize power generation was developed. Then, in [11], an application for a controller based on passivity was presented.
Due to the complexity involved in the development of an integrated control model which satisfies energy demand while also maintaining the desired hydraulic level, in [5] a new control system is designed by decomposing the dynamic system into subsystems, allowing it to tackle several problematics, including the management of the hydraulic head. To reach a desirable performance in the CLS, in [12,13] the implementation of metaheuristic algorithms were proposed in order to optimize the controller parameters.
Another trait as important as the management of the hydraulic head is trajectory tracking, in which a compromise between ensuring the desired output and the discharge evacuation through landfills are sought after. Several works exist within the literature that are aimed to tackle this problem, especially those that describe the IOC technique as an alternative for referenced tracking, in which a passive system is constructed through the feedback of the output, ensuring its stability, and optimizing its response [14][15][16]. However, despite its good performance, this technique is not suitable in the presence of parametric variations within the system [17]. Furthermore, in [18] an integrator is added to the control scheme to gain robustness to uncertainties. Another challenge when implementing IOCs in the modeling of hydropower systems involves the management of external perturbations modifying the input. With the purpose of solving this problem, in [19] a state observer is proposed to handle these perturbations. In [20] an application of both an adaptable fuzzy controller and IOC is proposed with great success in terms of reference tracking in the presence of external dynamics. Furthermore, in [21], a combination of sliding mode controllers and neural IOC was presented, resulting in a robust scheme even in the face of parametric variations and external perturbations, although modeling of hydropower systems is not implemented. While the previously mentioned works have shown good results in their own right, these drive away from the necessity of ensuring optimal trajectory tracking in hydraulic systems, as identifying external perturbations are required. To tackle this problem, many researchers have proposed the use of artificial intelligence algorithms for prediction. A first approximation is presented in [3], where a chaotic model was described for flow simulation subjected to environmental variations. A better approach is presented in [22], where a Nonlinear Auto-regressive Neural Network (NN) with exogenous input was applied with good results, while in [23] predictions on the flow tendency and its remainder was also implemented to compensate the resulting output, thus reducing prediction error. Due to the previous, NN are attractive for predictions, particularly for a task known as neural forecasting [24]. To allow efficient neural forecasting, Feature Engineering (FE) is a crucial step in the machine learning pipeline [25]. However, there are some issues involved in this methodology, such as loss of information and the presence of noise (spurious data). In [26] several data input approaches based on Nearest Neighbor criterion where presented. Then, in [27,28] some techniques based on averaged noise filter for time series are presented. In [29] some criterions aimed to measure statistics anomalies are presented; however, the use of these techniques requires measuring distances and also need the implementation of several data transformations to normalize data [30] and correct the bias related to historical data [31].

Contribution and Organization
In this research, a robust IOC scheme with an RNN is applied to the HP system. The control is adjusted through the Particle Swarm Optimization (PSO) algorithm [32], and the data for the neuronal predictor is worked with FE. The objective of this work is to track a reference to hydraulic head and mechanical power with an unknown flow rate depending on the weather. The main contribution of this paper is to incorporate an accumulative gain for Inverse Optimal algorithm to attenuate the parametric variations and include a NN to forecasting the external disturbance. Also, to train the NN, the historical data is pre-processed using FE methods.
To achieve the objective presented, it is necessary to obtain a discrete-time model of a reservoir, parsed in Section 2. Then, the Inverse Optimal control with an accumulative action is described in Section 3. Subsequently, the NN structure is deduced and the FE methods are applied to perform the forecasting in Section 4. The general scheme is simulated in Section 5. After that, we present the discussion, results, and their interpretation in Section 6. Finally, the paper is concluded in Section 7.

Dynamical Model
Consider the balance flow equation, that describes the variation rate of hydraulic head [4][5][6][7] where, Q in is the unknown input flow. Q out = q s + ∑ n i=1 Q i , with q s as the flow out of spillway and, Q i as the flow in the duct i, for i = {1, 2, . . . , n}. Besides, A is the transverse area of the reservoir and H is the hydraulic head level.
Then, an expression to represents the dynamical contribution for flow rate, is obtained from the second Newton law with the sum of interactions ∑ F = A c Hρg + A c L i ρg h b L i − A c h i ρg, this describes the contribution of the reservoir, the gravity interaction from the duct i, and the opposite force by the pressure in the turbine i, respectively. The mass in each duct is m = A c L i ρ, where A c is the transverse section of each duct, L i is the longitude of conduct, ρ and g are the water density and the acceleration by gravitational interaction, respectively. Notice that, the model contemplates that the reservoir must always maintain a minimum level of water. This level has to be greater than the diameter of the cross-section to avoid variations in mass from (2). Whereas,v i =Q i A c correspond to the derivative of velocity water. Then the flow dynamics is depicted as followṡ Consider also, the energy balance proposed by the Bernoulli equation where, ∆P = P i − P 0 is the pressure difference, with P i as the pressure in the turbine i and P 0 the atmospheric pressure. ∆H = H + h b − h i is the difference between the hydraulic head and the turbine level, with h b as the level from the turbine under the reservoir. The Figure 1 shows the considered architecture. Also, ∆v 2 = v 2 0 − v 2 i is the difference between the square of the top and turbine flow velocities. The velocity at the top is assuming v 0 ≈ 0.
The pressure in each turbine is represented as a water column of h i length, this is P i = P 0 + h i ρg. Then, from (4), we can see that To apply the Inverse Optimal Control algorithm, the system (1)-(3) is discretized where, the head at the turbine i is being G i (k) the i-th gate opening percentage. Finally, the turbines mechanical power generated is computed as follows [8] Once the model is presented, first, the algorithm to control the plant is described, and then to this application case is designed.

Inverse Optimal Control
Inverse optimal control is a technique, which consists in designing a stabilizing feedback control law, which is based on an a priori known Control Lyapunov Function (CLF), that ensures the optimization of a cost functional [15].
Considering the affine in the control input of a discrete-time nonlinear system where x 0 = x(0) and x(k) ∈ R n is the state of the system at time k, the output y(k) ∈ R m . Also, the mappings h : R n → R m , and j : R n → R m×m are smooth. Under assumption h(0) = 0.
To establish the target system, lets consider the following definition Definition 1 (Passivity). The dynamic system (10) is said to be passive for a supply rate w(y(k), u(k)) = y T (k)u(k) if there exists a non-negative function S(x(k)) with S(0) = 0, defined as the storage function, such that for all u(k) which is named the passivation inequality [14].
For passive systems, asymptotic stability can be achieved using the output feedback setting the storage function S(x(k)) > 0 as a Lyapunov function V (x(k)).
On the other hand, a purpose of this control scheme is to make (10) to a passivated system through the feedback input, this system is defined as follows with v(k) ∈ R m as the new input and α is a smooth function of the suitable dimension. Such that, the CLS (10), (12) is described by The Inverse Optimal Control is formulated as follows: consider the passivity condition where P = P T > 0. And the control input as in (12), where the passifying law is where I m is the identity matrix, and If there exists P such that then, the system (10) with input (12) is globally asymptotically stable at the equilibrium point x(k) = 0. Moreover, is inverse optimal in the sense that the input control minimizes the meaningful functional given by [16] J ) − x(k) T Px(k)) ≤ 0 and the optimal function value is J * (x(0)) = V (x(0)).
For this application, two inverse optimal control schemes are applied. First, to tracking a level reference from the hydraulic head using the spillway as a control input, and then to tracking a mechanical power using the gates to manipulate the flow rate for each turbine.
To apply the desired algorithm, it is necessary to know the discrete model on its totality (6)-(7); in particular, for the hydraulic level, it is necessary to know the flow rate Q in (k) occurring at every day. With this in mind, a NN is proposed to predict the flow rate. The proposed general scheme is shown in Figure 2. Once the flow rate occurring during the day is predicted, the IOC is applied to the continuous plant model by knowing the level and flow rate references.

Inverse optimal control Reservoir
Predictor Figure 2. Control-forecasting scheme to reservoir.
It is supposed that A, A c , L i , and h b are all known, although this may have slight measurement variations with regard to their real dimensions, this under the supposition that, in case of a landslide, some of their magnitudes may become altered. To compensate for these deviations, an error cumulative action is proposed for each of the outputs that are required to be controlled.
where H r (k) and Q i,r (k) are the hydraulic head and flow rate references, respectively. Also, i H (k) is the level error accumulation and i Q (k) is the flow rate error accumulation. These magnitudes correspond to the discretization of integral action for continuous systems [18]. Then, for the hydraulic head subsystem (6)- (16), the estimation of input flow rate Q in (k), and the reference H r (k) are assume knowing. To control q s , this estimation is adding, then, the control law results wherej the matrix P H is found through PSO algorithm. Finally, as a new control input to the passivated subsystem is choose since the subsystem output is defined as a linear combination of the tracking error and the error accumulation.
In the case of the mechanical power, by knowing the reference P r (k), we calculate Q i,r (k) from H(k) using (9) by imposing the same flow circulation in all turbines, this is Q i,r = Q i+1,r for i = {1, 2, . . . , n − 1}. Thus, with a known Q i,r , we use h i (k) as a pseudocontroller to calculate the opening of each gate G i using (8). Thus, the Inverse Optimal Controller used for tracking the reference in subsystem (7)- (17) is wherej Matrix P Q is found using the PSO algorithm [32]. Also, the new input for the control of this passive subsystem is Thus, from (8), using head at the turbine as in (21) we get the control law for the opening which represent the opening rate G i ∈ (0, 1], on each gate so that the pressure on each turbine i correspond to a water column h i in meters.
To find the matrix P H and P Q , we used the PSO algorithm describe in follows.
The algorithm describes the social behavior of birds and fish. It relies mostly on the basic principles of self-organization. The algorithm starts with a population of random position, the control scheme in Figure 2 is evaluated under a minimization criterion J, the best positions are selected from the proposed criterion. Then, the positions and velocities of the particles are recalculated so that the objective function decreases, this process is repeated until a stopping criterion is met. The flowchart is depicted in Figure 3.
To update the population, the follow equations are presented where V i (t) is the velocity of particle i at time t, K i (t) is the position of particle. The first element in velocity equation is the inertia term, next the cognitive component, and finally, the social component. Besides, w(t) is the inertia coefficient, c 1 and c 2 are the personal and social acceleration coefficient, respectively. r 1 and r 2 are two uniformly aleatory numbers between (0, 1]. L i (t) is the personal best position of particle i, and B(t) is the social best position. Additionally, the inertia coefficient is damping in each iteration by.
Then, the matrix P H and P Q used in (19), (20), (22), and (23) are constructed from B at final iteration as follows . From that, the dimension of each particle is six.

Forecasting
The inverse optimal control algorithm by itself is not robust [17], since it is limited to the fact that, through an inverse model, the system behaves passively, to later stabilize it. For this reason, an integral gain is added to each output to attenuate the parametric variations of the system, however, in the presence of external disturbance, such as Q in , the results may not be as expected.
To reduce the effect of flow velocity at the input, a NN is designed. this NN is based on the nature of the historical data. To obtain a better representation of the signal, a series of treatments are applied to reduce the possible effect of noise on sensors, human error, or data loss. Treatment is as follows: 1. Missing values imputation using Nearest Neighbor method [26]. [27,28]. 3. Yeo-Johnson transformation to minimize variance in prediction [31]. 4. Identification and imputation of atypical data by the quantile method [29]. 5. Offset and scaling [30].

Discrete low-pass filtering
When obtaining a prediction, the inverse of scaling and offset must be applied to the result.
Once we have our pre-processed data is obtained, we construct the NN by taking the following criteria into account: 1. The signal auto-correlation.

Hypothesis 1.
The auto-correlation indicates the number of sufficient inputs of previous data that have less relation with each other. In this way it does not fall into information redundancies.
2. Order estimation of the system that generates the signal.

Hypothesis 2.
If the signal is a solution of a differential equation of order n, then, n − 1 feedback outputs are sufficient.

Signal frequency.
Hypothesis 3. The number of hidden layers in the NN is directly proportional to the frequencies of the signal to be represented.
4. Finally, the number of data to be predicted is considered, which is the number of outputs of the mentioned network.
These hypotheses are proposed and applied in this work but we plan to develop them in-depth as the main contribution in a future work.
To calculate the auto-correlation, we define a set of elements that constitute the signal as s(a, b) = {Q in (a), Q in (a + 1), . . . , Q in (b)}, with 0 < a < b. Then, the auto-correlation using the x previous elements is calculated as follows: where, µ x and σ x are the mean and standard deviation for s(k − x, k), whereas µ and σ stand for the mean and standard deviation for s(k − x, k + 1).

Tracking Control
Let's take the Inverse Optimal Control scheme shown in Figure 2. As reference, we use the hydraulic head and mechanical power described by the following equation H r = 60 + 5 cos(0.05t) P m,r = 0.8g 60 + 5 cos(0.05t) 1300 + 100 sin(0.2t) sin(t/365) . If the hydraulic system has three turbines, and the control reference for each flow is known, for an efficiency of η = 0.8, we have that: Assuming a trained NN, the control gains are found through the PSO algorithm. The stopping criteria is to reach 100 iterations. The objective function is where T 1,Q and T 2,Q are the transitory and stationary time for the flow rate. T 1,H and T 2,H are the times analog to flow rate for hydraulic head. J Q and J H are the objective functions associated to tracking error from flow rate, and hydraulic head, respectively. The corresponds objective function are computed as being, w 1,Q and w 2,Q are constants to weigh the transient and steady state for Q i . On the other hand, the objective function of the hydraulic head is with w 1,H and w 2,H as constants to weigh the transient and steady state for H. In other matters, the constants used in (27) and (28) are described in Table 1. Furthermore, on the PSO algorithm, we use the parameters described in Table 2. Besides, the range search is between [−1, 5000]. The algorithm is simulated using the code provided by the Yarpiz team [34].

Inverse Optimal Control including Parametric Variations
Further, at May-2013 a landslide is simulated, this causes an increase in the hydraulic head of 10 m, a decrease of the reservoir area, and obstruction in duct area of 50%. It is expected that these parametric changes are attenuated by the cumulative dynamic proposed in this work. The simulated system includes a discrete model of the hydropower plant (6)-(7) and the Inverse Optimal Controller (18)-(24) with cumulative action (16)-(17).

Inverse Optimal Control including External Disturbances
To eliminate the undesired effect of increasing the hydraulic head, knowing the flow that serves as the input of the hydropower system is required. Thus, it is needed to construct a neural predictor to estimate the flow at each day k, to do this, we first process the data needed to train such predictor. The historical data for Q in is shown in the Figure 4 and its attributes are depicted in Table 3. From this table, it may be noticed that the minimum value of the flow is rate negative; however, as the flow rate must be positive, these negative values are considered outliers.  The results of inputting the missing values and applied the discrete low-pass filter, is displayed in Figure 5, after, the data is studied from a frequency point of view, then, a Yeo-Johnson transform is applied to discriminate the outliers values, this is considered all values 1.5(q 3 − q 1 ) times under q 1 or upper q 3 , where q 1 and q 3 are the first and third quartile of the data filtered, respectively. A comparative histogram is shown in Figure 6. The parameter to transform the data is calculated from the maximum likelihood method, this is λ = 0.12217. Once the obtained data is clean, based on the attributes shown in Table 4, the negative values suggest an offset on the sensor reading of at least 0.7090. This value must be found by testing the model recursively on the plant. A first compensation value of 0.8000 is proposed. The input for the transformed data is scaled so that its response values are −1 ≤ x ≤ 1 with mean zero, and they are added to time in the form of a sine/cosine signal with a period of one year to the NN input. These inputs are shown in Figure 7. For the design of the NN, we find the auto-correlation of the clean data. In Figure 8 it is shown that with x = 182 the first minimum is reached, thus this value is chosen as the number of inputs for the previous days. To know the number of output elements yr that are feedback to the NN, we start from the supposition that this signal may be generated by knowing its actual value, as well as its derivative, this supposing a 2nd order system; this means that a delay of the output is enough to feedback, y r = 1. Furthermore, it is supposed that, due to its periodicity, the NN only requires a hidden layer to be applied. The number of neurons before the output layer is directly proportional to the desired resolution, thus, it is proposed to use 30 of them. The number of neurons on the input layer is given by 5(x + yr + dt), where dt stands for the number of elements considered to model time. The scheme of the applied NN is shown in Figure 9. The historical data taken to train the NN contemplates the hydraulic head monitored through January-2006 to April-2011. On the other hand, data from April-2011 to December-2014 is utilized to test the prediction model.

Model: Technical Considerations
Due to the control robustness, a non-linear model is described simplifying the losses due to evaporation, friction, leaks, etc. These are absorbed by the cumulative action of the IOC.
The model contemplates that the reservoir must always maintain a minimum level of water. This level has to be greater than the diameter of the cross-section to avoid variations in mass from (2).
The power generated by the turbines is calculated as a stationary state. This is assuming an adequate speed regulation in the governor.
The model includes a spillway as a control input which is used as an alternative to releasing water when the levels of the dam exceed the reference. This work does not contemplate policies of release by spillways that limit the outflow. A general multipurpose study is necessary for the proper management of the water diverted by the system. This analysis should consider the need to satisfy specific demands, such as agriculture, navigation, human accents, among others. To know the destination of the water released by the spillway, as well as the limits of the river to avoid risks and losses in the aforementioned aspects.

Tracking Control
At the end of the metaheuristic algorithm, the best cost was 0.0128 with position B = [6.4409, −1.0000, 15.7979, 50,000.00, 21.654, 2835.136]. The Figure 10 shows a comparison between hydraulic head reference and Inverse optimal control output without forecasting (dotted line), Figure 10a. Also, the Figure 10b shows a comparison between the flow rate reference (discontinuous line) and the output without forecasting (dotted line). From this, it is observed that there is a period where the hydraulic head remains below the reference. This is because the permitted limit is greater than that provided by the input flow, which suggests a dry season.

Inverse Optimal Control including Parametric Variations
By simulating the controller without prediction data, we get a root mean squared error of 0.0145. Thanks to the integral action it is possible to handle perturbations Q in , however, the hydraulic head is above the permitted limit, which could be undesirable. The landslide is marked with a star in Figure 10. Figure 10a shows a comparison between hydraulic head references and IOC output with accumulative action (dotted line). Furthermore, Figure 10b shows a comparison between the flow rate reference (dashed line) and the output (dotted line). It is observed that the controller manages to adjust the response of the system to the desired reference. A prediction block is included to improve this response.

Inverse Optimal Control including External Disturbances
The response of the applied NN is shown in Figure 11. Data suggest that the filtering criterion is very conservative. There is a presence of values that provoke an unlikely derivative in the signal (this by considering a Bayesian criterion); however, these values are not eliminated with the noise filtering. These results motivate the use of a filter that limits the derivative, which can also update itself as time progress. We simulated the operation process of the dam by using the flow rate historical data corresponding to 365 days of operation, starting from July 2012. The resulting prediction outputted by the NN (Figure 9) is used to calculate the inverse Yeo-Johnson function, giving as a result what is shown in Figure 12. However, after discarding and replacing the spurious data with the Yeo-Johnson data transformation ( Figure 6) it is possible that the resulting signal has a distribution when said transformation is applied once more. In future works, it is planned to explore the idea of a recursive data transformation-elimination process to minimize spurious data. With the inclusion of the NN for prediction we achieve a root mean squared error of 0.0143. The Figure 10a shows a comparison between hydraulic head reference and Inverse Optimal Control output with forecasting (continuous line), Figure 10a. Also, the Figure 10b shows a comparison between the flow rate reference (dashed line) and the IOC output using forecasting (continuous line). The landslide is also simulated in May-2013, causing the same variations of 50%. The control response is depicted in Figure 13, the evolution of q s (k) is in Figure 13a and the gate of a turbine is shown in Figure 13b. The power generates is shown in Figure 14, where it is compared with the reference signal (dashed line). In here, there are two important observations 1. The output in the hydraulic head should always be below the reference since it is considered that the controller input q s is always positive; this is, the flow is always unidirectional, Figure 10. 2. Besides handling parametric variations (such as the case of landslides), the proposed Integral Action Inverse Optimal Controller scheme is also able to properly manage external perturbations; in fact, it is observed that the output of the controller does not differ significantly from the reference, Figure 10.

Tracking Control
The implementation of the Inverse Optimal Control technique for tracking along with the application of the PSO algorithm, the model manages to minimize the tracking error, thus demonstrating it to be a good alternative for adjusting to a time-varying reference. It is been observed that, when modeling a system with abundant flows, a negligible margin of error of up to 1.7% is observed. Following the recommendation proposed in the discussion, it is suggested to establish upper/lower reference limits that allow one to keep the head level below the reference, this is a safety measure.

Inverse Optimal Control including Parametric Variations
In the event of an unforeseen event (such as a landslide), it is been shown that the Inverse Optimal Control scheme is capable of following the reference with a minimum deviation, thus ensuring that the systems level is maintained close to the reference while also warranting a supply which satisfies the current energy demand despite a reduction in its storage capacity.

Inverse Optimal Control including External Disturbances
By applying Feature Engineering techniques to the system's historical data and subsequently feeding the neural network with this information, an excellent approximation is achieved in the training stage and a minimum error in the predictions.