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

### PID controllers as data assimilation tool for hydraulic model

*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)):

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

*Q*represents correction flow that reduces the water level discrepancy between observations and model. It is calculated using the PID controller's theory.

_{corr}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*), as stated in Equation (6). In many real-world problems the simulation time step

_{obs}*Δt*is much shorter than the observation time step

*Δt*. Therefore, in the period between two existing observations, the ‘observed’ water level

_{obs}*Z**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

_{obs}*Z**, multiplier

_{obs}*C*, named as attenuation factor (Figure 3), is used. In Equation (7),

*t*is time of the previously available observation,

_{prev_obs}*t*is current simulation time,

*t*is observation time, and

_{obs}*Δt*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.

_{obs}*error*using the

*control variable*. The control variable used to reduce this error is lateral inflow

*Q*(

_{corr}*t*) presented by Equation (8).

PID parameters are: *K _{p}* – proportional gain factor used to multiply the current error value,

*K*– integrative gain factor used to add the influence of error history, and

_{i}*K*– 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

_{d}*K*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

_{p}*K*avoid problems of an unstable system, but time needed for reaching the setpoint increases and, practically, makes system unable to reach the goal. Therefore,

_{p}*K*, 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

_{i}*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*, and

_{i}*K*) should be estimated. Therefore, simple tuning procedures are necessary to complete this step.

_{d}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*and

_{p}*K*which showed the best results.

_{i}*K*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

_{d}*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*can affect the system's stability. This is a solid reason to use the proposed phases in sequential tuning procedure.

_{d}### 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

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

*Z*represents sample from the observed water level time series,

_{j,obs}*Z*represents sample from the simulated water level time series,

_{j,sim}*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

*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.where

*e*is error time series for the

_{i}*i*-th observation location and

*M*is the number of those locations.

#### 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.where

*M*is the number of assimilation points,

_{A}*t*is the time period of the simulation when assimilation is performed, and

_{sim}*Q*is the correction flow.

_{corr}#### Assimilation time ratio

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.

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

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^{9}m^{3}. 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^{3}/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^{9}, with average ‘uncertain’ inflow which is larger than ‘true’ inflow for 104.7 × 10^{6} m^{3}/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*= 7 days.

_{sim}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*, and

_{i}*K*) 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.

_{d}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}*n*is integer from . In Phase 2,

*K*is fixed to the value that gave the best DA quality indicators in Phase 1, while

_{p}*K*is changed using the values of 10

_{i}*() and derivative gain is switched off. In the last, Phase 3, the*

^{n}*K*is changed using the values of 10

_{d}*(),*

^{n}*K*is fixed to the same value as used in Phase 2 and

_{p}*K*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.

_{i}## 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.

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*smaller than 100.

_{p}*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*at 100,

_{p}*RMSE*and

*maxError*jump to values around 2.5 and 5, respectively, while increasing the

*K*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,

_{p}*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*is increased to 1,000.

_{p}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

^{−3}and 10

^{3}. Integrative gain is coupled with the proportional gain factor value of

*K*= 10. Figure 9 shows changes of DA quality indicators with the increase of

_{p}*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*values below 100. In these cases,

_{i}*maxError*is below 0.5 for assimilation points and below 0.6 for the validation point.

*K*set to 10 gives by far the best results for the

_{i}*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*makes the model unstable (in Figure 8 values of

_{i}*RMSE*and

*maxError*are limited to 100). Total correction volume increases by increasing

*K*, but not significantly as was the case when only

_{i}*K*was used. Increasing

_{p}*K*shows the increase in

_{i}*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*is set to 10. In this case,

_{i}*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

^{−3}to 10

^{4}(Figure 10), coupled with constant values of

*K*= 10 and

_{p}*K*= 10. Increasing the value of

_{i}*K*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

_{d}*de/dt*in Equation (8) have extremely low values. Multiplication of these values by

*K*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

_{d}*K*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.

_{d}Figure 11 shows the levels of the model simulation coupled with the PI controllers, where proportional and integrative gains’ factors, *K _{p}* and

*K*, 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

_{i}^{6}m

^{3}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

^{6}m

^{3}. 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

^{6}m

^{3}, 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.

## 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 (

*K*,_{p}*K*, and_{i}*K*). 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._{d}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

*,*

*In:*Principles of Fermentation Technology