Control theory-based data assimilation for open channel hydraulic models: tuning PID controllers using multi-objective optimization

Water resources management for different purposes (flood protection, hydropower production, water supply, and inland navigation) is under growing pressure due to climate change, population growth, and high urbanization rate (Serdarević & Šabanović 2018; Valipour et al. 2020; Maja & Ayano 2021). Optimal daily use and control of water resources require reliable decision support tools. These tools must provide reliable forecasts of the hydraulic data such as stage and flow hydrographs (Yang et al. 2021). Reliable hydraulic data can be used as a decision-making input eventually leading to a reduction of hydroclimatic extremes and better water resources planning strategies (Huang & Fan 2021; Nie et al. 2021) and can improve hydropower plant operation considering ecological risks (Wen et al. 2021).

Forecasting of the hydraulic data can be successfully conducted using a model-driven approach that relies on physically based hydrologic and hydraulic models. However, physically based models strongly depend on the initial and boundary conditions used for forecasting (e.g., estimated reservoir inflow), which are affected by numerous uncertainties (Bozzi et al. 2015; Ocio et al. 2017). Traditionally, initial conditions in hydraulic models are estimated using the steady-state values of hydraulic data (stage and flow). To improve the model-driven forecasting, initial conditions should be estimated in a hot-start manner, meaning that an appropriate sample from the unsteady process has to be provided (stage and flow across the whole domain at one-time sample). Therefore, observed data, available for the short previous period, are merged with physically based models using data assimilation (DA). DA enables updating of model state according to observed data for better representation of the real conditions. This results in improved initial conditions at the end of the assimilation window. The improved initial conditions are then used for forecasting (Vrugt et al. 2006).

Different DA methods have been applied to solve problems in hydrological and hydraulic modeling where the ensemble Kalman filter (EnKF) method is widely used (Evensen 1994, 2003). This method is successfully applied for improved flood forecasting on large-scale domains (Madsen et al. 2003; Neal et al. 2007; Li et al. 2014; Barthélémy et al. 2017; Jafarzadegan et al. 2021). Along with these applications, EnKF was also applied for improved urban flood forecasting and urban drainage system management (Lund et al. 2019; Kim et al. 2021; Palmitessa et al. 2021). In most of the 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) or streamflow observations (Thirel et al. 2010; Dumedah & Coulibaly 2014; Randrianasolo et al. 2014; Sun et al. 2015). Recently, remotely sensed data, such as water levels and streamflows from satellite images, are used for DA of river hydraulic models (Garambois et al. 2020; Khaki et al. 2020; Li et al. 2020; Mauro et al. 2020; Musuuza et al. 2020; Pujol et al. 2020; Annis et al. 2021). Along with the satellite data, crowdsourced data are used (Mazzoleni et al. 2017, 2018).

All these DA applications in hydrologic–hydraulic modeling show high applicability of the standard DA methods, such as EnKF. However, EnKF and similar DA methods are computationally expensive and often struggle to perform within a reasonable time frame (Madsen & Skotner 2005). This makes them less applicable for solving practical problems related to water resources management. Therefore, simplified tailor-made and time-effective DA methods suitable for solving some specific problems are used (Madsen & Skotner 2005; Hansen et al. 2014; Fava et al. 2020).

For DA application in 1D river hydraulic models, the problem of computational cost can be solved using indirect DA like control theory-based data assimilation (CTDA). In CTDA, water levels are assimilated by adding/subtracting correction flows (Rosić et al. 2017; Milašinović et al. 2018, 2019). These flows are calculated through Proportional–Integrative–Derivative (PID) controllers (Åström & Hägglund 1995), which are significantly time-efficient, when compared to the EnKF method (Milašinović et al. 2020). With its potential for practical and time-efficient application to water-related management problems, it is necessary to further investigate its properties.

