A Discontinuous Hamiltonian Approach for Operating a Dam-Reservoir System in a River

We present a new optimal control model for designing operation rules of a dam-reservoir system in a river. The system is assumed to be optimized so that the basic operation purpose of the system and downstream aquatic ecosystem are reasonably balanced. Starting from the water balance dynamics of a reservoir, we show that finding the optimal operation policy of the system ultimately reduces to solving a Hamilton-Jacobi-Bellman (HJB) equation having a discontinuous Hamiltonian. This discontinuity comes from a physical constraint in modeling the water balance dynamics and a functional shape of the performance index to be optimized. We obtain an exact continuous viscosity solution to the HJB equation. We also present numerical schemes for discretization of the HJB equation, which are the local Lax-Friedrichs scheme and its Weighted Essentially Non-Oscillatory (WENO) counterpart. Both schemes can approximate the exact solution, the latter being more accurate and efficient. A computational example of a hypothetical dam-reservoir system receiving a seasonal inflow discharge is finally demonstrated.


Introduction
Nowadays, operating a dam-reservoir system in a river requires controlling the outflow discharge so that the balance between the operation purpose of the system and the downstream ecological conditions is concurrently optimized. The outflow discharge should be controlled so that a sufficient amount of water resources are provided, while too large or too small discharge critically affects aquatic ecosystems in the dam-downstream reaches (1) . In some cases, water environment and ecosystem in the reservoir have to be concurrently considered when designing an operation policy of a dam-reservoir system as well (2) . Optimization of a dam-reservoir system is therefore a complex engineering problem.
Optimal control theory (3) provides an effective way for dynamically optimizing dam-reservoir systems in rivers. Recently, Yoshioka and Yoshioka (4) considered a stochastic control model for operating dam-reservoir systems, and mathematically and numerically analyzed its optimality equation, which is often referred to as the Hamilton-Jacobi-Bellman (HJB) equation (3) . There exist related studies on optimal control of dam-reservoir systems focusing on industrial applications such as the hydropower generation (5)(6) .
A drawback of the research of Yoshioka and Yoshioka (4) is that the admissible set of the controls, which collects admissible rules of the outflow discharges, is somewhat artificially modified for well-posedness of the optimization problem. This modification can improve regularity of the problem, while loses some physical aspects of the water balance especially when the reservoir water volume is sufficiently large or small. Another drawback is that the numerical method they used seems to be convergent but has only low-order accuracy. Overcoming the above-mentioned two drawbacks would contribute to establishment of a mathematically rigorous as well as practical model for optimization of dam-reservoir systems in rivers. This is the motivation of this paper.
The objectives and contributions of this paper are formulation of a new simple optimal control model of a dam-reservoir system and its numerical simulation. The water balance in a reservoir is described as an ordinary differential equation having the outflow discharge as the control variable. The performance index to be optimized contains the deviation of the outflow discharge from a prescribed target value and a penalty incurred when the water volume in the reservoir is too large or too small. Solving the optimal control problem reduces to solving an HJB equation having a discontinuous Hamiltonian.
The discontinuity partly comes from the absence of the artificial modification of Yoshioka and Yoshioka (4) . Our model having the discontinuity is more complicated but is instead more physically reasonable. We present an exact non-smooth steady weak solution of the HJB equation. The discontinuity of the Hamiltonian harmonizes with the non-smoothness of the solution as a continuous viscosity solution (7) . We demonstrate that the local Lax-Friedrichs finite difference scheme, a widely-used numerical scheme (8) , and that its higher-order counterpart (9) can approximate the solution and that the latter is more accurate and efficient. A demonstrative computational example of a hypothetical dam-reservoir system receiving a seasonal inflow is finally presented. In this paper, technical details like the proofs of the propositions are omitted. Instead, brief sketches of the proofs of the statements are provided.

