Control theory-based data assimilation for hydraulic models as a decision support tool for hydropower systems: sequential, multi-metric tuning of the controllers

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.

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.

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 − Z_model) and the setpoint of the variable (observed water level − Z_obs), as stated in Equation (6). In many real-world problems the simulation time step Δt is much shorter than the observation time step Δt_obs. 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), t_{prev_obs} is time of the previously available observation, t is current simulation time, t_obs is observation time, and Δt_obs 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.

PID parameters are: K_p – proportional gain factor used to multiply the current error value, K_i – integrative gain factor used to add the influence of error history, and K_d – 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 K_p 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 K_p avoid problems of an unstable system, but time needed for reaching the setpoint increases and, practically, makes system unable to reach the goal. Therefore, K_i, 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 (K_p, K_i, and K_d) 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 K_p which produced the best assimilation results according to assimilation quality indicators. In Phase 3, derivative gain is added (PID controller) to the values of K_p and K_i which showed the best results. K_d 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 K_d can affect the system's stability. This is a solid reason to use the proposed phases in sequential tuning procedure.

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).

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

Assimilation time ratio AssimTRatio [/] represents the ratio of period in which process error exceeds the threshold to the total simulation time t_sim. 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.

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.

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 × 10⁹m³. 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 m³/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 × 10⁹, with average ‘uncertain’ inflow which is larger than ‘true’ inflow for 104.7 × 10⁶ m³/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 Δt_obs = 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 t_sim = 7 days.

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 (K_p, K_i, and K_d) 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 K_p using the values of 10ⁿ, where n is integer from ⁠. In Phase 2, K_p is fixed to the value that gave the best DA quality indicators in Phase 1, while K_i is changed using the values of 10ⁿ (⁠⁠) and derivative gain is switched off. In the last, Phase 3, the K_d is changed using the values of 10ⁿ (⁠⁠), K_p is fixed to the same value as used in Phase 2 and K_i 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.

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.

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 K_p. 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 K_p 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 K_p at 100, RMSE and maxError jump to values around 2.5 and 5, respectively, while increasing the K_p 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 K_p is increased to 1,000.

Small values of the RMSE and maxError for the values of K_p 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 K_i are changed between 10⁻³ and 10³. Integrative gain is coupled with the proportional gain factor value of K_p = 10. Figure 9 shows changes of DA quality indicators with the increase of K_i. 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 K_i values below 100. In these cases, maxError is below 0.5 for assimilation points and below 0.6 for the validation point. K_i 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 K_i makes the model unstable (in Figure 8 values of RMSE and maxError are limited to 100). Total correction volume increases by increasing K_i, but not significantly as was the case when only K_p was used. Increasing K_i 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 K_i 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.

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 K_d is included with increasing values from 10⁻³ to 10⁴ (Figure 10), coupled with constant values of K_p = 10 and K_i = 10. Increasing the value of K_d 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 K_d 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 K_d 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 11 shows the levels of the model simulation coupled with the PI controllers, where proportional and integrative gains’ factors, K_p and K_i, 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 × 10⁶ m³ 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 × 10⁶ m³. 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 × 10⁶ m³, 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 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.

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 (K_p, K_i, and K_d). 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.

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