Abstract

Increasing renewable energy usage puts extra pressure on decision-making in river hydropower systems. Decision support tools are used for near-future forecasting of the water available. Model-driven forecasting used for river state estimation often provides bad results due to numerous uncertainties. False inflows and poor initialization are some of the uncertainty sources. To overcome this, standard data assimilation (DA) techniques (e.g., ensemble Kalman filter) are used, which are not always applicable in real systems. This paper presents further insight into the novel, tailor-made model update algorithm based on control theory. According to water-level measurements over the system, the model is controlled and continuously updated using proportional–integrative–derivative (PID) controller(s). Implementation of the PID controllers requires the controllers’ parameters estimation (tuning). This research deals with this task by presenting sequential, multi-metric procedure, applicable for controllers’ initial tuning. The proposed tuning method is tested on the Iron Gate hydropower system in Serbia, showing satisfying results.

HIGHLIGHTS

  • Uncertainty of the boundary and initial conditions affects model-driven forecasting.

  • Data Assimilation is used to overcome these problems.

  • Research presents potential of using novel, tailor-made, PID controllers-based data assimilation method for river hydraulic models update.

  • Method could be used as a decision-support tool for hydropower systems control.

  • Sequential, multi-metric tuning procedure has been introduced.

INTRODUCTION

The global tendency for reduction of greenhouse gas emissions switches focus from fossil to renewable energy sources. This tendency puts hydropower systems, as the most important manageable renewable source, in the first place, considering its availability and cost-effectiveness (Eurelectric 2011). Hydropower production depends on numerous dynamically changing factors, like extreme hydrological events (floods and droughts), environmental constraints, price fluctuation on energy markets, and conflicting goals of the stakeholders in charge. To use the water optimally, engineers and operators in charge of hydropower systems require different decision support tools (Ehteram et al. 2018). In many cases, these tools must provide forecasts of the hydrological data (water levels and river flows) on a daily basis. Hydrological forecasting is a challenging task, addressed by numerous researchers (e.g., Wu & Chau 2013; Taormina & Chau 2015; Fu et al. 2020). For the purpose of forecasting, monitoring systems and physically based models are used in model-driven forecasting procedures.

The quality of model-driven forecast directly depends on boundary and initial conditions. However, they are often affected by numerous sources of uncertainty (Bozzi et al. 2015; Ocio et al. 2017). These problems eventually lead to poor forecast results and inadequate hydropower system control operations. To overcome these problems, data assimilation (DA) methods are widely used for improving hydrological-hydraulic forecasting (Vrugt et al. 2006).

Various DA methods are applied in hydrological-hydraulic modeling. One of the widely used methods is ensemble Kalman filter (EnKF) (Evensen 1994, 2003). This method is often used for improved river flood forecasting (Madsen et al. 2003; Neal et al. 2007; Li et al. 2014; Barthélémy et al. 2017) but also for forecasting flows and overflows in urban drainage systems (Lund et al. 2019). In many hydrological-hydraulic applications of DA, measured water levels from monitoring systems are used to update the model (Romanowicz et al. 2006; Hostache et al. 2010; Jean-Baptiste et al. 2011; Rakovec et al. 2012). Some authors have successfully used streamflow observations too (Thirel et al. 2010; Dumedah & Coulibaly 2014; Randrianasolo et al. 2014; Sun et al. 2015). In more recent researches, there has been a tendency to use remotely sensed hydrological data, such as satellite images (Matgen et al. 2010; Giustarini et al. 2011; Yoon et al. 2012; García-Pintado et al. 2013; Andreadis & Schumann 2014; Garambois et al. 2020) or even crowdsourced data (Mazzoleni et al. 2017, 2018; Annis & Nardi 2019).

All these studies show high applicability of the standard DA methods, such as EnKF, in hydrological-hydraulic modeling. On the other hand, in most real-world applications, EnKF, and similar DA methods, struggle to perform within a reasonable time frame (Madsen & Skotner 2005), which makes them not favored by engineers and operators in charge of water systems. Additionally, mathematical apparatus, which is relatively complex, means these methods avoided in real-world applications. Owing to the necessity for improved water resources management and control, many decision support tools are required, especially for improved short-term forecasting. This enables improved (near) real-time operations in various water resources (Schwanenberg et al. 2015; Ahmad & Hossain 2019). This also makes DA an important decision support tool. Therefore, problems with real-world applications of standard DA methods are under the spotlight, creating a challenging task. To overcome these problems, researchers are searching for simplified and fast DA methods. Many of them are trying to develop tailor-made assimilation techniques, suitable for solving some specific problems. Madsen & Skotner (2005) developed a cost-effective assimilation method for river models. This method uses the simplified version of Kalman filtering, which certainly provides reduction of computational time, but requires the predefinition of the parameters used in a pre-processing phase. All these applications of ensemble-based DA methods or its derivatives use the direct approach for model updating. In hydraulic models, this means that water levels (stages) are directly corrected according to measured values. Sometimes, this procedure can cause ‘shocks’ in the model, disrupting the continuity equation (mass balance). Therefore, some researchers are trying to use a different and indirect model updating procedure. For example, Fava et al. (2020) presented a fast DA methodology for improved flood prediction by hydrological models in urbanized areas, using the rainfall input corrections instead of direct water level update. The methodology for indirect water level update in urban drainage systems, presented by Hansen et al. (2014), uses an approach of adding/subtracting correction flow to the system in the junctions where observations are available. This approach, with all details in the Methodology, has narrowed the scope of its application.

The general approach of the indirect model update procedure (Hansen et al. 2014) was the basis for the introduction of DA methodology for (near) real-time river model update where correction flows are calculated using control theory (Rosić et al. 2017a, 2017b; Milašinović et al. 2019). Development of this novel, control theory-based DA method and EnKF benchmark test was presented and verified by Milašinović et al. (2020). The study has proved the applicability of proportional–integrative–derivative (PID) controllers as a DA tool. Due to its potential as a time-efficient method, compared to the widely used EnKF method, it can be used as decision support tool for near real-time control of river hydropower systems. However, the study underlined the need of further investigation in the field of controllers’ tuning.