Water Balance Dynamics
We consider water balance dynamics of a reservoir in a dam-reservoir system. We assume that the system is created in a river, receives the inflow discharge from the upstream river reach, and releases the outflow discharge into the downstream river reach. See, Fig. 1.
The time is denoted as 0 t ³ . The water volume in the reservoir is the state variable of our model, and is denoted as t V at t . We assume that the range of the water volume is the compact set is the reservoir volume. The inflow discharge is denoted as t Q at t . Similarly, the outflow discharge is denoted as t q at t . Considering the classical mass conservation principle, the water balance of the reservoir is described by an ordinary differential equation as subject to an initial condition 0 represents the aggregated components other than the inflow an outflow discharges, such as the evaporation from the water surface, direct rainfall, and infiltration (4) . We assume that t Q and f are smooth functions of their arguments for the sake of simplicity of analysis.

Admissible Range of the Control
We assume that the inflow discharge is uncontrollable, while the outflow discharge is controllable through operating the dam. The outflow discharge is thus considered as the control variable of the model. Its admissible range has to be carefully specified considering the physical principle (4) . That is, the reservoir volume must satisfy the condition t V ÎW for all 0 t ³ . The admissible range

( )
A A V = of the outflow discharge at the reservoir volume We assume that the upper-bound q is sufficiently large. This assumption is justified when the reservoir does not receive large inflows like those observed at severe floods. Based on the physical consideration, we must have Then, we can set the admissible range A at the boundary points 0, Consequently, the admissible range ( ) A V is compact for each V Î W , and is not continuous at 0, . This discontinuity naturally arises from the physical consideration of the water balance. The past research avoided this discontinuity issue by artificially regularizing the admissible range near 0, V V = (3,10) . On the other hand, we do not use any regularizations in this paper.

Performance Index
The performance index is the index to be minimized through dynamically optimizing the outflow discharge. We assume that the operator of the dam wants to control the outflow discharge so that it does not deviate from the prescribed possibly time-dependent target value [0, ] t q q Î . In addition, we assume that the operator wants to place the water volume in a prescribed possibly time-dependent range [ , ] where 0 T > is the terminal time, 0 d > is the discount rate as the inverse of the characteristic time scale of the operation preference of the operator, 0 w > is the weight parameter, and c is the lower-semicontinuous step function given by In the performance index (7), the first term measures the deviation between the targeted and controlled outflow discharges. The second term penalizes the water balance dynamics such that the water volume does not belong to the desired range t W . The value function : [0, ] T F´W ® R is defined as the minimized performance index with respect to the path of the out flow discharge ( ) 0 t t q ³ : Note that each path ( ) 0 t t q ³ has to be measurable and its range is given by ( ) In addition, it has to be chosen so that the system (1) admits a unique continuous solution. Our goal is to find the control path ( ) 0 t t q ³ to achieve the value function F . The optimal outflow discharge as a minimizer of the right-hand side of (9) is denoted as * t q at each 0 t ³ .

HJB Equation
The HJB equation is the degenerate parabolic equation that governs the value function F . Solving the HJB equation directly leads to the optimal outflow discharge to achieve the value function. Based on the dynamic programming principle (3) , the HJB equation is derived formally as a terminal value problem: subject to the homogenous terminal condition 0 The Hamiltonian H is convex with respect to the third argument. In addition, it is discontinuous at the boundary points The discontinuity at the boundary points comes from the discontinuity of ( ) A v at these points and therefore has a physical background. On the other hand, the discontinuity at the interior points is due to the discontinuity of the penalization in the performance index (7). Handling a HJB equation having a discontinuous Hamiltonian should be with a special care as discussed in Calder (11) . Especially, it is nontrivial whether we can solve and/or approximate solutions to such an HJB equation in a classical point-wise sense. This topic is addressed in the next section.
The maximizer of the right-hand side of (12) is denoted as , which is actually the optimal outflow discharge as a function of the observation . This is why we employ the symbols with an abuse of notations. This * q is calculated as By (10) through (13), we can find the optimal outflow discharge if we could find the value function.
Finally, notice that our HJB equation does not require any boundary conditions at Heuristically, this is explained by the fact that the drift coefficient t v q f + is non-negative at 0 v = and non-positive Namely, characteristics curves are outgoing on both the boundaries, meaning that no information comes from outside the domain W .