PID controllers’ output strongly depends on their parameters. Therefore, determining the controllers’ parameters (tuning of the controllers) is one of the most important steps for any application (Ziegler & Nichols 1995). Various algorithms are proposed for tuning the PID controllers in industrial process applications (Bansal et al. 2012; Borase et al. 2020). Most of these methods use heuristic, trial-and-error approaches for tuning where a single objective (e.g., integral absolute error) is used to evaluate the quality of controllers’ tuning. Recently, optimization-based tuning algorithms have been applied where controller performance measures are used as the objective that has to be minimized (Altinten et al. 2008; Santos Coelho 2009; Gharghory & Kamal 2013). Additionally, when controller's performance is evaluated using multi-metrics, the tuning procedure can include the multi-objective optimization algorithms (Chiha et al. 2012). All these tuning procedures address the application of the PID controllers in traditional, standard conditions such as the control of different industrial processes.

Still, tuning of the multiple, interdependent PID controllers used in a nonstandard application (such as DA for river hydraulic models), where multi-metric performance evaluation is used, has to be thoroughly investigated. Previous research addressing this task (Milašinović et al. 2021) used a heuristic approach for initial multi-metric tuning of the PID controllers. The procedure can be very efficient when there are several controllers applied. It provides homogeneously tuned controllers (all controllers have the same, equal parameters’ values) which will improve the model's DA performance indicators. However, increasing the number of controllers, positioned in different locations along the river, could significantly reduce the performance of equally tuned controllers, making the application of trial-and-error tuning procedure(s) unjustified.

Hence, there is a necessity for an automatic tuning procedure for the PID controllers. Accordingly, the aim of this research is to present improved, multi-objective optimization-based tuning of the controllers used as the DA tool for 1D open channel hydraulic models. This automatic tuning uses DA performance indicators as the optimization objectives in a two-stage procedure. Stage 1 provides initially adjusted controllers and estimates ranges for PID controllers’ parameters using the manual multi-metric tuning approach (Milašinović 2020; Milašinović et al. 2021). This paper introduces Stage 2, which presents further tuning of PID controllers’ parameters using the multi-objective optimization with Nondominated Sorting Genetic Algorithm II (NSGA-II) (Deb et al. 2002). The optimization test objectives included different combinations of DA performance indicators. The NSGA-II produced the range of optimal PID parameters for each individual controller.

CTDA for open channel flows is based on updating water levels in the hydraulic model. For practical reasons, a brief explanation of the PID controllers’ implementation in 1D hydraulic models is presented in this paper. The CTDA method, in detail, can be found in previous research addressing this topic (Milašinović et al. 2020, 2021).

Here, x (m) is the space coordinate, t (s) is time, A (m²) is the cross-sectional area, Q (m³/s) is the discharge, q (m/s) is the lateral inflow, Z (m) is local water surface elevation (water level), g (m/s²) is the acceleration due to gravity, R (m) is hydraulic radius calculated as the ratio between the cross-sectional area and the wetted perimeter, β (/) is the velocity distribution coefficient, n (m^−1/3 s) is Manning's roughness, which is calculated according to Costabile & Macchione (2012), and B (m) is the top width of the cross-section (water surface width).

In Equation (3), index i represents the spatial location (ith cross-section), B (m) is the top width of the cross-section (water surface width), t in superscript denotes current time, Δx (m) is spatial resolution, and Δt (s) is temporal resolution.

In Equations (5)–(7), following variables are used:

e (m) – process error (difference between measured and modeled water levels),

Z*_obs (m) – measured water levels (* denotes that measured values can be temporally interpolated when it is necessary),

Z_model (m) – computed water levels (from the hydraulic model),

Δt_obs (s) – the period between two observations,

t_{prev_obs} (s) – time of the last observed data,

C (/) – attenuation factor [0–1] used to reduce the impact of the interpolated water levels (it is assumed that these values are less reliable than really observed data),

K_p (/) – proportional gain factor,

K_i (/) – integrative gain factor, and

K_d (/) – derivative gain factor.

PID controllers used as DA tools are applied on each location in the model where observed (measured) water-level data are available. These locations are named assimilation locations. For each controller, the optimal values of PID gain factors (K_p, K_i, and K_d) are needed. The values of gain factors are computed using the tuning procedure.