Assigning the values to the controllers’ parameters (tuning of the controllers) is one of the most important steps in any application of the PID controllers. Previous research used trial-and-error tuning rather than any systematical tuning procedure. For further and more general applications of this novel DA method, proper controllers’ tuning procedure has to be proposed. Therefore, this research addresses tuning of the controllers by presenting sequential and heuristic multi-metric tuning procedure. The goal is not to determine a global optimal solution on the controllers’ parameters, but to present a simple, engineering procedure to estimate the initial values and possible range of the controllers’ parameters preserving the model's stability. The presented procedure shows satisfying results, and thus can be considered for initial tuning of the controllers.

MATERIALS AND METHODS

Overview of the methodology

Standard DA procedure (e.g., EnKF) corrects the state variable (e.g., water level) based on observation data, a model's sensitivity on state variables and compared model uncertainty and observation uncertainty. When model uncertainty is bigger than observation uncertainty, DA procedure puts more trust in observations, correcting the model output towards observed data, and vice versa. This approach is derived using the assumption that there is no dominant source of uncertainty. On the other hand, in many engineering problems, some simplified approaches can, or must be used.

Control theory-based DA is developed using the assumption that system state observations (water levels) are with much lower uncertainty than models and can be used for direct, deterministic model update (Hansen et al. 2014). Therefore, comparison of the model and observation uncertainty is skipped, which allows significant speed-up of the assimilation process (Milašinović et al. 2020). The second used assumption is that dominant sources of the model uncertainty are boundary and initial conditions, while model parameters are properly determined during model calibration. In most cases, inflow data are not directly measured, have high uncertainty, especially during flood events, and are with poor time resolution. Additionally, there is lack of the input data (e.g., rainfall) for large water systems due to inability to implement monitoring stations, with good spatial resolution, across an entire catchment area. Therefore, initial and boundary conditions are inadequate and could be blamed for discrepancy between calculated and observed water levels. Eventually, these problems lead to a bad estimation of reservoir state and hydropower system operation rules.

Considering these two assumptions, a novel control theory-based DA with simplified model update procedure is developed. This DA procedure treats the problem (water level discrepancy) as the lack or excess of water amount in the model. Hence, the model update procedure is operated by adding or subtracting the water from the model at the locations where observations are available. This reduces the water level discrepancies, improving the initial conditions for future forecasting. The process of adding/subtracting water from the model is controlled by PID controllers (hence it is named control theory-based DA). Figure 1 represents the general model update procedure in control theory-based DA.

Figure 1

Algorithm for control theory-based DA procedure.

Figure 1

Algorithm for control theory-based DA procedure.

PID controllers as data assimilation tool for hydraulic model

Model-driven forecasting requires physically based models. In this paper, diffusion wave hydraulic model is used, as proposed by Milašinović et al. (2019, 2020). This model is based on one-dimensional (1D) Saint-Venant's Equations (1) and (2). Diffusion wave model for 1D open channels (DiffW1D) uses a simplified dynamic equation (Equations (2)) derived from the full Saint-Venant's equations by neglecting convective acceleration (second term in Equation (2)):
formula
(1)
formula
(2)
Here, x is space coordinate, t time, A cross-sectional area, Q discharge, q lateral inflow, Z local water surface elevation (water level), g acceleration due to gravity, R hydraulic radius calculated as the ratio between cross-sectional area and wetted perimeter, β velocity distribution coefficient, and n is Manning's roughness, calculated according to Costabile & Macchione (2012). Model domain discretization is presented in Figure 2(a). Numerical DiffW1D model uses a staggered numerical scheme where water levels and flows are calculated in alternating cross sections, as presented in Figure 2(b). Numerical model of the diffusion wave is given by Equations (3) and (4):
formula
(3)
formula
(4)
where index i represents spatial location (i-th cross section), B is the top width of the cross section (water surface width), t in superscript denotes current time, Δx spatial resolution, and Δt is temporal resolution.
Figure 2

(a) Model domain discretization and (b) numerical scheme using DiffW1D model (Milašinović et al. 2020).

Figure 2

(a) Model domain discretization and (b) numerical scheme using DiffW1D model (Milašinović et al. 2020).

Procedure for adding/subtracting water to/from the model is based on adding/subtracting correction flows, in a form of lateral inflow, at the location where observations are available (assimilation point). Correction flows are calculated based on the difference between observed water levels and simulated water levels (error). Conversion of the errors to correction flows is conducted using the PID controllers. Correction flows are implemented as (fictive) lateral inflow elements in the continuity equation (Equation (3)) of the hydraulic model. A new form of continuity equation is presented by Equation (5).
formula
(5)
where Qcorr represents correction flow that reduces the water level discrepancy between observations and model. It is calculated using the PID controller's theory.

PID controller is a control loop feedback mechanism, where input in the next step is a function of the previous output (Astrom & Hagglund 1995). This mechanism is often used for real-time-control (RTC) of different process (e.g., RTC of hydraulic structures in urban drainage systems (Schütze et al. 2004), RTC of irrigation canals (Bonet et al. 2019), or RTC in water resources (Schwanenberg et al. 2015)). In this paper, the mechanism is used for on the fly hydraulic model update.

The PID controller's input is named as error e(t), which is calculated as the difference between current value of the process variable (simulated water level − Zmodel) and the setpoint of the variable (observed water level − Zobs), as stated in Equation (6). In many real-world problems the simulation time step Δt is much shorter than the observation time step Δtobs. Therefore, in the period between two existing observations, the ‘observed’ water level Z*obs is calculated using linear interpolation. It can be assumed that the accuracy of interpolated ‘observed’ level decreases as the time interval from the last observation increases. To reduce the influence of the interpolated water levels Z*obs, multiplier C, named as attenuation factor (Figure 3), is used. In Equation (7), tprev_obs is time of the previously available observation, t is current simulation time, tobs is observation time, and Δtobs is observation time step. Attenuation factor in Equation (7), will gradually turn off the PID controller (Figure 3) in periods between two measurements. This means that, as the model progress forward in time, in periods without measurements, smaller weight is given to the errors calculated using interpolated water levels.

Figure 3

Error attenuation factor (periodical, discontinuous, function) as a measure of uncertainty (minimal value of the attenuation factor depends on specific values of simulation time step Δt and observation time step Δtobs) (Milašinović et al. 2020).

Figure 3

Error attenuation factor (periodical, discontinuous, function) as a measure of uncertainty (minimal value of the attenuation factor depends on specific values of simulation time step Δt and observation time step Δtobs) (Milašinović et al. 2020).

