Model based sustainable management of regional water supply systems

Sustainable management of water resources and a safe supply of drinking water will play a key role for the development of the human prosperity in the following decades. This paper presents an optimal control approach for the management of the total water resources in a fast developing region under a critical water shortage. A model of the water allocation system is developed which considers both surface and groundwater resources and the distribution network. The spatially distributed groundwater model is considered as a reduced model of the complex 3D Finite Element model. Hence the optimization problem, which is formulated as large scale structured non-linear programming problem, can be solved in an appropriate computation time. The performance of the proposed concept is demonstrated by close to reality optimization scenarios.


INTRODUCTION
The sustainable management of the water resources and a safe supply of drinking water will play a key role for the development of the human prosperity in the following decades.The fast growth of many cities puts a large pressure on the local water resources, especially in regions with arid or semiarid climate.A current research project at Fraunhofer IOSB and AST aims to investigate ways for economic and sustainable use of the available water resources in the region of the capital of China, Beijing [12].
A main issue of the project is to develop components for a model-based decision support system (DSS), which will assist the local water authority in management, maintenance and extension of the water supply system at hand [11].This paper deals with the derivation of suitable management strategies for a mid (till long) term horizon based on assumptions for future environmental and socio-economic conditions, which are provided by other modules of the DSS.The general structure of the proposed optimal control DSS module is shown in Fig. 1.An overview about several DSS concepts and implementations is provided in [5,6].A challenge of the given problem is the large area of the water supply system.The water management has to consider the total water resources of five river basins with an area of 16,800 km² as well as a large groundwater storage in the plain with an area of 6,300 km².The main portion of the annual precipitation (85 %) in this semi-arid region is falling from June to September leading to a highly uneven distribution throughout the year.The formerly abundant groundwater resources have been overexploited over the last decades resulting in a strong decline of the groundwater head (up to 40 m).Five reservoirs are important for the management of the surface water in the considered area, where the two largest account for about 90 % of the total storage capacity of roughly m³.The water is distributed to the customers using rivers and artificial transport ways (channels, pipes) of a total length of about 400 km.
A common approach for policy generation in this field of application is to use mathematical programming techniques based on a dynamic model of the essential elements of the water allocation and distribution system [5].The great impact of the groundwater storage for the supply system at hand requires a more detailed description of the groundwater flow dynamics compared to other known DSS implementations.Therefore a 3D Finite-Element model of the plain region has been developed.However, a direct integration of this 3D model into an optimal control framework is not possible due to its computational costs.A trajectory based model reduction scheme is proposed, which guarantees a very fast response of the DSS in combination with a specially tailored non-linear programming algorithm.
The paper is organized as follows: In section 2 the essential models (surface and groundwater models) and the model reduction approach are described.While the formulation of the optimal control problem is subject of section 3, the numerical solution of the large scale structured non-linear programming problem is described in section 4. First results of the optimal water management approach are presented in section 5.

II. WATER ALLOCATION MODEL
The water allocation model can be divided into the surface water model and the groundwater model.The parts of the water allocation model are described in the sequel.

A. Surface water model
The surface water model has to comprise all important elements for the allocation, storage and distribution of water within the considered region.The intended field of application of the decision support system under development embraces management and upgrading strategies for the mid and long term range.This implies a significant simplification for the process model, because the retention time along the different transport elements (river reaches, channel or pipelines) is less than the desired minimum time step for decision of one day.Therefore a simple static approach for flow processes is sufficient and the use of sophisticated models for the dynamics of wave propagation (like e. g.Saint-Venant-Equations) with respect to control decisions is avoided.The surface water system is described as a directed graph.The edges characterize the transport elements and introduce only a variable for the discharge.The nodes represent reservoirs, lakes, points of water supply or extraction and simple junction points.Every node constitutes a balance equation involving the edges linked with and possibly the storage volume.The sole nonlinearity results from modeling the evaporation from the water surface of the storages (volume-area-curve), which is described by a piecewise polynomial approach.Channels with a very low slope are modeled as water storage.The level dependent upper bound for the channel outflow is derived from a steady state level-flow relation like e.g.Chezy-Manning friction formula and is directly added as constraint to the optimization problem.At the moment the surface water model consists of 53 nodes and 57 edges.