PID controllers’ tuning procedure is assessed using an adequate performance indicator (Åström & Hägglund 1995; Bansal et al. 2012; Borase et al. 2020). The tuning procedure tends to minimize the difference between a process value (modeled water levels) and a setpoint (observed, measured water levels), oscillations in process error, and settling time. When PID controllers are used as a DA tool, four tuning indicators are proposed in the previous research (Milašinović et al. 2021). Root mean square error (RMSE) represents the difference between a process value and a setpoint, the 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 a 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 of excellent assimilation performance (perfectly tuned PID controllers).

Assimilation time ratio AssimTRatio (/) represents the ratio of the period in which process error exceeds the threshold to the total assimilation time t_sim. When AssimTRatio ratio is 0, it shows excellent DA performance (best value, perfectly tuned controllers), while value 1 shows poor DA performance (worst value), where assimilation fails to correct the model toward measurements. To evaluate the AssimTRatio, the error time series has to be transformed into the error duration curve (Milašinović et al. 2021). As there are several observation locations where assimilation quality has to be estimated, AssimTRatio is calculated as the average for all locations

The threshold value used to estimate AssimTRatio can vary and it should be carefully selected. Using a very rigid threshold (e.g., less than 0.001 m) can lead to the underestimation of DA performance according to the AssimTRatio. The threshold value must be selected according to the conditions where CTDA is applied. Some of the conditions that have to be considered to select an appropriate threshold value are average depth and threshold ratio, water-level measurement uncertainty, measurement location (e.g., river bay or open shore), and wind effects for large-scale applications. River water-level measurement on the river and data acquisition is commonly conducted using some types of sensor (e.g., pressure sensor). In that case, water-level measurement uncertainty can be up to 10 cm, depending on the sensor location and velocity (Larrarte et al. 2021). Considering this limitation of the water-level sensors, the error threshold value used in this research is set to 0.05 m.

A manual, sequential tuning procedure for PID controllers used as a DA tool presented in previous research (Milašinović et al. 2021) provides acceptably well-tuned controllers and the possible range of the controllers’ parameters preserving the model stability. The procedure also gives a set of homogeneously tuned controllers where each one has the same values for PID parameters (K_p, K_i, and K_d). Even though this tuning approach gives satisfactorily good results in most cases, DA performance can significantly drop by increasing the number of assimilation locations (the number of controllers applied). Solving this issue requires an automatic tuning procedure. This paper presents an improvement of the previous method in a two-stage tuning procedure (Figure 2). Here, multi-objective optimization is used to improve the initial tuning solution provided by the manual procedure. Generally, multi-objective optimization can be conducted using two types of methods: scalarization and Pareto (Cui et al. 2017; Gunantara 2018). Methods using scalarization would require assigning the weights for each objective function and combining them into the single-objective function. Additionally, this approach requires the normalization of the objective functions. Evaluating different weight combinations would require extensive analysis which would change the focus of this research. Therefore, here presented multi-objective optimization problem is solved using the Pareto-based optimization algorithm, Nondominated Sorting Genetic Algorithm (NSGA-II; Deb et al. 2002), which is widely used for solving multi-objective optimization tasks for water resources (Artina et al. 2012; Darvishi & Kordestani 2019; Gao et al. 2019; Wang et al. 2019; Gaur et al. 2021). Here, this algorithm is used to improve the initial tuning solution provided by the manual procedure.

A sequential tuning procedure where all controllers will have the same values of PID parameters is considered as Stage 1. Obtained parameters are then used as the initial population for Stage 2, where parameters for each controller are improved using the NSGA-II multi-objective optimization algorithm. The result of Stage 2 is a set of heterogeneous controllers which can be additionally re-tuned from time to time using the previous parameters as initial populations.

