Model based control for run-of-river system . Part 1 : Model implementation and tuning

Optimal operation and control of a run-of-river hydro power plant depends on good knowledge of the elements of the plant in the form of models. River reaches are often considered shallow channels with free surfaces. A typical model for such reaches use the Saint Venant model, which is a 1D distributed model based on the mass and momentum balances. This combination of free surface and momentum balance makes the problem numerically challenging to solve. The finite volume method with staggered grid was compared with the Kurganov-Petrova central upwind scheme, and was used to illustrate the dynamics of the river upstream from the Grønvollfoss run-of-river power plant in Telemark, Norway, operated by Skagerak Energi AS. In an experiment on the Grønvollfoss run-of-river power plant, a step was injected in the upstream inlet flow at Årlifoss, and the resulting change in level in front of the dam at the Grønvollfoss plant was logged. The results from the theoretical Saint Venant model was then compared to the experimental results. Because of uncertainties in the geometry of the river reach (river bed slope, etc.), the slope and length of the varying slope parts were tuned manually to improve the fit. Then, friction factor, river width and height drop of the river was tuned by minimizing a least squares criterion. The results of the improved model (numerically, tuned to experiments), is a model that can be further used for control synthesis and analysis.


Introduction 1.Background
There are in total over 1500 hydro power plants operating in Norway today, with a total capacity of more than 28 GW.Hydropower plants had an annual production of 129 TWh in 2013, which constituted some 96 % of the production of electricity.Most hydropower production comes from large hydropower plants with water reservoirs.Small hydropower plants is a generic term for plants that have small production capacity.There are in total operating over 1100 such smaller hydropower plants in Norway today and production from these accounts for about 7% of the total hydropower production (Vytvytskyi et al., 2015).Skagerak Energi is operating two run-of-river hydro power plants in river Tinnelva flowing out of lake Tinnsjøen in Eastern Telemark.The uppermost, Årlifoss, has a 16.2 m water fall, and provides inlet flow to the pondage of Grønvollfoss; the Grønvollfoss power plant is studied here.Grønvollfoss is situated 5 km downstream from Årlifoss, and has a water fall of 22.6 m.The installed power in Grønvollfoss is 2x16 MW (i.e. two turbines), and the average annual production is 172 GWh.
It is of interest to control the operation of the Grønvollfoss power plant in such a way that the level in front of the dam is as high as possible in order to maximize the available power, while variation in the level is minimized to avoid loss of water flowing over the dam.The level depends both on the power production/turbine flow and on the upstream inflow of water.For control studies, it is imperative to have available a simulation model of sufficient quality.Therefore, it is of interest to develop a dynamic model for such a run-of-river system in a suitable form for control studies, which is solved efficiently with sufficient numerical accuracy, and which gives good representation of the real power plant.

Previous work
The dynamic model for such run-of-river system has an important role in control and optimization.In previous work (Vytvytskyi et al., 2015;Vytvytskyi, 2015), a model for the Grønvollfoss power plant was developed using the Saint Venant equations, was discretized using the staggered grid scheme, and was tuned to experiments.A case study of the run-of-river power plant using the Kurganov-Petrova second order central upwind scheme is described in Sharma (2015), with more details about the development of this scheme in Kurganov and Petrova (2007).Discretization of the Saint Venant equations using the staggered grid scheme has been studied in Farina et al. (2011).The Quick and Superbee algorithms are more complex, but better algorithms for solving the Saint Venant equations (Johansen et al., 2012).Even better are adaptive WENO schemes (Johansen et al., 2012).In Xu et al. (2012), the problem of estimating the levels of three cascaded river reaches by using the Saint Venant model has been described.

Overview of paper
In this work, the finite volume method with staggered grid (Vytvytskyi et al., 2015;Vytvytskyi, 2015;LeVeque, 2002) is compared with the Kurganov-Petrova central upwind scheme (Sharma, 2015;Kurganov and Petrova, 2007), and the chosen scheme compared to experiments on the system with the model tuned by changing the river bed (slope), the friction factor, the width and the height drop of the river.The paper is structured as follows: Section 2 consists of a system description of the run-of-river power plant together with a description of the experiment on the power plant.Section 3 describes formulation and discretization of the model.Model simulations with comparison of the discretization schemes are presented in Section 4. Experimental results and their application for tuning of the model are provided in Section 5. Finally, conclusions are given in Section 6.