Exact Viscosity Solution
Finding an exact solution to a HJB equation is an interesting mathematical research topic in itself. In addition, such a solution can be utilized for checking accuracy of numerical schemes. Our HJB equation is exactly-solvable under the steady condition, which is the large T limit with time-independent coefficients: where the dependence on t has been omitted. The steady problem is valid if the environmental conditions vary much slowly than the state and control variables. It is crucial to define the meaning of solutions to the HJB equation because its solutions are non-smooth in general, despite the equation is a "differential" equation. Set the upper-and lower-semicontinuous envelopes of the Hamiltonian H with respect to the first argument as * H and * H , respectively (3) . We show that the following definition of continuous "viscosity" solutions, which are appropriate weak solutions to degenerate elliptic equations, harmonizes with our optimal control problem. Notice that the discontinuity of H is effectively considered in the definition and the possible non-smoothness of viscosity solutions are not directly incorporated into it. For more details of viscosity solutions, see the past research (3)(4) .
A function ( )

C Y Î W is a viscosity solution if it is a viscosity sub-solution as well as a viscosity super-solution.
According to Definition 1, viscosity solutions must be continuous, but need not be continuously differentiable. In this way, we can consider non-smooth solutions to the HJB equation (16). Although not discussed in this paper, viscosity solutions to the original time-dependent HJB equation (10) can also be handled in a similar manner (3) .
We present a candidate of exact viscosity solutions to the steady problem (16). Without any loss of generality, we set 1 V = under a suitable normalization. Assume the following condition on the model parameters: The first one (19) means that the target discharge equals the inflow discharge, which is the baseline operation policy in some existing dams. The second one (20) means that the inflow and outflow discharges are dominant in the water balance dynamics and there exists no effective lower-bound of the outflow discharge. The third one (21) is satisfied if the discount rate is sufficiently small (the operator has a long perspective on the dam operation) and the penalization is moderately small. Under the assumptions (19) through (21), we find the following candidate of exact viscosity solutions: It is straightforward to see that the candidate (22) is continuous on W , but is not continuously differentiable at the interior points , v W W = .
Considering Definition 1, we should check the viscosity solution property at the four points 0, , ,1 v W W = because it is non-trivial to verify the property on the boundary and at the points where the solution is non-smooth. The proposition we want to show here is the following one:

Proposition 1
The function 0 F is a viscosity solution in the sense of Definition 1.

(Sketch of the proof of Proposition 1)
At , v W W = , there exists no test function for viscosity sub-solutions. The viscosity super-solution property is directly verified with the help of the lower-semicontinuity of c . At the boundary points 0,1 v = we can again directly verify the viscosity solution property using the ordering properties of the derivatives between the viscosity sub-and super-solutions and their test functions with the help of (14) and (15). As an example of the viscosity super-solution property, we see Similar (but reversed) inequalities apply to viscosity sub-solutions. Checking the viscosity properties except for at 0, , ,1 v W W = is trivial by the smoothness of 0 F . □ Fig. 2 shows the exact solution 0 F and its associated optimal outflow discharge * q . The latter is calculated using the formula (13). The non-smoothness of both 0 F and * q is clearly visualized in the figure. We see that the solution 0 F is decreasing and increasing in [0, ] W and [ ,1] W , respectively, and vanishes in [ , ] W W . The optimal outflow discharge * q equals the target value q in [ , ] W W , while it is determined in a state-dependent manner in the other parts of W . This * q is discontinuous at the interior points , v W W = in particular. A remark on the definition of viscosity solutions is placed at the end of this sub-section. The Hamiltonian is discontinuous at the boundary points 0, v V = as discussed above. It is possible to define a "constrained" viscosity solution to the presented optimal control problem (3) . The main difference between the above-presented and this definitions of viscosity solutions is that the latter only requires the one-sided inequality at the boundary points: In (27), the admissible range of q is still given by (2) even at the boundary points. This formulation is not used in this paper, but is certainly another candidate to handle the HJB equation. In this case, the super-solution property is required only inside W . In this way, we can avoid the discontinuity but loses the equality at 0, v V = .