Previous research (Milašinović et al. 2021) introduced four indicators (RMSE, maxError, AssimTRatio, and CorrVol) for assessing the quality of DA performed using the PID controllers. These indicators are suitable for assessing the multi controllers’ performance when the setpoint is nonstationary (setpoint in CTDA is stage hydrograph). Indicators are used to represent the average discrepancy between model and observed data (RMSE), overshooting (maxError), and settling time (AssimTRatio). A CorrVol indicator is method-specific and describes the ‘intensity’ of the controllers’ intervention in the model. These indicators have to be minimized in Stage 2 of the tuning procedure. Some of these indicators could be opposed, while some of them could be of the same nature (minimizing one leads to the minimization of the other). To identify the indicators suitable for automatic tuning of the controllers, different combinations of the objective combinations are used (Table 1). Combining two objectives in multi-objective optimization makes it a lot easier to analyze the behavior of each one and to identify (potentially) dominant objectives or those which are not suitable for automatic tuning procedure.

A two-stage tuning procedure for PID controllers in the CTDA method is tested for the transboundary hydropower system Iron Gate 1 (IG1) on the Danube River shared between Serbia and Romania. Everyday operations of hydropower plants (HPPs) require short-term forecasting of water levels in the river section, which is affected by hydropower plants. For improved forecasting (improved initial conditions), DA is required.

The hydraulic model for the Danube River section is 170 km long and has the hydropower plant IG1 representing downstream boundary conditions (Figure 3). The domain is discretized by 189 cross-sections, with an average distance of 900 m between Jaroslav Černi 2019 and 2020. The river model is developed using the diffusion wave model with one upstream inflow, neglecting all tributaries (lack of data). This supports the assumption that unreliable boundary conditions (inflow data) are a dominant source of uncertainty in the model. 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 locations (PID controllers are applied at these locations), while Orsova (closest to the downstream boundary) is used as a validation point.

A 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 solid line in Figure 4). The total inflow volume for the 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 circle markers in Figure 4). The 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 arbitrary 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 generate ‘true’ water levels). This also represents the case of sampling the inflow data on a poor timescale, which results in big discrepancies. The altered inflow is set as an upstream boundary, with a total volume of 2.268×10⁹ m³, with an average ‘uncertain’ inflow which is larger than ‘true’ inflow of 104.7×10⁶ m³ (Figure 4). At the downstream boundary, the measured outflow and stage (dashed line and red dashed line with square markers are shown in Figure 4) are used. Stage hydrograph is used to generate the true state, while outflow hydrograph is used as a 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 the computation time step is, in most cases, smaller than the observation time step. The 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. On the other hand, in most real-world cases, HPP operations are controlled by the discharge (when a reservoir exists), so water-level (stage) and stage–discharge curve are not easy to determine. Therefore, the outflow hydrograph is used as a downstream boundary condition in the model exploitation phase. This also induces big discrepancies between observed water levels and model results. As it can be seen in Figure 6, the outflow hydrograph reflects the daily power production in real-world HPP operations.

The multi-objective tuning procedure is evaluated on the iron gate hydropower system test case (Figure 3). In Stage 1 of the two-stage tuning, the procedure was determined that the derivative gain of the PID controller has no impact on the improvement of DA quality indicators (Milašinović et al. 2021). Increasing the value of derivative gain factors K_d (Figures 11, 12, and 13 in.supplementary materials), for this case study, results in abnormal water-level fluctuations and leads to model instability. Therefore, this parameter is set to 0. Using the estimated values of proportional and integrative gain factors (K_p and K_i) from Stage 1, the initial population representing PID controllers’ parameters is created for Stage 2. Further tuning is provided using the multi-objective optimization algorithm NSGA-II for different combinations of two-objective functions (Table 1). The goal of the multi-objective tuning procedure (Stage 2) is to find the range of optimal values for K_p and K_i parameters for each assimilation location in the case study (each PID controller individually). Here, NSGA-II is applied using the population size of 100. Ranges of the K_p and K_i are bounded by 0 at the lower bound and 100 at the upper bound, according to the results obtained in Stage 1 (see Supplementary Material). Initial values of the PID parameters (initial population) are set to 10 (both K_p and K_i) for each of the five assimilation locations.