The PID controller tends to reduce error using the control variable. The control variable used to reduce this error is lateral inflow Qcorr(t) presented by Equation (8).
formula
(6)
formula
(7)
formula
(8)

PID parameters are: Kp – proportional gain factor used to multiply the current error value, Ki – integrative gain factor used to add the influence of error history, and Kd – derivative gain factor used to adapt control to current trends in error change. Proportional gain produces an output based only on the current value of the error. High values of Kp cause big variations in controllers’ output that can make the system unstable (extremely big water level oscillations causing the occurrence of physically impossible values of water depths). Low values (towards zero) of Kp avoid problems of an unstable system, but time needed for reaching the setpoint increases and, practically, makes system unable to reach the goal. Therefore, Ki, along with error integral, forms integrative gain, which is used to collect error history and its duration, to minimize them over time. This gain can significantly reduce time needed for reaching the setpoint. In some cases, in highly dynamic systems with rapid changes, the derivative component is included to estimate the error trend (what will be the error in the near future). However, the derivative component is sensitive in systems with high measurement noise and can enhance the controllers’ instability. Hence, proper tuning of the PID controllers’ gain factors must be provided. Finding optimal values of the controllers’ parameters is a challenging task. In common applications of the PID controllers this is often conducted using some of the simple, heuristic methods (e.g., Ziegler & Nichols 1995; Skogestad 2004). However, these tuning methods, systematically presented by Ang et al. (2005), are mostly physical systems-oriented (in other words, they are proven for automated control of the real physical systems). Recommendations in these methods are derived from the specific field of applications and are unreliable when the PID controllers’ theory is used in hydraulic modeling. Additionally, they are single metric or just visually based rather than multi-metric based.

Rather than using simple, heuristic tuning methods, this procedure can be done by optimization algorithms (Ou & Lin 2006; Altinten et al. 2008; Solihin et al. 2008; dos Santos Coelho 2009; Wang & Wei 2009; Chiha et al. 2012). The main problem in optimization algorithm usage (e.g., genetic algorithm) is the difficulty to find the global optimal solution (in reasonable time) when the number of variables increases (multiple controllers). The solution depends on the initial assumption on the variables’ values (Deng et al. 2015). To prevent this problem, appropriate initial values of the controllers’ parameters (Kp, Ki, and Kd) should be estimated. Therefore, simple tuning procedures are necessary to complete this step.

This paper presents sequential, heuristic, multi-metric tuning procedure for initial tuning of the controllers. It is a simple tuning procedure and convenient to start with. Hence, it should not be omitted. Further, fine tuning can be performed using the results of sequential multi-metric procedure as an initial point. Procedure is based on gradual increase of the gain factors of each component independently in three phases. Phase 1 is used to estimate the value of the proportional gain only (P controller) that provides fastest approaching to the prescribed level, with minimal overshooting or oscillations. In Phase 2, integrative gain is coupled (PI controller) with the value of proportional gain factor Kp which produced the best assimilation results according to assimilation quality indicators. In Phase 3, derivative gain is added (PID controller) to the values of Kp and Ki which showed the best results. Kd is increased to analyze if there is any improvement by adding the derivative gain in the controllers. Phases have to be performed in this order, because application of derivative gain as the first or coupling it first with proportional gain (PD controller) is not recommended (suggested by many researchers, e.g., Stanbury et al. (2017)). Additionally, almost 80% of PID controllers’ applications have the derivative part switched off while the other 20% use it coupled with proportional and integrative gains (Ang et al. 2005), because inappropriate tuning of Kd can affect the system's stability. This is a solid reason to use the proposed phases in sequential tuning procedure.

Indicators for DA performance assessment

To assess the effect of the multi-metric tuning procedure, DA performance assessment (DA quality) indicators must be established. The aim of the PID controllers’ tuning is to minimize difference between process value and the setpoint, oscillations in process error and settling time. For the application of PID controllers as a DA tool, four indicators are used in this paper: root mean square error RMSE (to represent difference between process value and the setpoint), amplitude of the process error maxError (to represent the oscillations), assimilation time ratio AssimTRatio (to represent the settling time), and total correction volume CorrVol, which is method-specific indicator, a product of PID controller-based assimilation (it represents the intensity of the controllers’ ‘intervention’ into the hydraulic model).

Root mean square error

Usage of the root mean square error RMSE as DA quality indicator depends on the number of locations where direct model update is conducted (number of assimilation points). In this research, mean value of RMSE indicators for all observation locations is used (Equation (9)).
formula
(9)
where Zj,obs represents sample from the observed water level time series, Zj,sim represents sample from the simulated water level time series, N is the number of time samples in time series, and M is the number of locations used for calculation of RMSE indicator.

The lower bound of this metric is 0, which represents the ideal case showing the excellent assimilation performance.

Amplitude of the process error

During the assimilation process, the model is updated towards the measured values, but big oscillations of simulated water level are present, especially at the beginning of the assimilation process. Therefore, amplitude of water level oscillations (error amplitude) maxError is used to estimate DA quality (Equation (10)). Since there are several observation locations where DA performance has to be estimated, mean value of error amplitudes is used. The lower bound of this metric is 0, and DA performance is better when maxError tends to this value.
formula
(10)
where ei is error time series for the i-th observation location and M is the number of those locations.

Total correction volume

Total correction volume CorrVol represents the total volume of water added or subtracted by controllers at the assimilation points. This indicator is method-specific, applicable only in PID controller-based assimilation. It represents the estimation of controllers’ intervention into the model. When confidence in boundary conditions is high, the total correction volume tends to be zero. In this research, it is evaluated as a sum of total correction volumes added/subtracted from the model at all assimilation points (Equation (11)). This indicator also has the lower boundary set at 0.
formula
(11)
where MA is the number of assimilation points, tsim is the time period of the simulation when assimilation is performed, and Qcorr is the correction flow.

Assimilation time ratio

Assimilation time ratio AssimTRatio [/] represents the ratio of period in which process error exceeds the threshold to the total simulation time tsim. When AssimTRatio ratio has value 0, it shows excellent DA performance (best value), while value 1 shows poor DA performance (worst value), where assimilation fails to correct the model towards measurements. To evaluate the AssimTRatio, error time series has to be transformed into the error duration curve (Figure 4). The threshold used in this paper is set to 0.05 m. As there are several observation locations where assimilation quality has to be estimated, AssimTRatio is calculated as the average of assimilation times for all locations.