System description
A typical run-of-river power plant is depicted in Figure 1.As described in Vytvytskyi et al. (2015); Vytvytskyi (2015), this type of hydropower plant has a high volumetric flow rate, which passes through the turbine and spins the rotor.The turbine transfers the kinetic energy from water to rotational mechanical power, which will then be transferred to electric power in the generator.After the turbine, the water returns to the river via a short tailrace.The Grønvollfoss power plant is situated in river Tinnelva and is the second power plant in this river.Thus, this power plant does not have a dedicated reservoir; instead the outlet flow from the Årlifoss power plant is an input flow to the pond in front of the Grønvollfoss dam.The Årlifoss plant is situated ca. 5 km upstream in the river.

Geometry
The following information about the system was provided by Skagerak Energi AS or was found from maps on the Internet (Vytvytskyi et al., 2015;Vytvytskyi, 2015).All of the provided parameters are approximations.The water level before the Grønvollfoss dam is around 17.1 meters above the river bed (144.5 m.a.s.l: meters above sea level).1The height difference between the two power plants ( Årlifoss and Grønvollfoss) is around 17.5 meters.The length of the river between the two power plants is L = 5 kilometers, and the width of the river is on average w=180 m.The slope of the river bed and the level of the water surface are schematically shown in Figure 2. As seen, there is a steeper drop of the river bed half way down the river reach, ca.2.5 km down streams from Årlifoss.
The bending and width variation of the river can be seen from the map, screen dump of which is shown in Figure 3.