Numerical Computation
We could find the exact solution 0 F ; however, its applicability is limited because it applies only to the steady case under the assumption (21) of the model parameters.
We develop two numerical schemes for numerical discretization of the HJB equation (10). The first scheme is the local Lax-Friedrichs finite difference scheme as a simple and widely-used scheme for solving degenerate elliptic and hyperbolic problems (8) .
The details of the scheme are omitted here because full discretization in our case is a bit long. An important theoretical fact is that, by Oberman (12) , the scheme can approximate steady viscosity solutions to HJB equation under a sufficiently fine temporal-resolution with the classical forward-Euler scheme. This is because the scheme is monotone (discretized system inherits the degenerate elliptic property in a discrete sense), stable (numerical solutions do not diverge), and consistent (A formal convergence is guaranteed).
A drawback of the local Lax-Friedrichs scheme is that it has at most first-order accuracy, and is possibly too much diffusive. The latter implies that the non-smoothness of viscosity solutions to our HJB equation may be smeared out. The scheme with the third-order WENO reconstruction (9) is therefore developed as well. WENO reconstructions have been one of the most effective methods for sharpening numerical solutions of degenerate elliptic and hyperbolic equations. A WENO reconstruction means that the computed numerical solutions at each computational vertex is nonlinearly reconstructed considering some local smoothness indicators. This enables us to improve accuracy of the local Lax-Friedrichs scheme. Recent developments of WENO reconstructions and related methods are surveyed in Shu (13) . We used a one-sided discretization (9) at both the boundary points since WENO reconstructions use non-compact computational stencils.
The WENO reconstruction used in this paper is the third-order one developed by Jiang and Peng (8) . An (Blue) and the associated optimal control * q (Red).
advantage of this WENO reconstruction is its simplicity and high implementability. Using a higher-order WENO reconstruction is theoretically possible, but may not be so practical because of the discontinuity of the Hamiltonian H around which accuracy of numerical schemes will degrades. Nevertheless, we show that the Lax-Friedrichs scheme with the WENO reconstruction generates numerical solutions more accurate than that without it. A series of numerical experiments are carried out to see computational accuracy of the local Lax-Friedrichs scheme and its WENO counterpart. We use the exact solution 0 F for their verification. The domain W is uniformly divided into N cells with 1 N + vertices. The distance between each successive vertex is 1 / N . The time step to discretize the temporal term of (10) is set as 0.5 / N and the terminal time is set as the sufficiently large value 500 T = so that each numerical solution is sufficiently close to be time-independent at the initial time The classical forward Euler scheme is used for the temporal integration of the spatially semi-discretized equation. Table 1. summarizes the maximum error between the exact solution and the numerical solution for each N . Here, the error is defined as the maximum absolute value between the exact and numerical solutions among all the vertices. We see that both the local Lax-Friedrichs scheme and its WENO counterpart have first-order computational accuracy and that the latter has smaller errors. Fig. 3 shows that employing the WENO reconstruction can reduce artificial dissipation and more sharply captures non-smoothness of the exact solution. It is remarkable that the WENO reconstruction does not achieve third-order accuracy. This may be due to the discontinuity of the Hamiltonian, which potentially degrades regularity of solutions and thus accuracy of numerical schemes. Using a locally stretched computational mesh around the discontinuity may resolve this issue. This is beyond the scope of this paper and will be addressed in future.
Finally, the computational time of the local Lax-Friedrichs scheme is almost equal to or larger than the half of that of the WENO counterpart. Under a fixed set of model parameters, assume that we uniformly divide the domain W into 1 N + vertices as in the computational example above. Then, because of using an explicit temporal integration method, we need ( ) O N temporal steps due to the Courant-Friedrichs-Lewy condition. Therefore, using the local Lax-Friedrichs scheme requires approximately 2 CN operations with 0 C > a positive constant. The computational example implies that the WENO reconstruction, despite it is more computationally expensive, requires only approximately ( ) 2 2 2 / 2 0.5 C N CN = operations to achieve the same level of accuracy. This estimate implies that the WENO reconstruction, despite it is not truly third-order accurate in the presented case, is more efficient than the original local Lax-Friedrichs scheme.