Figure 4

Transformation of the error time series into the error duration curve for determining the AssimTRatio indicator.

Figure 4

Transformation of the error time series into the error duration curve for determining the AssimTRatio indicator.

Test case

The presented methodology is developed for the purpose of near real-time control of transboundary hydropower system Iron Gate 1 on the Danube river shared between Serbia and Romania. Everyday operations of the hydropower plant (HPP) require short-term forecasting of water levels (and river flows) in the vicinity of HPP. For improved forecasting (improved initial conditions) DA is required.

The test case of the hydropower system consists of hydropower plant Iron Gate 1 (IG1) and a 170 km-long upstream Danube river section (Figure 5), the lower half of the total reservoir's length. The domain is discretized by 189 cross sections, with the average distance of 900 m between. The river model is developed using the DiffW1D model with one upstream inflow, neglecting all tributaries. Tributaries are neglected due to lack of data, further supporting the assumption that dominant uncertainty source in the model is the inflow data. Six observation locations are available across the river section: Nera (132 km from IG1), Golubac (100 km from IG1), Dobra (74 km from IG1), D. Milanovac (47 km from IG1), Dubova (25.7 km from IG1), and Orsova (10.3 km from IG1). The first five locations are used as assimilation points (PID controllers are set at these locations), while Orsova is used as a validation point. The test case used for the analysis of the proposed DA methodology contains real river geometry.

Figure 5

Iron Gate hydropower system with 170 km-long Danube river section coupled with PID controllers for continuous model update.

Figure 5

Iron Gate hydropower system with 170 km-long Danube river section coupled with PID controllers for continuous model update.

Synthetic test case scenario is used for analysis in this paper. Observed water levels at six locations are generated using the ‘true’ inflow hydrograph (black line in Figure 6). The total inflow volume for a given 7 days is 2.163 × 109m3. Then, the ‘true’ inflow is altered and changed to ‘uncertain’ inflow, or model driving inflow (dashed line with black circles in Figure 6). Difference between the ‘true’ and the ‘uncertain’ inflow is a sum of two sinusoidal functions with the amplitude of 1,000 m3/s and frequencies of π/(10,000 s) and π/(5,000 s), which are arbitrarily selected. The scenario is adopted to emulate extreme conditions where model driving inflow (‘uncertain’ inflow) significantly differs from the ‘true’ inflow (in this study used to emulate ‘true’ water levels). This also represents the case of sampling the inflow data on poor timescale, which results in big discrepancies. It is also used to show the potential of using PID controllers as a DA method in a more dynamic environment (note that ‘true’ inflow is more dynamic than model driving inflow). The altered inflow is set as the upstream boundary, with total volume of 2.268 × 109, with average ‘uncertain’ inflow which is larger than ‘true’ inflow for 104.7 × 106 m3/s. At the downstream boundary, the measured outflow and stage (dashed line and dashed line with square markers in Figure 6) are used. Stage hydrograph is used to generate the true state, while outflow hydrograph is used as model driver. Water level observations are obtained with Δtobs = 1 h, while computation time step Δt is set to 60 s, simulating real-world conditions where computation time step is, in most cases, smaller than observation time step. Total simulation time of the test case scenario is tsim = 7 days.

Figure 6

Boundary conditions: true inflow, uncertain upstream inflow hydrograph (model driving inflow), downstream outflow hydrograph, and downstream stage hydrograph.

Figure 6

Boundary conditions: true inflow, uncertain upstream inflow hydrograph (model driving inflow), downstream outflow hydrograph, and downstream stage hydrograph.

Stage (water level) hydrograph or stage-discharge curve (rating curve) should be implemented as the downstream boundary condition for proper open channel hydraulic modeling. But, in most real-world cases, HPP operations are controlled by the discharge (when reservoir exists), not water level (stage) and stage-discharge curve is not easy to determine. Therefore, the outflow hydrograph is used as the downstream boundary condition. As can be seen in Figure 6, the outflow hydrograph reflects the daily power production in real-world HPP operations. Thus, this supports the assumption that boundary conditions create, dominantly, the model's uncertainty.

To analyze the influence of PID controllers’ gains, sets of parameter values are used. Parameters are increased step-by-step by the order of magnitude. A set of identical controllers with the same values of parameters (Kp, Ki, and Kd) are used for all assimilation points. Tuning each controller separately requires thorough, independent, investigation, including multi-objective optimization, which will be the subject of a forthcoming research. This paper is focused on determining the initial values and bounds of these parameters, and good to start with.

Analysis is carried out in three phases. In Phase 1, only proportional gain is analyzed (integrative and derivative are switched off), by changing the Kp using the values of 10n, where n is integer from . In Phase 2, Kp is fixed to the value that gave the best DA quality indicators in Phase 1, while Ki is changed using the values of 10n () and derivative gain is switched off. In the last, Phase 3, the Kd is changed using the values of 10n (), Kp is fixed to the same value as used in Phase 2 and Ki is fixed to the value that gave the best results in Phase 2. The set of parameters’ values could be done with more samples and it would probably show an even better solution. The goal of the presented sequential tuning procedure is not to determine the global optimal solution on the controllers’ parameters but to estimate the initial values and bounds preserving the model's stability. Hence, the adopted values are good enough to present the tuning method potential.

RESULTS AND DISCUSSION

The base case scenario is the model's free run, without DA. Figure 7 shows the water levels obtained using the model in free run mode. The water levels obtained show the model's inability to reach the true state, which is most notable at the stations closer to the upstream boundary (Nera and Golubac). This originates from the bad boundary conditions being used and the fact that these stations are more affected by the upstream boundary condition. Stations closer to the downstream boundary are dominantly under the influence of the hydropower plant which practically mitigates all upstream model driving due to big reservoir volume. Using such a model for HPP control operations will lead to numerous problems. Hence, DA is required.

Figure 7

Comparison of the observed water levels and water levels obtained by the model free run simulation.

Figure 7

Comparison of the observed water levels and water levels obtained by the model free run simulation.