The NSGA-II-based tuning procedure found optimal, nondominated values of objective functions – Pareto fronts (Figure 5) and connected ranges of optimal values for K_p (Figure 6) and K_i (Figure 7) at each assimilation location. Ranges of optimal values for K_p and K_i are represented by the min, mean, median, and max values (Figures 6 and 7 and Table 2 in Supplementary Material).

When ranges of optimal values for K_p are analyzed with RMSE and AssimTRatio, objectives combination OC3 results in the narrowest ranges (Figure 6). These ranges (integer boundaries) are 51–61 for Nera assimilation location, 46–48 for Golubac, 65–70 for Dobra, 31–37 for D. Milanovac, and 61–65 for Dubova. Similarly, optimal values for K_i have the narrowest range when objective combination OC3 is used: 32–33 for Nera assimilation location, 69–72 for Golubac, 5–10 for Dobra, 25–28 for D. Milanovac, and 41–42 for Dubova. When it comes to objective values, for this OC3 objective combination (Figure 5(c)), it is obvious that the NSGA-II algorithm created the Pareto front where the range of objectives’ values varies around 0.031 for AssimTRatio and around 0.021 m for RMSE. Practically, it signalizes that these two objectives are not conflicted but have the same nature (reducing one of them leads to reduction of the second one and vice versa) and a global optimum can be found.

Generally, PID controllers’ output can be modified to reduce the max value of the error process (maxError). This reduction is generated by ‘slowing down’ the controller, which results in increasing the settling time (here described by the AssimTRatio indicator). This means that objectives RMSE and AssimTRatio are opposed to the maxError indicator. When RMSE and AssimTRatio are combined with maxError (objectives combinations OC1 and OC4), ranges of optimal values for PID parameters are wider (Figures 6 and 7). In these cases, K_p values are 44–61 (OC1) and 27–42 (OC4) for Nera assimilation location, 55–61 (OC1) and 34–76 (OC4) for Golubac, 61–70 (OC1) and 41–68 (OC4) for Dobra, 37–44 (OC1) and 39–56 (OC4) for D. Milanovac, and 59–68 (OC1) and 21–37 (OC4) for Dubova. K_i values are 32–34 (OC1) and 32–34 (OC4) for Nera assimilation location, 51–61 (OC1) and 46–62 (OC4) for Golubac, 10–16 (OC1) and 10–18 (OC4) for Dobra, 7–21 (OC1) and 19–25 (OC4) for D. Milanovac, and 41–43 (OC1) and 40–42 (OC4) for Dubova. Combining RMSE or AssimTRatio with the maxError indicator provides slight improvement of objectives when it is compared to values of the indicators obtained in Stage 1. When OC1 objective combination is used, RMSE varies around 0.021 m, which is pretty much the same when it is compared with the value of this indicator obtained in Stage 1, while maxError is in the range of 0.073–0.077 m comparable to the results from Stage 1. The usage of OC4 objective combination (combining the AssimTRatio and maxError indicators) shows improvement in the values for AssimTRatio. Stage 1 (see Supplementary Material) provided the solution of the PID parameters where AssimTRatio was 0.075. This means that process error (difference between measured and modeled water levels) is below the threshold for more than 92.5% of the assimilation period. When the tuning procedure is extended into Stage 2 using the NSGA-II, this indicator is reduced by 50% that is approximately 0.031–0.032. It means that the process error is below the threshold for more than 96.8% of the assimilation window. The lack of any significant improvement in DA performance indicators (only AssimTRatio displays a noticeable improvement) can be explained by the short period used for assimilation. In addition, a small number of assimilation locations are used. In this case, it is possible to manually tune PID controllers and reach some optimal values. The more complex case, with a longer period for assimilation and tuning of the controllers, will be the subject of future research.