Demonstrative Application Example
The presented mathematical model is finally applied to finding the optimal operation policy of hypothetical dam-reservoir system and environmental conditions. We set the following parameter values that are time-independent: (m 3 /s). We assume the relationship (19). We set (4)  . In addition, we set 50 w = . Finally, we set the time-periodic inflow discharge assuming two peaks in each year (summer and winter) as follows: where the unit of the time t in (29) is (day). Under this problem setting, the operator wants to have the same inflow and outflow discharges, but there exists a positive minimum value of the admissible range of the outflow discharge. This minimum value can be considered as the minimum flow discharge serving as an environmental flow to maintain functions of the downstream aquatic ecosystems (14) . By (29), the inequality t Q q < holds true in some time intervals. Therefore, at least the trivial control policy t t q Q = is not admissible in this hypothetical case.
The local Lax-Friedrichs scheme with the WENO reconstruction is utilized in the numerical simulation. The temporal integration is again carried out with the forward Euler scheme. The domain W is uniformly divided with 401 vertices and the time increment for the temporal integration is set as 0.01 (day). Under the above-specified computational setting, the numerical solution is sufficiently close to be time-periodic in the time interval [0,365] (day). Here, recall that we start the computation from the terminal time t T = to the initial time 0 t = . Figs. 4 and 5 show the computed value function F and the optimal outflow discharge Q (m 3 /s) in the time interval [ ] 0, 365 (day), demonstrating a time-periodic state-dependent optimal dam operation policy. The operator of the dam-reservoir system can find the optimal outflow discharge * q by referring to each observed ( ) Fig.  5. The computational results demonstrate that the employed numerical scheme can handle our HJB equation in a stable manner. The presented numerical simulation is under a simplified condition, but its applications to more realistic problems would not encounter severe technical difficulties.

Conclusions
An optimal control model for designing balanced operation rules of a dam-reservoir system in a river was established. The associated HJB equation having the discontinuous Hamiltonian was derived and its exact viscosity solution was found. The local Lax-Friedrichs scheme and its WENO counterpart were presented and were applied to discretization of the HJB equation. The computational results suggested that both schemes can deal with the HJB equation in a stable manner and the WENO counterpart is more accurate and efficient. An application example of the model to the hypothetical dam-reservoir system demonstrated its potential applicability to real problems arising in environmental engineering.
There remain many issues to be addressed. First issue is considering the stochastic nature of the inflow discharge, which was assumed to be deterministic in this paper for the sake of simplicity of analysis. There exist at least two options to resolve this issue. The first one is using the regime-switching modeling that can handle high and low inflow regimes in the simplest case (15) . The other is using some stochastic differential equation governing the inflow discharge 16) . The second issue is to consider dynamics of aquatic ecosystems in the dam-downstream. This will  require adding some population dynamics of aquatic species to the model (17) . In any case, the total number of the state variables of the model increases. The other issue would be how to consider river environmental management based on the presented model. In real problems, several stakeholders, such as a dam operator and fishery cooperatives, are involved in river environmental management. Such a problem can be formulated as a differential game among many players. From a computational standpoint, we can examine the other higher-order schemes for solving the HJB equation (18) . Currently, we are planning to extend the presented model to address the above-mentioned remaining issues.