Initial tuning analyzed in this paper starts with Phase 1, where only proportional gain is used. Figure 8 shows the changes in DA quality indicators by increasing the Kp. Results are presented for five assimilation points (black bars) and one validation point (Orsova, gray bars). Please note that CorrVol is zero at validation point since there is no controller here to add/subtract the correction flow. RMSE values are pretty much the same for the values of Kp smaller than 100. RMSE values vary around 0.3 m for assimilation points and around 0.37 m for validation point. The value of maxError has the same trend, varying around 0.63 m for assimilation points and 0.77 m for validation point. Setting the value of Kp at 100, RMSE and maxError jump to values around 2.5 and 5, respectively, while increasing the Kp to 1,000 makes the model unstable. This means that the controller's operations cause the model to become first unstable (large water levels oscillations) and then to produce physically impossible solutions. In such situations, RMSE and maxError are assigned a value of 100, reflecting the poor performance of DA. Also, the CorrVol presented in Figure 8 increases for a few orders of magnitude when Kp is increased to 1,000.

Figure 8

Change of the DA quality indicators by increasing the Kp (Ki = 0, Kd = 0).

Figure 8

Change of the DA quality indicators by increasing the Kp (Ki = 0, Kd = 0).

Small values of the RMSE and maxError for the values of Kp below 100 could lead to the wrong conclusion that values around 10 can produce satisfying results. Therefore, AssimTRatio has also to be analyzed. This indicator shows that process error (water level difference) is not below the threshold set at 0.05 m for more than 97% of simulation time. This indicates that using the proportional gain only is not a good option for the application of the controllers in the hydraulic model. Therefore, Phase 2 is required.

In Phase 2, integrative gain is switched on. Values of Ki are changed between 10−3 and 103. Integrative gain is coupled with the proportional gain factor value of Kp = 10. Figure 9 shows changes of DA quality indicators with the increase of Ki. RMSE and maxError values are better than using the proportional gain only. RMSE is less than 0.25 for assimilation points and less than 0.3 for the validation point, for listed Ki values below 100. In these cases, maxError is below 0.5 for assimilation points and below 0.6 for the validation point. Ki set to 10 gives by far the best results for the RMSE and maxError. RMSE is 0.02 m for the assimilation points, and 0.07 m for the validation point. Value of maxError is 0.074 m for assimilation points and 0.187 m for the validation point. Further increase of Ki makes the model unstable (in Figure 8 values of RMSE and maxError are limited to 100). Total correction volume increases by increasing Ki, but not significantly as was the case when only Kp was used. Increasing Ki shows the increase in CorrVol but with a tendency to reach the maximum of the total correction volume. When AssimTRatio is analyzed, by far the best results are obtained when Ki is set to 10. In this case, AssimTRatio is 0.075 for assimilation points, which means that process error has the value over the threshold for only 7.5% of the simulation time. This is an extreme improvement from Phase 1. This value is 0.48 for the validation point, which is bigger than for the assimilation points, yet much better than in Phase 1.

Figure 9

Change of DA quality indicators by increasing Ki (Kp = 10, Kd = 0).

Figure 9

Change of DA quality indicators by increasing Ki (Kp = 10, Kd = 0).

Even though Phase 2 shows significant improvement in DA quality indicators, Phase 3 is used to determine whether the addition of derivative gain could give some further improvement. Factor Kd is included with increasing values from 10−3 to 104 (Figure 10), coupled with constant values of Kp = 10 and Ki = 10. Increasing the value of Kd up to 1,000 shows no improvement in the DA quality indicators. All DA performance indicators have the same values as obtained with the best solution from Phase 2. This can be explained by the fact that properly tuned PI controllers reduce process error to values near 0, so derivatives de/dt in Equation (8) have extremely low values. Multiplication of these values by Kd gives low addition in correction flow thus, practically, having no impact on the assimilation process. Effects of the derivative gain can be seen only if Kd has extreme values, such as 10,000. It increases the impact to the correction flow. In this case study, usage of the derivative gain creates problems, because the model becomes too sensitive to derivative gain output, resulting in the model's instability. This can be explained by the fact that derivative gain output (and total correction flow accordingly) creates highly dynamic input to the model (extremely different from the natural events). These events cannot be propagated by the model, which creates the model's instability. This is not unexpected considering other applications of the PID controllers (80% of the applications omit the derivative gain). Therefore, the multi-metric tuning procedure used in this case study shows that adding the derivative gain gives no benefits. This outcome is case specific, and it should not be taken for granted. Instead, the main contribution of the sequential tuning procedure is showing the potential to narrow the search space for optimal controllers’ parameters. It could improve fine tuning of the controllers using the multi-objective optimization, which will be analyzed in a forthcoming research.

Figure 10

Change of DA quality indicators by increasing Kd (Kp = 10, Ki = 10).

Figure 10

Change of DA quality indicators by increasing Kd (Kp = 10, Ki = 10).

Figure 11 shows the levels of the model simulation coupled with the PI controllers, where proportional and integrative gains’ factors, Kp and Ki, are set to 10. In this case, the model will reach the true state in the first half of a day, which gives the improved water level state at the end of the 7th day needed for the short-term forecasting. Figure 12 shows the correction flow hydrographs for all five assimilation points. When correction flow hydrograph at assimilation point 1 (Nera) is integrated to estimate the correction volume at this point, the total volume of 132.1 × 106 m3 is obtained. Total volume difference between true inflow hydrograph (Figure 6) and the one used to drive the model in the assimilation process is 104.7 × 106 m3. This shows that controller(s) located near the uncertain inflow can reconstruct the ‘missing’ inflows (flow hydrograph – Nera in Figure 12). The rest of the total correction volume CorrVol is the result of usage of inappropriate downstream boundary condition in simulation mode (which happens often). Controllers react to water level difference between true and simulated state even when simulated levels are the product of bad boundary condition. This is even more visible at downstream controllers. They are dominantly under the influence of the hydropower plant operations which makes them more sensitive to the faults in downstream boundary condition. Correction volumes for each station separately (Figure 12) show significant difference between Dubova station and the others. At this station the controller subtracted 829.8 × 106 m3, which is enormous compared to other stations. This indicates that downstream boundary condition strongly affects the results. Therefore, PID controllers have the potential to be used for detection of problems which are the result of bad boundary conditions applied, beside estimated uncertain inflows. This can be a very significant decision support tool, and it will be analyzed in a future research.

Figure 11