The application of CorrVol indicator in the tuning procedure shows a significant impact on the multi-objective optimization results. In each objective combination where CorrVol appears, there is a significant deterioration of the assimilation results (Figure 8). In OC2 combination, optimization algorithm found Pareto front where CorrVol indicators vary between 0.4×10⁹ and 1×10⁹ m³. Minimization of the CorrVol indicators increases RMSE values eight times that are approximately 0.021–0.181 m compared to results when OC1 and OC3 objectives are used. In this case (OC2), maxError and AssimTRatio are used as independent indicators (not used for optimization). Minimizing the CorrVol also increases AssimTRatio and maxError 10 times approximately (Figure 8). The same trend occurs in OC5 and OC6 combinations. The increase of the RMSE, AssimTRatio, or maxError ranges, when these indicators are used as the objectives combined with CorrVol, provides the results of very low quality (the model is unable to reach observed data). Minimization of the CorrVol in the tuning procedure leads to a significant increase in all other DA performance indicators, no matter if they are used as the objective function or as an independent indicator (Figure 8). The best results, in objective combinations where CorrVol appears as one of the objectives, are obtained when optimal values of the PID parameters are represented by the min value from the range. Using any other optimum representative (mean, median, or max value) leads to a significant decrease in the DA quality (Figure 8). This indicates that CorrVol should not be used as the objective for tuning of the PID controllers.

The CTDA method assumes that unreliable boundary conditions dominantly create a discrepancy between modeled and measured water levels. This means that the water volume which has to be added/subtracted to the model in the CTDA method (in order to keep the model and measurements assimilated) is predetermined and it cannot be minimized. Therefore, CorrVol cannot be utilized as the objective in multi-objective optimization because it leads to significant deterioration of all other indicators. Obtained results (for different objective combinations) show that the NSGA-II algorithm used in a two-stage tuning procedure can provide satisfying results only when appropriate criterion functions are used in multi-objective optimization. Results show that process error during the assimilation window varies between −0.1 and 0.1 m (Figure 9) for controllers optimized using objective combinations OC1, OC2, and OC4, considering even validation station Orsova (no PID controller at this location).

Even though the CorrVol should not be used as the objective for tuning of the controllers, it should not be neglected as the DA performance indicator. When the PID controllers are well-tuned, CorrVol can be calculated for each assimilation location separately. Spatially distributed CorrVol can be used to identify the location of the biggest discrepancy between the model and measurements. For the Iron gate synthetic test case analyzed in this paper, this approach shows that most of the correction flow volume is generated at the downstream section, which is near the hydropower plant boundary condition (Figure 10). Locations closest to the downstream boundary condition (Dubova and Orsova) have a greater process error range than other stations. This information can help in adjusting the downstream boundary condition. Besides this, the assimilation location nearest to the upstream boundary conditions (Nera station) shows a significant discrepancy between modeled and measured data. This is the consequence of the unreliable inflow estimation, and it indicates the necessity for improved inflow forecasting. In theory, CorrVol could be dominantly distributed at the assimilation locations in the middle of the model domain. It can indicate that some significant tributary is neglected in the model, and this internal boundary condition should be better implemented.

Therefore, the pre-processing phase should be performed where CorrVol should be considered to identify the locations of dominant uncertainty sources and eventually reduce it. After this identification, appropriate model reconfiguration should be conducted. For example, if the downstream boundary condition dominantly creates uncertainty in the model, better rating curves have to be estimated, or it would require a representation of the downstream boundary in more detail (detailed equation representing hydropower plant instead of a rating curve). Upstream boundary condition uncertainty is mainly driven by unreliable inflow forecasting. To reduce this reservoir inflow uncertainty, for example, more reliable hydrologic models should be used, or even data-driven inflow forecasting should be considered.

Finally, even though CTDA has been proved as the time-efficient method compared to ensemble-based DA (Milašinović et al. 2020), the multi-objective tuning procedure slows it down. Using iterative tuning (in the optimization process) increases computational time (depending on the hardware configuration) which reduces the possibility of using this tuning approach on a daily basis. This could be considered as the main disadvantage of the proposed tuning procedure.