B. Groundwater water model
The second important water resource is the groundwater in the considered area that is modeled by a spatially distributed finite element groundwater model.The governing equation for groundwater flow is Darcy's law [7] describing slow streams through unconfined aquifers.It denotes h [L] the hydraulic head (which corresponds to the groundwater level) and kf [L/T] the hydraulic conductivity that governs the hydrogeological properties of the soil.Combining Darcy's law with mass conservation yields the partial differential equation (1) which is a diffusion equation.

 
In (1) S0 [1/L] denotes the specific storage coefficient.The terms on the right hand side of (1) summarize all sources and sinks that coincide with the time dependent groundwater exploitation due to industry, households and agriculture (Qexpl [1/T]) and recharge e. g. due to precipitation and irrigation (Qrech [1/T]) in .
(1) is an initial-boundary value problem which has to be solved numerically for h in the 3 dimensional model domain .The groundwater model ( 1) is implemented within the FEFLOW, which is a Finite Element (FEM) software [4] specialized on subsurface flow.The initial condition is h (,t0) (groundwater surface) at the initial time t0.The inflow/outflow is described by Dirichlet boundary conditions, i.e. h () at the boundary  and by well boundary conditions, that define a particular volume rate into or out of .The advantage of the latter one is that they are scalable.The 3D FEM model consists of more than 150,000 nodes (cf.Fig. 2), distributed on 25 layers.Huge computational costs result from this high resolution.The simulation of 5 years needs ~15 Minutes on an Intel Core 2 Duo CPU (2.5 GHz).This is the motivation for model reduction (see next subsection).The main task with respect to the groundwater model is the parameterization of the large-scaled model covering an area of 6,300 km².On the one hand the time independent soil parameters k f , S 0 have to be estimated and generalized for the whole domain  by a (small) set of measured values.On the other hand the source / sink terms Q expl , Q rech have to be calculated time dependent.For these calculations time dependent maps of precipitation and water demand are needed.The water demand is splitted into the three user groups households, industry and agriculture (see [8] for details).This parameterization issue is supported by powerful geographical information systems (GIS).

C. Reduced groundwater water model
For the optimization of the water allocation system the full 3D Finite Element model (with > 150,000 nodes) is not very suited due to the mentioned large computational time.As for the optimization task a prediction of the hydraulic head (groundwater level) at a set of representative and fixed points is sufficient, an input-output model (e.g. a linear state space model) with considerably smaller order n (e.g.n < 50) than the original FEM model has to be derived.
Methods for model reduction of those large scale systems have gained increasing importance in the last few years [9].Two main classes of methods for model reduction can be identified, namely methods based on singular value decomposition (SVD) and Krylov based methods.SVD based methods are suited for linear systems and nonlinear systems of an order n < 500 (e.g. balanced truncation for linear systems, proper orthogonal decomposition (POD) for nonlinear systems).Most of these methods have favourable properties like global error bounds and preservation of stability [9].Krylov based methods are numerically very efficient as only matrix multiplications and no matrix factorization or inversion are needed.Hence they are also suited for large-scale systems.Unfortunately, global error bounds and preservation of stability cannot be guaranteed.Hence actual research is focused on the development of concepts which combine elements of SVD and Krylov based methods [9].
All of these approaches have in common that they aim to approximate the state vector x with respect to a performance criterion, e. g. minimize the deviation between original system and reduced system for a given test input.As a black box input-output model would be sufficient for our purposes there is no need to approximate the whole state vector x.Furthermore, the dimension n of a reduced model which approximates the whole state space vector x would be in most cases n > 100.With this dimension, for the given optimization problem the solution time would be unacceptably high (~ hours).Last but not least the use of the commercial software FEFLOW also prevents the application of e. g. a Krylov based method as no model representation (e. g. state space model) is provided by the software.
The basic idea of the proposed model reduction method is sketched in Fig. 3.We assume the existence of a reference scenario which means that the time dependent input parameters uref(t) of the FEM groundwater model (especially groundwater exploitation Qexpl and recharge Qrech) are determined for the whole optimization horizon.In practical cases these reference scenarios are mostly available or can be generated by plausible assumptions.Hence the task consists in the derivation of a model which approximates the behavior of the full FEM model in the case that the input parameters u differ from uref(t).This model is gained by identification techniques: Test signals (e.g.steps) are added to the reference input uref(t) (dimension p) and the corresponding deviations from the reference output yref(t) (dimension q) are identified.Doing this separately for every component of the input-/output vectors u and y, we finally merge the (pq) single input-single output (SISO) models to a multi input-multi output (MIMO) model.For the groundwater model, the input parameters are e.g. cumulated (e.g.spatially integrated) exploitation of certain regions or cumulated exploitations of large well fields.The output parameters of the groundwater model are the hydraulic head at representative points ("observation wells").In our application 13 input and 13 output parameters have been defined by the users: The inputs consist by 9 counties and 4 wellfields, the 13 output parameters are 12 observation wells and the mean hydraulic head of the whole area of the water supply system.As the slow stream groundwater flow can be interpreted as diffusion process (cf.( 1)) only nearby located input and output parameters (e.g.regions/wellfields and the correspon¬ding observation wells) have some correlation and a SISO model with these input-/output combinations can be gained.Due to this physical reason the number of relevant SISO models is relatively small and hence the resulting MIMO model of relatively low dimension (n < 50) which is appropriate for the optimization problem.This proposed approach can be called trajectory and identification based model reduction.A similar approach (for a nonlinear large scale system) is found in [10].