Comparison of the observed water levels and water levels obtained by the model coupled with PI controller (Kp = 10, Ki = 10).

Figure 11

Comparison of the observed water levels and water levels obtained by the model coupled with PI controller (Kp = 10, Ki = 10).

Figure 12

Correction flows added (positive values) and subtracted (negative values) at five assimilation points (Nera, Golubac, Dobra, Donji Milanovac, and Dubova) by PI controllers (Kp = 10, Ki = 10) and flow hydrograph downstream of Nera (true inflow corrected by Nera PI controller).

Figure 12

Correction flows added (positive values) and subtracted (negative values) at five assimilation points (Nera, Golubac, Dobra, Donji Milanovac, and Dubova) by PI controllers (Kp = 10, Ki = 10) and flow hydrograph downstream of Nera (true inflow corrected by Nera PI controller).

Figure 12 also shows that some controllers supply the high correction flows at the beginning of the assimilation. Since the simulation model needs some time to respond to those corrections, during this short period the integrative part will accumulate the error causing water flow to overshoot. To overcome this problem at the beginning of the assimilation, hot start for integral error should be adopted as a common practice. This would increase the impact of the integrative gain from the assimilation start.

CONCLUSIONS

This paper presents further insight into the application of control theory-based data assimilation for 1D river hydraulic models used as decision support tool for hydropower production control. PID controllers are applied as a DA tool, in a form of simple (fictive) lateral inflow elements that reduce the water level difference between the model and the observations. The methodology is applied on the real-world case study, HPP Iron Gate with a 170 km-long Danube section with a synthetic test case scenario. The general goal of this analysis is to show the capability of PID controllers to properly update the model state and produce the improved initial conditions for the forecasting. Sequential, multi-metric tuning procedure for the controllers is introduced. It is performed in three phases, considering the DA quality indicators: RMSE, maxError, CorrVol, and AssimTRatio. In Phase 1, controllers are used in a form of P controllers. Multi-metric assessment of the tuning shows that that this type of controller is not able to produce satisfying results based on all indicators. In Phase 2, controllers are upgraded by adding the integrative gain (PI controllers). This upgrade produced highly improved results of model updating, based on all indicators. It also showed that tuning of the controllers’ parameters must be carefully carried out. Inadequate tuning of the controllers can result in model instability. In Phase 3, PI controllers are upgraded to the full form of PID controllers by adding the derivative gain. Multi-metric assessment can show the necessity of adding the derivative gain. If the metrics show no improvement in Phase 3, derivative gain can be switched off. According to the results obtained in this research, the following specific conclusions can be derived:

  • PID controllers can be applied as a DA tool for 1D open channel hydraulic model update.

  • Usage of only one DA quality indicator, such as RMSE, cannot guarantee the proper implementation (tuning) of the controllers.

  • Other DA quality indicators should be used as well, especially AssimTRatio, in combination with RMSE.

  • Multi-metric assessment of controllers’ parameters, used in this research, cannot result in the global optimum (Kp, Ki, and Kd). This procedure provides parameters’ values good to start with, and this procedure should not be omitted, because it can significantly narrow the search space for finding the global optimal solution.

  • More precise tuning (fine tuning) should be conducted using the multi-objective optimization algorithms. The presented sequential tuning procedure should be used to determine initial values of the controllers’ parameters.

Based on the results and previous specific conclusions, control theory-based DA justifies its application in large-scale open channel hydraulic models. However, full application of this methodology requires further testing on the influence of numerical model complexity on the system's performance and procedure for fine tuning of the controllers (each controller individually and all parameters simultaneously), which will be the subject of forthcoming research.

DATA AVAILABILITY STATEMENT

All relevant data are included in the paper or its Supplementary Information.

REFERENCES