Instead of running the optimization algorithm every day (or every hour, depending on the timescale), PID controllers should be tuned in a pre-processing phase. This means that controllers should be re-tuned from time to time using a longer assimilation window (more data can be obtained in the period between two tunings) before it is applied for (near) real-time DA. In this case, the performance of the PID controllers should be monitored in (near) real-time using DA quality indicators with assigned trigger values. When some of the DA quality indicators drop below the assigned trigger value, the tuning procedure should be restarted. In the period between two tunings, PID controllers should use parameters’ values obtained when the last tuning was conducted.

This paper presents in detail the tuning procedure required to estimate PID controllers’ parameters when they are used as the DA tool for 1D hydraulic models. Relying on the results of the previous research where CTDA was developed with a manual procedure of adjusting all PIDs with the same parameters (Stage 1), the tuning procedure is now extended with the multi-objective optimization (Stage 2), so each PID will have its own optimal parameters (proportional gain factor K_p and integrative gain factor K_i) according to local conditions at assimilation location. The presented Stage 2 tuning is performed using the multi-objective optimization algorithm NSGA-II with different combinations of two objectives. The objectives are selected from DA performance indicators introduced in previous research: RMSE, maxError, AssimTRatio, and CorrVol.

The output from Stage 2 provides a set of heterogeneously tuned controllers, where a range of optimal values is given for each PID parameter (K_p and K_i) at all assimilation locations. Based on the results obtained in this research, the following specific conclusions can be derived:

• Multi-objective optimization using the NSGA-II algorithm can be successfully used for tuning of the PID controllers used as DA tools for 1D hydraulic models.

• The success of the multi-objective optimization-based tuning procedure depends on the selection of performance indicators used as the optimization objectives.

• CorrVol, as the method-specific DA performance indicator, cannot be used as the tuning objective, because it is predetermined. The minimization of this objective significantly increases all other DA performance indicators.

• Instead of using the CorrVol in the tuning process, this parameter should be used for identifying the dominant uncertainty sources. Based on this identification, further improvements in the hydraulic model can be conducted.

Specific conclusions from this research show that the CTDA method provides good DA results if all tuning procedure steps are completed. Considering the main outputs from this research and previous research on the application of CTDA, a few more steps should be clarified to complete the framework for the implementation of the PID controllers as the DA tool for 1D open channel flow models. As it is derived from this paper, the two-stage tuning procedure should be restarted occasionally, especially when model inputs significantly differ from the one used for the previous tuning. This procedure should be also done for an extended assimilation window to show the full capacity. This would require a definition of a straight-forward algorithm for re-tuning the PID controllers used for (near) real-time DA. The completed framework should result in improved initial conditions (as the main output when the CTDA method is used), which can be used for further forecasting of the hydraulic data. However, reliable near-future forecasting to help decision-making in water resources management should be considered as the ultimate goal. This means that the CTDA method improves only one aspect of the forecasting (initial data). Forecasting results are still affected by unreliable boundary conditions in the forecasting window (e.g., unreliable inflows). Therefore, an extension of the DA into the forecasting window should be also considered as the next step in CTDA application, which will be a subject of forthcoming research.

The main goal of this paper was not to provide the best solution for PID parameters used as a DA tool for the Iron gate case study. To demonstrate the two-stage tuning procedure and identify the crucial indicators that can be used as the objectives in tuning multi-objective optimization, a synthetic test case scenario was used. Using the proposed two-stage tuning procedure will help in further real-world applications of the presented tailor-made data assimilation method applied in the study: Hydroinformatic system for Iron Gate Hydropower Plant done by Jaroslav Černi Water Institute and University of Belgrade – Faculty of Civil Engineering, during 2019–2020, for Serbian Electricity Company (Published in Serbian).

This work was also supported by the Serbian Ministry of Education, Science and Technological Development (Agreement No. 451-03-68/2022-14/200092 and Agreement No. 451-03-68/2022-14/200122).

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

The authors declare there is no conflict.