III. WATER RESOURCES MANAGEMENT AS OPTIMAL CONTROL PROBLEM
The water resources allocation problem is formulated as a discrete-time optimal control problem: The state variables x are the volume content of the reservoirs and channels with small slope and the states of the reduced groundwater model.Control variables u are the discharge of the transport elements as well as the water demand of the customers.The uncontrollable inputs z are the direct precipitation and the potential evaporation for the reservoirs and the flow of rivers entering the considered region, which is derived by means of rainfall-runoff-models.The number of time steps within the optimization horizon is denoted by K .
The process equations (3) consist of the balance equations of the storage nodes and the reduced groundwater model.The balance equations of the non-storage nodes are formulated as general equality constraints (4).The objective function (2) contains the goals of the water management, which are primarily the fulfillment of the customer demand, the compliance with targets for the reservoir and groundwater storage volume and the delivery of water with respect to environmental purposes.Therefore quadratic terms are formulated, which penalize the deviations from desired values, like e.g. for the demand deficit of the demand node j: is the demand of and j u is the discharge delivered to the customer.While this term applies for every time step within the optimization horizon, other terms are formulated only for the final point of the horizon, like e.g. for the desired volume content of the reservoirs.
The inequality constraints (6) follow from the technical capabilities of the water distribution system and rules to guarantee safe operation, which are simple bounds for the control variables: as well as constraints for the reservoir volume v x : and the hydraulic head hydr h of the observation wells: With respect to the practical applicability selected parts of the inequality constraints (6) can be relaxed in order to avoid infeasible optimization problem with respect to unrealistic management demands.

PROBLEM
The optimal control problem is numerically solved as large scale structured non-linear programming problem (NLP): The optimization variables are the state and control variables of the several stages in time: The state equations of the process model are incorporated as equality constraints: The advantage of this problem formulation with optimization variables is the special sparsity structure with a block-diagonal Hessian-matrix and blockbanded Jacobian matrices ( K -number of time steps, n - number of state variables, m -number of control variables), which follows from the fact, that the process equations as coupling element between adjusting stages are linear in the state variables and the large amount of control variables (number of edges in the network description) this approach is not promising for this special field of application.
For the numerical solution of large scale non-linear programming problems interior point (IP) solver have become popular during the last years because of their superior behavior for NLPs with many inequality constraints.In this approach the objective function is expanded by adding barrier terms for the inequality constraints: The solution of the original NLP (11) results from the subsequent solution of ( 14) with a decaying sequence of 0   .The identification of right active set with its combinatorial complexity is avoided.The state of the art nonlinear interior point solver IPOPT is used for the application at hand [1].The interface for multistage optimal control problems of the optimization solver HQP [2], which provides an efficient way for problem formulation along with routines for the derivation of the g h    , , J by means of automatic differentiation [3], is used and coupled to IPOPT.
A typical water management problem (horizon of five years, discretization of one month) has about 8000 optimization variables, 5500 equality constraints and 7200 inequality constraints.The numerical solution takes approximately 60 iterations and a calculation time of 10 seconds on an Intel Core 2 Duo CPU (2.5 GHz).Table 1 shows the linear dependency of the numerical solution effort from the number of time steps within the optimization horizon.