Ahmad
S. K.
&
Hossain
F.
2019
A web-based decision support system for smart dam operations using weather forecasts
.
Journal of Hydroinformatics
21
(
5
),
687
707
.
Altinten
A.
,
Ketenanlioglu
F.
,
Erdogan
S.
,
Hapoglu
H.
&
Alpbaz
M.
2008
Self-tuning PID control of jacketed batch polystyrene reactor using genetic algorithm
.
Chemical Engineering Journal
138
(
1–3
),
490
497
.
Andreadis
K. M.
&
Schumann
G. J.-P.
2014
Estimating the impact of satellite observations on the predictability of large-scale hydraulic models
.
Advances in Water Resources
73
,
44
54
.
http://dx.doi.org/10.1016/j.advwatres.2014.06.006
.
Ang
K. H.
,
Chong
G.
&
Li
Y.
2005
PID control system analysis, design, and technology
.
IEEE Transactions on Control Systems Technology
13
(
4
),
559
576
.
Annis
A.
&
Nardi
F.
2019
Integrating VGI and 2D hydraulic models into a data assimilation framework for real time flood forecasting and mapping
.
Geo-Spatial Information Science
22
(
4
),
223
236
.
https://doi.org/10.1080/10095020.2019.1626135
.
Astrom
K. J.
&
Hagglund
T.
1995
PID Controllers: Theory, Design and Tuning
, 2nd edn.
Instrument Society of America
,
Research Triangle Park, NC
,
USA
.
Bonet
E.
,
Gómez
M.
,
Yubero
M. T.
&
Fernández-Francos
J.
2019
Gorosobo simplified: an accurate feedback control algorithm in real time for irrigation canals
.
Journal of Hydroinformatics
21
(
6
),
945
961
.
Bozzi
S.
,
Passoni
G.
,
Bernardara
P.
,
Goutal
N.
&
Arnaud
A.
2015
Roughness and discharge uncertainty in 1D water level calculations
.
Environmental Modeling and Assessment
20
(
4
),
343
353
.
Chiha
I.
,
Noureddine
L.
&
Borne
P.
2012
Tuning PID controller using multiobjective ant colony optimization
.
Applied Computational Intelligence and Soft Computing
1
,
1
7
.
Costabile
P.
&
Macchione
F.
2012
Analysis of one-dimensional modelling for flood routing in compound channels
.
Water Resources Management
26
(
5
),
1065
1087
.
Deng
Y.
,
Liu
Y.
&
Zhou
D.
2015
An improved genetic algorithm with initial population strategy for symmetric TSP
.
Mathematical Problems in Engineering
2015
(
3
),
212794
.
https://doi.org/10.1155/2015/212794
.
dos Santos Coelho
L.
2009
Tuning of PID controller for an automatic regulator voltage system using chaotic optimization approach
.
Chaos, Solitons and Fractals
39
(
4
),
1504
1514
.
http://dx.doi.org/10.1016/j.chaos.2007.06.018
.
Ehteram
M.
,
Mousavi
S. F.
,
Karami
H.
,
Farzin
S.
,
Singh
V. P.
,
Chau
K.-W.
&
El-Shafie
A.
2018
Reservoir operation based on evolutionary algorithms and multi-criteria decision-making under climate change and uncertainty
.
Journal of Hydroinformatics
20
(
2
),
332
355
.
Eurelectric
.
2011
Renewables Action Plan (RESAP)
.
Evensen
G.
1994
Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics
.
Journal of Geophysical Research
99
(
C5
),
10143
.
http://doi.wiley.com/10.1029/94JC00572
.
Fava
M. C.
,
Mazzoleni
M.
,
Abe
N.
,
Mendiondo
E. M.
&
Solomatine
D. P.
2020
Improving flood forecasting using an input correction method in urban models in poorly gauged areas
.
Hydrological Sciences Journal
65
(
7
).
http://dx.doi.org/10.1080/02626667.2020.1729984
.
Fu
M.
,
Fan
T.
,
Ding
Z.
,
Salih
S. Q.
,
Al-Ansari
N.
&
Yaseen
Z. M.
2020
Deep learning data-intelligence model based on adjusted forecasting window scale: application in daily streamflow simulation
.
IEEE Access
8
,
32632
32651
.
Garambois
P. A.
,
Larnier
K.
,
Monnier
J.
,
Finaud-Guyot
P.
,
Verley
J.
,
Montazem
A.-S.
&
Calmant
S.
2020
Variational estimation of effective channel and ungauged anabranching river discharge from multi-satellite water heights of different spatial sparsity
.
Journal of Hydrology
581
,
124409
.
https://doi.org/10.1016/j.jhydrol.2019.124409
.
García-Pintado
J.
,
Neal
J. C.
,
Mason
D. C.
,
Dance
S. L.
&
Bates
P. D.
2013
Scheduling satellite-based SAR acquisition for sequential assimilation of water level observations into flood modelling
.
Journal of Hydrology
495
,
252
266
.
Giustarini
L.
,
Matgen
P.
,
Hostache
R.
,
Monranari
M.
,
Plaza
D.
,
Pauwels
V. R. N.
,
De Lannoy
G. J. M.
,
De Keyser
R.
,
Pfister
L.
,
Hoffmann
L.
&
Savenije
H. H. G.
2011
Assimilating SAR-derived water level data into a hydraulic model: a case study
.
Hydrology and Earth System Sciences
15
(
7
),
2349
2365
.
Hansen
L. S.
,
Borup
M.
,
Møller
A.
&
Mikkelsen
P. S.
2014
Flow forecasting using deterministic updating of water levels in distributed hydrodynamic urban drainage models
.
Water (Switzerland)
6
(
8
),
2195
2211
.
Hostache
R.
,
Lai
X.
,
Monnier
J.
&
Puech
C.
2010
Assimilation of spatially distributed water levels into a shallow-water flood model. Part II: use of a remote sensing image of Mosel River
.
Journal of Hydrology
390
(
3–4
),
257
268
.
http://dx.doi.org/10.1016/j.jhydrol.2010.07.003
.
Jean-Baptiste
N.
,
Malaterre
P. O.
,
Dorée
C.
&
Sau
J.
2011
Data assimilation for real-time estimation of hydraulic states and unmeasured perturbations in a 1D hydrodynamic model
.
Mathematics and Computers in Simulation
81
(
10
),
2201
2214
.
http://dx.doi.org/10.1016/j.matcom.2010.12.021
.
Li
X.-L.
,
H.
,
Horton
R.
,
An
T.
&
Yu
Z.
2014
Real-time flood forecast using the coupling support vector machine and data assimilation method
.
Journal of Hydroinformatics
16
(
5
),
973
988
.
Lund
N. S. V.
,
Madsen
H.
,
Mazzoleni
M.
,
Solomatine
D.
&
Borup
M.
2019
Assimilating flow and level data into an urban drainage surrogate model for forecasting flows and overflows
.
Journal of Environmental Management
248
,
109052
.
https://doi.org/10.1016/j.jenvman.2019.05.110
.
Madsen
H.
,
Rosbjerg
D.
,
Damgård
J.
&
Hansen
F. S.
2003
Data assimilation in the MIKE 11 flood forecasting system using Kalman filtering
.
Water Resources Systems – Hydrological Risk, Management and Development
281
,
75
81
.
Matgen
P.
,
Montanari
M.
,
Hostache
R.
,
Pfister
L.
,
Hoffmann
L.
,
Plaza
D.
,
Pauwels
V. R. N.
,
De Lannoy
G. J. M.
,
De Keyser
R.
&
Savenije
H. H. G.
2010
Towards the sequential assimilation of SAR-derived water stages into hydraulic models using the particle filter: proof of concept
.
Hydrology and Earth System Sciences
14
(
9
),
1773
1785
.
Mazzoleni
M.
,
Verlaan
M.
,
Alfonso
L.
,
Monego
M.
,
Norbiato
D.
,
Ferri
M.
&
Solomatine
D. P.
2017
Can assimilation of crowdsourced data in hydrological modelling improve flood prediction?
Hydrology and Earth System Sciences
21
(
2
),
839
861
.
Mazzoleni
M.
,
When
U.
,
Alfonso
L.
,
Norbiato
D.
,
Monego
M.
,
Ferri
M.
&
Solomatine
D.
2018
Exploring the influence of citizen involvement on the assimilation of crowdsourced observations: a modelling study based on the 2013 flood event in the Bacchiglione catchment (Italy)
.
Hydrology and Earth System Sciences
22
(
1
),
391
416
.
Milašinović
M.
,
Zindović
B.
,
Rosić
N.
&
Prodanović
D.
2019
PID Controllers as Data Assimilation Tool for 1D Hydrodynamic Models of Different Complexity
. In:
Proceedings of the 5th International Conference SimHydro 2019
,
Nice, France
.
Milašinović
M.
,
Prodanović
D.
,
Zindović
B.
,
Rosić
N.
&
Milivojević
N.
2020
Fast data assimilation for open channel hydrodynamic models using control theory approach
.
Journal of Hydrology
584
,
124661
.
https://doi.org/10.1016/j.jhydrol.2020.124661
.
Neal
J. C.
,
Atkinson
P. M.
&
Hutton
C. W.
2007
Flood inundation model updating using an ensemble Kalman filter and spatially distributed measurements
.
Journal of Hydrology
336
(
3–4
),
401
415
.
Ocio
D.
,
Le Vine
N.
,
Westerberg
I.
,
Pappenberger
F.
&
Buytaert
W.
2017
The role of rating curve uncertainty in real-time flood forecasting
.
Water Resources Research
53
(
5
),
4197
4213
.
Ou
C.
&
Lin
W.
2006
Comparison between PSO and GA for parameters optimization of PID controller
. In:
2006 IEEE International Conference on Mechatronics and Automation, ICMA 2006
,
Luoyang, Henan, China
, pp.
2471
2475
.
Rakovec
O.
,
Weerts
A. H.
,
Hazenberg
P.
,
Torfs
P. J. J. F.
&
Uijlenhoet
R.
2012
State updating of a distributed hydrological model with ensemble Kalman filtering: effects of updating frequency and observation network density on forecast accuracy
.
Hydrology and Earth System Sciences
16
(
9
),
3435
3449
.
Randrianasolo
A.
,
Thirel
G.
,
Ramos
M. H.
&
Martin
E.
2014
Impact of streamflow data assimilation and length of the verification period on the quality of short-term ensemble hydrologic forecasts
.
Journal of Hydrology
519
,
2676
2691
.
http://dx.doi.org/10.1016/j.jhydrol.2014.09.032
.
Romanowicz
R. J.
,
Young
P. C.
&
Beven
K. J.
2006
Data assimilation and adaptive forecasting of water levels in the River Severn catchment, United Kingdom
.
Water Resources Research
42
(
6
),
1
12
.
Rosić
N.
,
Jaćimović
N.
,
Prodanović
D.
&
Stojanović
B.
2017a
Data assimilation for operational reservoir management on the Danube River
. In:
7th International Conference on Information Society and Technology ICIST 2017
,
Kopaonik, Serbia
, pp.
210
213
.
Rosić
N.
,
Prodanović
D.
,
Stojanović
B.
&
Obradović
D.
2017b
Near real time data assimilation of numerical simulation model for Danube river from Novi Sad to Iron Gate I, test results
.
Vodoprivreda
49
(
288
),
253
261
(in Serbian)
.
Schütze
M.
,
Campisano
A.
,
Colas
H.
,
Schilling
W.
&
Vanrolleghem
P. A.
2004
Real time control of urban wastewater systems – where do we stand today?
Journal of Hydrology
299
(
3–4
),
335
348
.
Schwanenberg
D.
,
Becker
B. P. J.
&
Xu
M.
2015
The open real-time control (RTC)-tools software framework for modeling RTC in water resources sytems
.
Journal of Hydroinformatics
17
(
1
),
130
148
.
Skogestad
S.
2004
Simple analytic rules for model reduction and PID controller tuning
.
Modeling, Identification and Control
25
(
2
),
85
120
.
Solihin
M. I.
,
Wahyudi
,
Kamal
M. A.
&
 Legowo