Experiment
To check and improve the accuracy of the run-of-river model, it is useful to compare the model and tune it against experimental data from the Grønvollfoss power plant (Vytvytskyi et al., 2015;Vytvytskyi, 2015).The experiment that was carried out consists in changing the inlet flow rate to the system Vin for some time (the    1. It should be noted that the water level ahead of the Grønvollfoss dam was decreased by 15 cm before the experiment (to 144.35 m.a.s.l.) to avoid overflowing the dam after increasing the inlet flow rate.In addition, the inlet volumetric flow rate specified indirectly by increasing or decreasing power production and thus inferring the flow rate from the power production, thus the values of the flow rates are approximate values.Also, increasing or decreasing the power production is not instantaneous, but takes place in steps 30 and 20 sec, respectively.The water level ahead of the Grønvollfoss power plant is measured by a level sensor every 100 ms.These values are filtered in a moving average filter using 300 past values, and the filtered value is logged every 5 s.

Model presentation
A common model for shallow water systems is given by the Saint Venant equations, consisting of two coupled partial differential equations.These equations describe the behavior of the river and consist of both mass and momentum balances, and are challenging to solve numerically.A common formulation of these Saint Venant equations for channel is as follows (Vytvytskyi et al., 2015;Vytvytskyi, 2015): with the boundary conditions: Here, A is the wetted cross section area, V is the discharge (the volumetric flow rate), g is the gravitational acceleration, S 0 is the bed slope, S f is the friction slope and w is the width of the river (rectangular channel).
In vector form, a standard formulation is (Sharma, 2015): where: The Saint Venant equations are hyperbolic, and the Finite Volume method is common for solving such hyperbolic equations (LeVeque, 2002).Two different schemes for solving this method are presented below.

The staggered grid scheme
The staggered grid is a scheme where the grid for the mass balance is shifted with regards to the grid for the momentum balance.The detailed development of this scheme is given in Vytvytskyi et al. (2015); Vytvytskyi (2015) and here the results of discretizing the Saint Venant equations using the Finite Volume method with staggered grid are presented.The equations for the water level above the river bed in each cell become: Here, h is the river depth, x is the cell size and depends on the number of cells N V .The equations for the volumetric flow rate in each cell are as follows: Here, θ is the slope of the river bed, Ā is the averaged wetted cross section area, φ is the averaged wetted perimeter, C is the averaged Chézy friction coefficient, which is defined in relation to the Strickler's friction factor and the hydraulic radius, Ṁi,i  ρ and Ṁo,i ρ are the momentum input and output flows per density and are defined as follows:

The Kurganov-Petrova central upwind scheme
The Kurganov-Petrova scheme is a second order scheme which is well balanced and maintains the positivity of the height.The detailed development of this scheme is given in Kurganov and Petrova (2007) and a case study of the run-of-river power plant with using this scheme is described in Sharma (2015).The result of discretizing the Saint Venant equations using the Finite Volume method known as the semi-discrete Kurganov-Petrova scheme is: where Ūj is the cell center average values, H j± 1 2 (t) is the central upwind numerical fluxes at the cell interfaces, defined as: , where a ± j± 1 2 are the one-sided local speed of propagation and are defined in relation to the water depth and the velocity of the water, B is the river bed height.Here, vectors U, F and S are written in the following form: .
It should be noted that z = h+B is the water surface height and q = Vi /w is the discharge per unit width of the river.

Operational data and parameters
Normally, Finite Volume methods are discretized in both time and space, and then solved recursively/by iteration.Here, the model is only discretized in space, leading to a number of ordinary differential systems.Thus, after the model for the run-of-river system is discretized, both schemes are implemented in MATLAB and solved using an ode solver.Model parameters are given in Table 2.

Simulations
Some basic simulations are done to compare the two methods for solving the Saint Venant model.All parameters are as given in Table 2.The simulation has the same scenario as in the experiment and is started from steady state, and after 10 min the inlet volumetric flow rate Vin is suddenly changed from the normal value, which is 120 m 3 /s, to 160 m 3 /s.This new flow rate is maintained for 15 min, and then the inlet flow rate is reduced back to 120 m 3 /s.The increase and decrease in the inlet flow rate is not instantaneous, and in the simulations these are ramped over 30 and 20 sec, respectively.
The results of the simulation of the model with the two methods are shown in Figure 4, where the top plot shows the water level ahead of the Grønvollfoss dam as a function of time.The bottom plot shows how volumetric flow rate in the last cell (ahead of the dam) changes with time.Here, the ode23 solver in MATLAB has been used; this is an implementation of an explicit second order Runge-Kutta method with variable step length4 and gives a balance in the temporal and spatial accuracy.
From Figure 4, the simulation results show that it takes approximately 11-12 min from the change in the inlet flow rate to the water level ahead of the dam starts to increase.Steady state is reached after more than one hour, when the oscillations have decayed.These oscillations describe a wave, which is initiated by the inlet flow jump, and with reflections from the Grønvollfoss dam back and forth to Årlifoss.
Figure 4 also shows that the results for the water level ahead of the dam from both schemes are almost the same.Only by zooming in on the plot for the water level ahead of the dam, shown in Figure 5, it can be seen that the finite volume method with staggered grid (red line) shows more noisy (small oscillations) results than the results from the Kurganov-Petrova central upwind scheme (blue line).The same can also be seen from zoomed plot of the volumetric flow rate in the last cell (Figure 6), where the central upwind method shows some numeric noise (small oscillations), but much less pronounced than the noise in the finite volume method with staggered grid.
Another ode solver in MATLAB reduces the noise in the simulation results: simulating the same two discretizations of the Saint Venant equations using the ode23t solver gives the results shown in Figure 7.This solver is an implementation of the trapezoidal rule using a 'free' interpolant and is used if the problem is only moderately stiff and a solution without numerical damping is needed5 .Figure 7 shows that with ode23t, the noise is only visible in the results for the volumetric flow rate, and only in the staggered grid scheme.For better comparison, more detailed plots are shown in Figure 8 and Figure 9, where the zoomed plots of the water level ahead of the dam and the volumetric flow rate in the last cell are presented, respectively.
In spite of better simulation results (less noise), the ode23t solver has one disadvantage: the simulation time is approximately 4 times longer with the ode23t than with the ode23 solver.
In summary, the Kurganov-Petrova central upwind scheme is better than the one with staggered grid, mainly because the results are less noisy -noisy, inaccurate results are not good for control design.However, there is also one more difference, apart from the noise, and this is the simulation time, which is also very important for a controller.With respect to simulation time, the (first order) staggered grid scheme is 3-4 times faster than the (second order) non-staggered grid scheme; the staggered grid method takes 3-4 sec to solve the model with N=200 grids, while the Kurganov-Petrova method takes 12-13 sec to solve the model over 200 grids.It should be noted that here, the Kurganov-Petrova central upwind scheme has been implemented utilizing array operations in MATLAB; if these are replaced by loop constructs, the simulation time increases to some 1.5 min (Sharma, 2015).
Another way to reduce the simulation time is decreasing the number of grid cells, but this would also lead to less accurate results.It is thus of interest to compare the results of simulating the model with the  10.As it is seen from figure, the simulation results with the number of cells equal to 100 (2-3 s), 200 (12-13 s) and 500 (790 s) is almost the same, while the central upwind scheme with 50 (3-4 s) cells shows a small deviation from the results with the higher number of cells.Thus, decreasing the number of cells to less than 100 leads to significantly less accurate results.Thus, choosing 100 grid cells seems to be a good compromise between accuracy and computation time.
5 Experiment and model tuning

Experimental data
The results of the experiment for the water level ahead of the Grønvollfoss dam are presented in Figure 11.As described in Vytvytskyi et al. (2015); Vytvytskyi (2015), the increase from 120 to 160 m 3 /s of inlet flow rate starts after 10 min in the figure.Then after around 15 more min, the inlet flow decreases to the original 120 m 3 /s, thus after 25 min in the figure.The water level ahead of the Grønvollfoss dam appears to start increasing at around 20.5 min, which is 10.5 min after the inlet flow rate was increased.Then oscillations are seen in the water level ahead of the dam.These oscillations describe a wave, which is initiated by the inlet flow jump, with reflections from the Grønvollfoss dam back and forth to Årlifoss.This process continues with damping until a new steady state is reached.The period time of each oscillation is around 21 min and the steady state is reached after more than 80 min.The jump of the inlet flow rate leads to an increase in the water level ahead of the dam of around 3 cm.From Figure 11, it should be noted that there is a decreasing trend in the oscillations (the average of maximum and minimum values of one oscillation becomes lower for every next oscillation).This trend can be caused by an outlet flow rate which is slightly larger than the inlet flow rate to the system, e.g.Vin = 120 m 3 /s but Vout = 122 m 3 /s.Again, remember that the flow rates are only given indirectly and somewhat inaccurately via the power production.

Manual tuning of the model
After choosing the Kurganov-Petrova central upwind scheme as the chosen method for solving the Saint Venant model of the run-of-river hydropower system, it is of interest to fit this model to the experimental results.Using the same parameters in Table 2 and the number of cells equals to 100, the central upwind scheme is compared with the experimental data in order to observe differences and make some assumptions for the tuning; see Figure 12.The figure shows that the simulation results differ in the value of the water level as well as the period of oscillations before the new steady state is reached.In addition, the experimental results have decreasing trend caused by the difference in the inlet and outlet volumetric flow rates (the outlet flow rate is slightly bigger after the step).This difference in the volumetric flow rates leads to less increase of the water level than in the simulations.The previous work (Vytvytskyi et al., 2015;Vytvytskyi, 2015) shows that varying the length of the river can be used to fit the period of the oscillations.However, the length of the river is probably one of the most accurate parameters of the river geometry and it is somewhat undesirable to tune this parameter.Thus, here we attempt to fit the oscillation period by varying other parameters such as the slope geometry of the river bed and initial values of the water level ahead of the Grønvollfoss plant: these parameters are less reliable than the river length.
After some simulation experiments, it turns out that the position of the steeper part of the river bed (middle part, see Figure 2) together with the height drop of the river can significantly affect the period of oscillations.An example of this is shown in Figure 13 and Figure 14, where Figure 13 shows two different river bed postulates together with the initial water levels, while Figure 14 compares the water level ahead of the dam for these two cases of the river bed geometry.Figure 13 shows that for the new river bed geometry, the steeper (middle) part is shifted 0.3 km to the right (downstream, closer to the Grønvollfoss plant).This means that now the length of the first part of river bed is 2.8 km, the length of the middle part is still 0.5 km and the length of the third part is 1.7 km.There is also small decrease in the height drop, from 17.5 m to 16.7 m.It should be noted that the initial water level ahead of the dam in previous simulation was 17.1 m, which is 0.4 m less than the height drop, but for the new river bed this difference is reduced to 0.1 m, which means that the initial water level ahead of the dam is 16.6 m.The result of these changes of the river bed geometry is the increasing of the period of the oscillation, as shown in Figure 14.After this manual tuning of the river bed slope geometry and the initial level, it is of interest to tune other model parameters in order to improve the model fit.The friction factor and river width are chosen as tuning parameters.The results of a first, manual tuning is shown in Figure 15.Here we use the new river bed described above, together with the Strickler friction factor equal to 32 m 1/3 /s and a river width equal to 195 m, as well as an outlet volumetric flow equal to 122.3 m3/s.As seen from Figure 15, the simulation results fit the experimental data quite well.

Least squares model fitting
It is of interest to see if the manually tuned parameters can be improved further using a least squares error method.Using the manually fitted results as initial values, this method can find optimal values for the friction term, the river width and the height drop based on minimizing some squared error.The least squares error method is easily implemented in MATLAB using the functions for optimization problems.As variables that need to tuned, we will use the friction factor, the river width and the height drop of the river slope.The modified geometry of the river bed (the length of each parts) and the outlet volumetric flow rate are maintained as in the last simulations above.
The results of the least squares model fitting is shown in Figure 16.The following values of the optimized parameters are found: the Strickler factor is approximately 31 m 1/3 /s, the river width is approximately 194 m and the height drop is approximately 16.7 m.

Conclusions
In this paper, the use of the Saint Venant equations for a run-of-river system has been explored.The Finite Volume method with staggered grid was compared with the Kurganov-Petrova central upwind scheme and the latter was determined as better for discretizing the Saint Venant equations; the Kurganov-Petrova algorithm is a high resolution scheme with a less noisy result.There is also a disadvantage with this central upwind scheme: it leads to longer simulation time than the staggered grid method; this problem can be reduced by reducing the number of cells to, say 100, and by using the ode23 solver.
An experiment on the system was carried out to enable the comparison of the model results with reality, and to tune model parameters for better model fitting.After some simulation it was found that the slope of the river bed (the height drop together with position of the steeper part of the river bed) can affect the period of oscillations in the simulations results.Changing the slope geometry is better than changing the length of the river, as done in [1, 2], since the river length is known relatively accurately.Thus, firstly this slope was manually tuned against the experimental data and after this the least squares algorithm was used to define the optimal values of the friction factor, the width and the height drop of the river.This tuned model shows quite good fitting of the model to the experimental results.
The resulting improved model of water transport in the run-of-river system presented here, can be further used for control synthesis and analysis for the run-ofriver system.

Acknowledgment
Ingunn Granstrøm from Skagerak Kraft AS, Norway, provided necessary information about the geometry and other details of the Grønvollfoss run-of-river power plant in Telemark, and helped to organize the experiment on the Grønvollfoss power plant together with a team from Skagerak Kraft AS.This is gratefully acknowledged.

Figure 3 :
Figure 3: Google map showing a part of the river Tinnelva between Årlifoss and Grønvollfoss power plants.3

Figure 12 :
Figure 12: Comparison of the experimental and central upwind scheme results.
Figure13: Two different river geometries, together with initial water levels.

Figure 16 :
Figure 16: The results of the least squares model fitting.

Table 1 :
Experiment table Inlet volumetric Outlet volumetric Time interval for , Initial level ahead flow rate, [m 3 /s] flow rate, [m 3 /s] increased flow, [min] of the dam, [m.a.s.

Table 2 :
Parameters of the river flow system.Figure 5: Zoomed plot of the water level ahead of the Grønvollfoss dam (ode23 solver).Figure 6: Zoomed plot of the volumetric flow rate in the last cell (ode23 solver).Figure 8: Zoomed plot of the water level ahead of the Grønvollfoss dam (ode23t solver).central upwind scheme with different number of cells; see Figure Figure 11: Plot of experimental result of the water level ahead of the Grønvollfoss dam.
Figure14: The difference in the water level ahead of the dam for two cases of river bed geometries.Comparison of the experimental results with the results from the central upwind scheme with modified river bed slope geometry, and manually tuned parameters for friction factor and river width.