APPROACH
The proposed concept for optimal water management is evaluated for three scenarios based on the same set of input data, which are a combination of historical rainfall measurements and customer demands reflecting the predicted development of population size and economic growth.The objective function contains four quadratic terms in order to penalize deviations from the desired final water level of largest reservoir in the system (Miyun reservoir), from the desired final average hydraulic head of the groundwater storage as well as the deficit of the customer demand separated into two groups for household/industry and agriculture.The deficit of the delivered water must be less than 5 % for domestic/industrial clients and less than 25 % for agricultural clients.For the third scenario it is assumed, that water can be transferred to the considered region up to an annual amount of 300 Mio.m³ starting from the third year within the optimization horizon.The overall water demand exceeds noticeable the natural sources.The reduced groundwater model has been derived under the assumption that this over-consumption was covered by the groundwater storage.This leads to a strong reduction of the average hydraulic head of the groundwater storage of about 9 meters during 5 years (scenario 1, see Fig. 5).Using the initial state of the Miyun reservoir and the final value of the reference trajectory for the average groundwater head as target in the objective function, for the base scenario (scenario 1) there are only small deviations from these values in combination with a minor demand deficit observable.In the second scenario the increase of the target value for the final groundwater hydraulic head at 5 m and a corresponding shift of penalty coefficients in order to keep this value lead to a better spreading of the overdraft over the different storages as well as to the customers.m³ the Miyun reservoir plays an important role for the long term management of the overall system.As can be seen, the years with aboveaverage precipitation produce only a medium rise of the water level.The second scenario, which attempts to reduce the decline of the average hydraulic head of the groundwater storage, results in an increased release of this reservoir compared to the base scenario, which corresponds to a change of the final water level from 142 m above sea level to 137 m.In the third scenario a part of this release is replaced by water from outside of the considered region, which reduces the water level decrease at about 2 m.The second reservoir is situated at a channel from the Miyun-reservoir to the city of Beijing and serves only as intermediate storage.The admissible range for water management is completely utilized by the optimal control approach.The shift in the management target causes a different operating strategy because water from the connected channel is taken to replace groundwater abstractions in the nearby regions.The change of the management target for the second and third scenario induce also an increase of the overall demand deficit from 0.1 % to up to 3.5 % (second scenario), which is nearly Fig. 5 shows the time plot of the mean groundwater level of the optimization scenarios.It is obvious that in scenarios 2 and 3 the aimed increase of the target value for the final groundwater hydraulic head at 5 m is nearly achieved.In Fig. 6 the impact of the different strategies to the exemplary input 1 (exploitation in a certain region) and exemplary output 5 (groundwater level at a defined observation point) can be studied.The exploitation is clearly decreased in scenarios 2 and 3 which correspond to an increase of the groundwater level at the observation point.Finally, from Fig. 7 it can be seen that the performance of the drastically reduced groundwater model is good, reflecting the fact that the original FEM model with more than 100,000 nodes has been reduced to a state space model with 36 states.

VI. CONCLUSIONS
In this paper an optimal control approach as a component of a decision support system (DSS) for the management of the total water resources (surface water, groundwater, external water resources) in a fast developing region under a critical water shortage has been presented.It has been proven that in spite of the large area which has to be managed and the corresponding complex surface and groundwater models the optimization problem could be solved in an appropriate computation time (~ minutes).This could be achieved by a drastic reduction of the complex groundwater model to a state space model of relatively low dimension (n < 50).The user (i.e.water allocation decision maker) is enabled to select from a number of predefined performance criteria as well as to assign constraints to the elements of the water allocation system in order to specify the management targets according to his/her needs.The performance of the proposed concept is demonstrated by close to reality optimization scenarios, whereby the benefit of a new strategic channel has been investigated with a planning horizon of 5 years.Actually the developed DSS component is used in a first version by the decision makers.Future work will focus on the application and adaptation of the developed concept and software for the water resources management of further regions with critical water shortage.

Figure 1 .
Figure 1.General Structure of the proposed optimal control DSS module.

Figure 2 .
Figure 2. Mesh of the 3D Finite Element groundwater model.

Figure 3 .
Figure 3. Reduced groundwater model as a linear state space model in combination with a pre-simulated reference scenario.
consists of eliminating the state variables from the optimization problem.The reduced dimension of the according non-linear programming problem with   Km optimization variables comes along with a loss of structure.The Hessian and Jacobian matrices are full.Because of the solution effort of

Figure 4 .
Figure 4. Water level of the largest reservoir in the water distribution system (Miyun-reservoir) as well as of a reservoir with seasonal storage capacity (Huairou-reservoir).

Figure 5 .
Figure 5.Time plot of the mean groundwater level in the optimization scenarios 1 -3.

Fig. 4
Fig. 4 shows the course of the water level for two reservoirs.Because of its capacity of

Figure 6 .
Figure 6.Time plot of input 1 (exploitation in a certain region) and output 5(groundwater level at a defined point) in the optimization scenarios 1 -3.

Figure 7 .
Figure 7. Time plot of output 5 of the full FEM model compared with the reduced model (optimization scenario 1).

Table 1 :
Numerical solution effort in dependence on the optimization horizon.