A.
2008
Objective Function Selection of GA-Based PID Control Optimization for Automatic Gantry Crane
.
Proceedings of the International Conference on Computer and Communication Engineering 2008
,
ICCCE08: Global Links for Human Development
.
Kuala Lumpur, Malaysia
,
pp.
883
887
.
Stanbury
P. F.
,
Whitaker
A.
&
Hall
S. J.
2017
Instrumentation and Control. In: Principles of Fermentation Technology
.
Elsevier
,
Oxford
,
UK
, pp.
487
536
.
Sun
L.
,
Nistor
I.
&
Seidou
O.
2015
Streamflow data assimilation in SWAT model using extended Kalman filter
.
Journal of Hydrology
531
,
671
684
.
http://dx.doi.org/10.1016/j.jhydrol.2015.10.060
.
Taormina
R.
&
Chau
K. W.
2015
ANN-based interval forecasting of streamflow discharges using the LUBE method and MOFIPS
.
Engineering Applications of Artificial Intelligence
45
,
429
440
.
http://dx.doi.org/10.1016/j.engappai.2015.07.019
.
Thirel
G.
,
Martin
E.
,
Mahfouf
J.-F.
,
Massart
S.
,
Ricci
S.
,
Regimbeau
F.
&
Habets
F.
2010
A past discharges assimilation system for ensemble streamflow forecasts over France – part 1: description and validation of the assimilation system
.
Hydrology and Earth System Sciences
14
(
8
),
1623
1637
.
http://www.hydrol-earth-syst-sci.net/14/1623/2010/
.
Vrugt
J. A.
,
Gupta
H. V.
,
Nualláin
B.
&
Bouten
W.
2006
Real-time data assimilation for operational ensemble streamflow forecasting
.
Journal of Hydrometeorology
7
(
3
),
548
565
.
http://journals.ametsoc.org/doi/abs/10.1175/JHM504.1
.
Wang
Y.
&
Wei
B.
2009
A new particle swarm optimization based auto-tuning of PID controller
.
Journal of Information and Computational Science
6
(
1
),
219
226
.
Wu
C. L.
&
Chau
K. W.
2013
Prediction of rainfall time series using modular soft computing methods
.
Engineering Applications of Artificial Intelligence
26
(
3
),
997
1007
.
Yoon
Y.
,
Durand
M.
,
Merry
C. J.
,
Clark
E. A.
,
Andreadis
K. M.
&
Alsdorf
D. E.
2012
Estimating river bathymetry from data assimilation of synthetic SWOT measurements
.
Journal of Hydrology
464–465
,
363
375
.
http://dx.doi.org/10.1016/j.jhydrol.2012.07.028
.
Ziegler
J. G.
&
Nichols
N. B.
1995
Optimum settings for automatic controllers
.
InTech
42
(
6
),
94
100
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY-NC-ND 4.0), which permits copying and redistribution for non-commercial purposes with no derivatives, provided the original work is properly cited (http://creativecommons.org/licenses/by-nc-nd/4.0/).