Flood Control with Model Predictive Control for River Systems with Water Reservoirs

Flood Control with Model Predictive Control for River

Systems with Water Reservoirs

Maarten Breckpot

; Oscar Mauricio Agudelo

; and Bart De Moor

Abstract: Many control strategies can be found in the literature for controlling river systems. One of these methods is model predictive control (MPC), and it has already shown its efficiency for set-point control of reaches and irrigation channels. This paper shows that MPC can also be used for flood control of river systems. The proposed controllers use the buffer capacity of water reservoirs in an optimal way when there is a risk of flooding, and they recover the used buffer capacity as fast as possible. The performance of the controllers is tested on a river system consisting of multiple channels, gates, and a water reservoir. One controller is used in combination with a Kalman filter, which estimates all the states of the river system on the basis of a very limited number of measured water levels. It was observed that the influence of this estimator on the control performance was minimal. DOI:10.1061/(ASCE)IR.1943-4774.0000577. © 2013 American Society of Civil Engineers.

CE Database subject headings: Floods; Kalman filters; Open channel flow; Optimization; River systems; Reservoirs. Author keywords: Control systems; Feedback control; Floods; Kalman filters; Open-channel flow; Optimization; River systems.

Introduction

Many different types of control strategies can be found in the liter-ature for controlling channels, such as proportional-integral, heu-ristic, predictive, and optimal controllers (Malaterre et al. 1998;

Burt et al. 1998;Litrico et al. 2006;van Overloop et al. 2005). The control strategy used in this work is model predictive control (MPC) (Rossiter 2003;Mayne et al. 2000). Model predictive con-trol originates from the process industry and is used in various applications from chemicals and food processing to automotive and aerospace applications (Qin and Badgwell 2003). Because MPC formulates the control problem as an optimization problem, it can minimize the deviations of water levels from their targets, and it can be used for set-point control. It can also be used for flood control because flood limits can be incorporated inside the optimi-zation problem. To achieve a high control performance, the control-ler requires an accurate model of the river system to be controlled. The standard way to model the dynamics of a channel is done through the Saint-Venant equations. Based on these equations, to-gether with the dynamics of hydraulic structures and junctions, a mathematical model can be derived for river systems. Because of the complexity of these models, controllers do not work directly with these equations but use approximate models. Models based

on identification techniques or simple integrator-delay models form a good approximation and have been used in combination with MPC for set-point control (Schuurmans 1997; van Overloop 2006;Wahlin and Clemmens 2006;van Overloop et al. 2010;Puig et al. 2009;Negenborn et al. 2009;van Overloop et al. 2008). How-ever, because these models describe the water levels of every chan-nel at only one point, they are not the right choice for flood control. Therefore, a model approximating the Saint-Venant equations along the entire river is used in this study. A similar model has al-ready been used in combination with MPC (Xu et al. 2011), but this work looked at only one channel and did not focus on flood control. Other works can be found on flood control, but they work with models describing the dynamics at a very limited number of points (Romera et al. 2011), do not take the nonlinearities of gates and hydraulic structures into account (Thai 2005), and none of them have water reservoirs, which should be used optimally.

Previous works have used a simplified conceptual model for the controller (Barjas Blanco et al. 2010;Breckpot et al. 2010). Re-cently, an approximate model of the Saint-Venant equations was used with a very fine spatial discretization for controlling a single channel (Breckpot et al. 2012a) and a two-channel system with one gate (Breckpot et al. 2012b). In this work, this method is extended to a larger river system where also junctions and a water reservoir are present. Also, a Kalman filter (Kalman 1960) is included in the control scheme because in practice only a limited number of water levels are measured, and MPC requires knowing the entire state of the system.

Contribution

This paper shows how MPC and a Kalman filter can be used for flood control of river systems containing water reservoirs. The authors are not aware of any other publication in which these tech-niques are applied to optimally use the buffer capacity of reservoirs for flood control on the basis of models approximating the river dynamics along the entire river. Because the proposed controllers work directly with approximations of the Saint-Venant equations describing the water-level dynamics along every channel, they can effectively be used for flood control.

1_{Ph.D. Fellow, Research Foundation—Flanders (FWO), Dept. of} Electrical Engineering-ESAT (SCD-SISTA), KU Leuven, Kasteelpark Arenberg 10, B-3001 Leuven, Belgium (corresponding author). E-mail: maarten.breckpot@esat.kuleuven.be

2_{Postdoctoral Researcher, Dept. of Electrical Engineering-ESAT} (SCD-SISTA), KU Leuven, Kasteelpark Arenberg 10, B-3001 Leuven, Belgium. E-mail: mauricio.agudelo@esat.kuleuven.be

3_{Full Professor, Dept. of Electrical Engineering-ESAT (SCD-SISTA)/} iMinds Future Health Dept. KU Leuven, Kasteelpark Arenberg 10, B-3001 Leuven, Belgium. E-mail: bart.demoor@esat.kuleuven.be

Note. This manuscript was submitted on June 20, 2012; approved on January 14, 2013; published online on January 16, 2013. Discussion period open until December 1, 2013; separate discussions must be submitted for individual papers. This paper is part of the Journal of Irrigation and Drainage Engineering, Vol. 139, No. 7, July 1, 2013. © ASCE, ISSN 0733-9437/2013/7-532-541/$25.00.

The main contribution of this paper is the formulation of the optimization problem. Inspired by a paper of Malaterre and Baume (1999), the authors did not work with the gate positions as control variables but with the gate discharges. This means that the gate equations are pulled out of the optimization problem. Because the gate dynamics are the main nonlinearities of river systems, the resulting optimization problem is less complex to solve. After find-ing the optimal gate discharges, a conversion is made based on the gate equations. The remaining model equations present in the op-timization problem are based on the Saint-Venant equations. The paper shows that, for the purpose of controller design, these equa-tions can be approximated accurately with one linear state-space model without losing control performance. The big advantage of this is that the total computation time needed by the controller is limited.

Another contribution lies in the formulation of the optimization problem such that the buffer capacity is used in an optimal way and the modification of its parameters to recover the used buffer capacity quickly.

Paper Outline

The paper starts by discussing the river system and the correspond-ing control problem, followed by showcorrespond-ing how the dynamics of a river system can be modelled. Then it is explained how MPC can be used for flood control as well as the use of the Kalman filter for estimating the states. This is followed by a discussion of the sim-ulation results when the proposed controllers are used for the test system. The paper ends with some concluding remarks.

Notation

Matrices are denoted with bold capital letters, e.g., X, whereas bold lowercase letters are used for vectors, e.g., x. Scalars or entries of vectors are not bold. The ith component of x is xi.kxk2Wdenotes xT_{Wx. Superscript (i) indicates the channel or the gate the variable} or parameter belongs to, and1n represents a vector of ones of di-mension n.

Problem Formulation

The goal of this study is to find a controller that can be used for set-point control of rivers while preventing them from flooding. For this the controller needs to use the available buffer capacity of the present water reservoirs in an optimal way during periods of heavy rainfall. The controller should also empty the reservoirs as fast as possible before focusing on set-point control again. Each channel has a safety limit and a flood limit. When the water levels risk vio-lating only their safety limits, the controller is allowed to use the buffer capacity of the water reservoirs up to its own safety limit. If the water levels risk violating also the flood limits, then all the available buffer capacity of the water reservoirs can be used. The controller should also take into account that only a very limited number of water levels are measured in practice. The controller can use the future disturbances to find the best control actions.

Fig.1gives a schematic overview of the river system used to test the control performance. Channels 1–4 form the main part of the river, whereas Channels 5 and 7 connect the water reservoir (Channel 6) with the river. The squares represent hydraulic gates that can be used to control the discharges locally. QinðtÞ is the discharge entering the river, whereas QoutðtÞ corresponds with the discharge leaving the river. When there is no flooding risk, the con-troller should keep the water levels of the first and the fourth chan-nels as close as possible to their set points. The discharge QinðtÞ

entering the river system is a disturbance signal. The controller can control the water levels with QoutðtÞ and the three gates. Every con-trol variable has an associated upper and lower limits, and the gates have a maximal rate of change constraint.

Modeling

Channel Modeling

The standard approach in the literature to model open channel flow is done by using the Saint-Venant equations (Chaudry 2008;Cunge et al. 1980;Sturm 2001;Litrico and Fromion 2009). The equations modeling the dynamics of a single channel without later inflow are the following two partial differential equations (PDEs):

∂Aðz; tÞ ∂hðz; tÞ ∂hðz; tÞ ∂t þ ∂Qðz; tÞ ∂z ¼ 0 ð1Þ ∂Qðz; tÞ ∂t þ∂z∂ Qðz; tÞ2 Aðz; tÞ þ gAðz; tÞ
_{∂hðz; tÞ} ∂z þ Sfðz; tÞ − S0  ¼ 0 ð2Þ where t = time variable; z = space variable; Qðz; tÞ = water dis-charge (m3=s); hðz; tÞ = water depth (m); Aðz; tÞ = cross-sectional flow area (m2); g = gravitational acceleration (m=s2); S₀ = bed slope; and Sfðz; tÞ = friction slope. Eq. (1) describes the conserva-tion of mass, and Eq. (2) the conservation of momentum. The resistance law Sfðz; tÞ is modeled in this study with the Manning relation (Chow 1959):

Sfðz; tÞ ¼

n2mannQðz; tÞjQðz; tÞj

Aðz; tÞ2Rðz; tÞ1=3 ð3Þ

where nmann = Manning coefficient (s=m1=3); Rðz; tÞ ¼ Aðz; tÞ= Pðz; tÞ = hydraulic radius (m); and Pðz; tÞ = wetted perimeter of the cross section (m). The channels discussed in this study have a trapezoidal shape with side slope S. The parameters of a trapezoi-dal channel can be seen in Fig. 2, where L is the length of the channel (m), and B is the bottom width (m).

Together with these PDEs, upstream and downstream boundary conditions are needed for each channel. The most simple one is where the upstream or downstream discharge is given, e.g., it

Fig. 1. River system to be prevented from flooding; arrows indicate general flow direction from upstream to downstream

Fig. 2. Parameters of trapezoidal channel

can be controlled directly with a pump or it might be a disturbance signal. Other possibilities are gates and junctions.

Gate Dynamics

The gates used in this study are underflow-vertical sluices. In gen-eral, the gate dynamics can be modeled as follows:

QgateðtÞ ¼ ~f½cðtÞ; hupðtÞ; hdownðtÞ
ð4Þ where QgateðtÞ = discharge controlled by the gate; cðtÞ = gate open-ing (m); hupðtÞ = water level on the upstream part of the gate; hdownðtÞ = water level on the downstream part; and ~f∶R3→ R is a nonlinear function. The function ~f used in this study is based on Lin et al. (2002) and Sepúlveda et al. (2009) and has the follow-ing form: QgateðtÞ ¼ CDðtÞcðtÞw ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ghupðtÞ q ð5Þ where CDðtÞ = discharge coefficient; and w = width of the gate (m). The equation used for CDðtÞ depends on the flow condition:

Free flow∶ CDðtÞ ¼ CC ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ αðtÞ p ð6Þ Submerged flow∶ CDðtÞ ¼ CC  βðtÞ − ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi βðtÞ2₋ 1 αðtÞ2− 1 ₂ 1 − 1 γðtÞ2  r ₁₌₂ 1 αðtÞ− αðtÞ ð7Þ

where CC = contraction coefficient; βðtÞ ¼ ½1=αðtÞ − 1
2þ 2½γðtÞ − 1
; αðtÞ ¼ CCcðtÞ=hupðtÞ; and γðtÞ ¼ hupðtÞ=hdownðtÞ. CC can vary from 0.598 to 0.74, but for most applications, 0.611 is usually assumed (Henderson 1966; Liggett and Cunge 1975). The gate is in free flow if hdownðtÞ is below the limit hdown;maxðtÞ:

h_down;maxðtÞ ¼CCcðtÞ 2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ_{αðtÞ½1 þ αðtÞ}
16 s ! ð8Þ Otherwise, the gate is in a submerged flow condition. The roles of hupðtÞ and hdownðtÞ interchange when hdownðtÞ becomes larger than hupðtÞ, and the discharge becomes negative.

Given these equations, the boundary conditions for two chan-nels i and j connected to each other with a gate m can be formulated as follows:

QðiÞðLðiÞ; tÞ ¼ QðjÞð0; tÞ ð9Þ if hðiÞðLðiÞ; tÞ ≥ hðjÞð0; tÞ

QðiÞðLðiÞ; tÞ ¼ ~f½cðmÞðtÞ; hðiÞðLðiÞ; tÞ; hðjÞð0; tÞ
ð10Þ else

QðiÞðLðiÞ; tÞ ¼ − ~f½cðmÞðtÞ; hðjÞð0; tÞ; hðiÞðLðiÞ; tÞ
ð11Þ This simply means that the gate discharge is equal to the discharge leaving the first channel and to that entering the second channel. Junctions

Junctions are places along the river where multiple channels coincide. At these points all the water levels of the channels should be equal, and the sum of the discharges at the end of the upstream

channels should be equal to the sum of the discharges at the begin-ning of the downstream channels. Applying this to the river system in Fig.3 results in the following equations:

hð1ÞðLð1Þ; tÞ ¼ hð2Þð0; tÞ ð12Þ hð1ÞðLð1Þ; tÞ ¼ hð3Þð0; tÞ ð13Þ Qð1ÞðLð1Þ; tÞ ¼ Qð2Þð0; tÞ þ Qð3Þð0; tÞ ð14Þ

Test Example

The river system to be controlled is shown in Fig. 1. Based on the equations and boundary conditions defined in the previous subsections, the mathematical model corresponding to this system consists of the following set of equations:

• The Saint-Venant equations [Eqs. (1) and (2)] to model the dynamics of every channel individually;

• The gate equations [Eqs. (9)–(11)] as boundary conditions between Channels 2 and 3, 5 and 6, and 6 and 7;

• Equations similar to Eqs. (12)–(14) for all channels ending in or starting from a junction: Channels 1, 2, and 5 and Channels 3, 4, and 7; and

• The following two last boundary conditions:

Qð1Þð0; tÞ ¼ QinðtÞ ð15Þ

Qð4ÞðLð4Þ; tÞ ¼ QoutðtÞ ð16Þ Table1contains the values of the parameters of the channels and the gates of the test example.

Discretization and Numerical Implementation

Because there is no analytical solution for the Saint-Venant equa-tions, the infinite dimensional variables will be approximated on a finite grid (Strelkoff and Falvey 1993). The partial derivatives are approximated with finite differences, whereas the θ-method, i.e., fðtjþ θΔtÞ ¼ θfðtjþ ΔtÞ þ ð1 − θÞfðtjÞ with θ ∈ ½0; 1
and Δt the integration step, is used for the time integration. The spatial

Fig. 3. Example of river system consisting of three channels coinciding in junction; arrows indicate general flow direction

Table 1. Parameters of Different Channels and Gates of Test Example

Parameter Channels Gates 1 2–4 5 and 7 6 1–3 L (m) 800 1,000 500 3,000 S₀ 0.0002 0.0002 0.0001 0.0001 nmann·ðs=m1=3Þ 0.02 0.02 0.02 0.02 B (m) 6 6 6 6 S 1 1 1 1 w (m) 6

grid used in this paper is a staggered grid (see Fig.4), with N water levels and Nþ 1 discharges for a given channel. The derivatives in Eq. (1) are approximated by the following [Qðzi; tjÞ ¼ Qji; and Qj_iþθ¼ θQj_iþ1þ ð1 − θÞQj_i]: ∂hj i ∂t ≃ hjþ1i − h j i Δt ð17Þ ∂Qj i ∂z ≃ Qjþθ_i − Qjþθ_i−1 Δz ð18Þ ∂Aj i ∂h ≃
∂A ∂h _jþθ i ð19Þ A similar approach is used for the terms∂Qðz; tÞ=∂t, Aðz; tÞ, ∂hðz; tÞ=∂z, and Sfðz; tÞ in Eq. (2). The advection term∂½Q2ðz; tÞ= Aðz; tÞ
=∂z is approximated with an upwinding approach:

∂ ∂z
Q2 A _j i ≃ 8 > < > : 1 Δz hQ2 A jþθ iþ1− Q2 A jþθ i i Qj_i <0 1 Δz hQ2 A jþθ i − Q2 A jþθ i−1 i Qj_i ≥ 0 ð20Þ In this way, the two PDEs describing the dynamics of a single channel are transformed into a system of nonlinear equations:

f½hðtjþ1Þ; hðtjÞ; qðtjþ1Þ; qðtjÞ
¼ 0 ð21Þ where hðtjÞ ¼ ½h1ðtjÞ; : : : ; hNðtjÞ
T; qðtjÞ ¼ ½Q1ðtjÞ; : : : ; QNþ1ðtjÞ
T; and f∶R4Nþ2→ R2Nþ1. The system of nonlinear equa-tions [Eq. (21)] for each channel, together with all the boundary conditions, needs to be solved for hðtjþ1Þ and qðtjþ1Þ of every channel given the values of these variables at time step tj. This can be done with Newton method. A discussion about the choice ofΔt and θ can be found in (Clemmens et al. 2005). In this paper, θ is set to 0.6.

Model Predictive Control

To find a controller that is able to prevent rivers from flooding, it is important that the controller considers the flood levels when look-ing for the best control actions. The controller should also take pos-sible actuator limits and rainfall predictions into account. Model predictive control is an optimization-based control strategy that can perfectly deal with these requirements. It uses a process model to predict the future process outputs within a specified prediction horizon. It determines the next inputs for the process by solving an optimization problem over this horizon, taking input and output constraints (e.g., flood limits), future disturbances (rainfall predic-tions), and the process model into account. Only the first element of the complete optimal sequence is applied to the process, the new current state of the system is measured or estimated, and the entire procedure is repeated.

Nonlinear MPC Approximate Model

Model predictive control solves at every time step an optimization problem in which the model equations are considered as equality constraints. However, the model equations discussed in the pre-vious section are too complex to be used directly inside the opti-mization problem. The main complexity comes from the nonlinear gate equations. Inspired by Malaterre and Baume (1999), these gate equations can be excluded from the optimization problem when we work with the discharge at the gates as control variables instead of the gate positions. Once the optimal gate discharges are found, a conversion step is needed to find the corresponding gate posi-tions. Under this control scheme, the boundary conditions for the river system are given by the following instead of Eqs. (9)–(11): QðiÞðLðiÞ; tÞ ¼ QðmÞgateðtÞ ð22Þ

QðjÞð0; tÞ ¼ QðmÞgateðtÞ ð23Þ The resulting model equations will be approximated with linear time-varying state-space models by linearizing the discretized Saint-Venant equations [Eq. (21)] around the linearization point xlinðkÞ, ulinðkÞ, and dlinðkÞ:

xðk þ 1Þ − xlinðk þ 1Þ ¼ AðkÞ½xðkÞ − xlinðkÞ
þ BðkÞ½uðkÞ − ulinðkÞ

þ DðkÞ½dðkÞ − dlinðkÞ
ð24Þ where xðkÞ ¼ ½hð1ÞTðtkÞ; : : : ; hðncÞTðtkÞ;qð1ÞTðtkÞ; : : : ;qðncÞTðtkÞ
T; nc = number of channels; uðkÞ ¼ ½Qð1ÞgateðtkÞ; Qgateð2ÞðtkÞ; Qð3ÞgateðtkÞ; Qð4ÞðLð4Þ; tkÞ
T; dðkÞ ¼ Qð1Þð0; tkÞ; and AðkÞ, BðkÞ, and DðkÞ = state-space matrices at time tk. Eq. (24) can be rewritten as

xðk þ 1Þ ¼ AðkÞxðkÞ þ BðkÞuðkÞ þ DðkÞdðkÞ þ βðkÞ ð25Þ where

βðkÞ ¼ xlinðk þ 1Þ − AðkÞxlinðkÞ − BðkÞulinðkÞ − DðkÞdlinðkÞ ð26Þ The linearization points xlinðkÞ, ulinðkÞ, and dlinðkÞ should be as close as possible to the values found by the optimizer based on these linear models. Therefore, dlinðkÞ is set equal to the predicted disturbance signal, and ulinðkÞ will be set equal to the optimal con-trol actions found by the optimizer in the previous time step. xlinðkÞ is found by using the nonlinear model equations to predict what the states are corresponding to ulinðkÞ and dlinðkÞ. After this prediction step, the linear models are derived and can be used in the optimi-zation problem.

Optimization Problem

The optimization problem that needs to be solved at every time step is a quadratic program (QP) of the following form:

min u;x;ξ;ζ XNP k¼1 kxðkÞ − rk2 Wþ X NP−1 k¼0 kuðkÞ − uðk − 1Þk2 R þ kξk2 Sþ sTξ þ kζk2Vþ vTζ ð27Þ s:t: xð0Þ ¼ ˆx₀ ð28Þ xðk þ 1Þ ¼ AðkÞxðkÞ þ BðkÞuðkÞ þ DðkÞdðkÞ þ βðkÞ ð29Þ

Fig. 4. Staggered grid structure used to discretize Saint-Venant equations

uðkÞ ≤ uðkÞ ≤ ¯uðkÞ; ∀ i ¼ 1; : : : ; nc∶ ð30Þ hðiÞðkÞ ≤ hðiÞmax;1þ ηðkÞξi ð31Þ hðiÞðkÞ ≤ hðiÞmax;2þ ηðkÞζi ð32Þ

ξ ≥ 0 ð33Þ

ζ ≥ 0 ð34Þ

where NP = prediction horizon; W ∈ Rnx×nx, R ∈ Rnu×nu, S ∈ Rnc×nc_{, and V ∈ R}nc×nc = four positive semidefinite diagonal

weighting matrices; nx = total number of water levels and dis-charges; nu= number of inputs; s ∈ Rncand v ∈ Rnc= two weight-ing vectors; r = set points for the states; ˆx0= current state estimate of the process; hmax;1= safety levels; hmax;2= flood levels; u and ¯u = operational limits on the inputs; ξ; ζ ∈ Rnc = two vectors of

slack variables (one slack variable for each channel); andηðkÞ ¼ 1=rk−1

c = time-dependent weight (rc> 1). The weighting matrices W, R, S, and V define the relative importance of the difference between the states and their set points, the changes of the control actions, and the two vectors of slack variablesξ and ζ. It can be shown that this QP has only one (global) solution (Nocedal and Wright 2000). In this study, quadprog of MATLAB (2011a) is used to solve the QPs. This controller will be called NMPC (nonlin-ear MPC).

The meaning of the safety limits and the flood limits was ex-plained in the “Problem Formulation” section. The flood limits and safety limits for the water levels are implemented as soft con-straints [Eqs. (31) and (32)]. The maximal violation of these limits for each channel along the prediction horizon are penalized in the objective function by the termskξk2Sþ sTξ and kζk2Vþ vTζ. A suf-ficiently large s and v will ensure that the constraints will only be violated when no feasible solution exists for the hard constrained optimization problem (Hovd and Braatz 2001). This means that the constraints are enforced as exact soft constraints. If the constraints cannot be prevented from being violated, the controller will min-imize these violations and hence reduce the flood risk. The quad-ratic termskξk2_Sand kζk2V are imposed to have a well-posed QP and are extra tuning parameters (Scokaert and Rawlings 1999). A time-dependent weightηðkÞ is used to penalize future constraint violations in the prediction horizon increasingly to avoid long-lasting constraint violations (Hovd and Braatz 2001).

Because the buffer capacity of the reservoir between the safety limit and the flood limit may not be used to keep the water levels of all the other channels below their safety limits, the diagonal element of S and the element of s corresponding to the reservoir are set higher than the elements corresponding with the other chan-nels. The diagonal elements of V and the elements of v are set higher than all the elements of S such that the controller will use the remaining buffer capacity for flood prevention of all chan-nels. If there is no flood risk, the controller needs to focus on set-point control of the water levels of Channels 1 and 4. This is achieved by choosing the elements in the matrix W corresponding to these water levels larger (but sufficiently smaller than the weights of the slack variables) than the entries associated to the water levels of the other channels. R influences the control effort of the different input variables.

uðkÞ and ¯uðkÞ are the operational limits on the inputs at time step k. These limits are time variable for the control actions corresponding to the three gate discharges. For the controllable

discharge at the end of the river, they are constant and assumed to be given. These limits are time variable because they depend on the rate of change constraintΔmaxof the gates and on the up-stream and downup-stream water levels at each gate at every time step. These limits can be calculated with an Algorithm 1. In this algo-rithm ngis equal to the total number of gates, upredis the sequence of control actions found by the optimizer in the previous time step, c_min;mis the lower limit on gate m, and c_max;mis the upper limit on gate m. The algorithm uses a function l to calculate the gate posi-tion corresponding to a desired gate discharge and upstream and downstream water levels. There is one important aspect of this function. There can be situations in which the position of a gate can remain unchanged from iteration j to jþ 1 to achieve a given gate discharge. However, for decreasing water levels, the gate open-ing can become larger than these water levels. The calculated umðj þ 2Þ and ¯umðj þ 2Þ will have the same value if the Δmax;m is not big enough to lower the gate between the maximal and min-imal water levels. This can happen for all future values of j and at the end the lower limit is equal to the upper limit for that control variable at every time step within the prediction window. This in-dicates that the value for the gate discharge is fixed and the con-troller will stop using this gate, losing in this way one degree of freedom. This uncontrollability can be avoided by limiting the maximal value returned by this function l with the maximal water level at the gate. This will make sure that umwill be different from ¯umat every time step, which will ensure the controllability of the gate. Internally, the function l solves Eqs. (5)–(7) iteratively with the bisection method.

Because the linear model used in the optimization problem is only an approximation of the nonlinear model and because of the dependence of the limits u and ¯u on the states of the system, the solution of the QP is not guaranteed to be the solution of the nonlinear optimization problem. Therefore, the prediction step, the linearization step, the algorithm, and the optimization step will be iterated over. First, the optimal control actions found in the previous time step are used with the nonlinear model to predict what the future states of the model will be within the prediction window. On the basis of this sequence of future states, a sequence of linear state space models is derived, which is used to find u and ¯u and in the optimization problem. This procedure is repeated with the new control actions until convergence or a maximum number of itera-tions is reached. Then, the controller looks for the gate posiitera-tions corresponding to the desired gate discharges and applies them to the river system. This conversion is done with the function l dis-cussed previously.

Algorithm 1. Calculation of the time-varying limits

Algorithm to find the time-varying limits uðkÞ ¼

½u1ðkÞ; : : : ; ungðkÞ

T _{and ¯uðkÞ ¼ ½¯u}

1ðkÞ; : : : ; ¯ungðkÞ

T_{. The}

function s selects the upstream and downstream water levels corresponding to gate m, whereas the function l calculates the gate position corresponding to the desired discharge for gate m.

% at time step tk ~x ¼ xðtkÞ cprev¼ ½cð1Þðtk−1Þ; : : : ; cðngÞðtk−1Þ
T for j¼ 1; : : : ; NPdo for m¼ 1; : : : ; ngdo ½hup; hdown
¼ sðˆx; mÞ if hup≥ hdown then

umðjÞ ¼ ~f½maxðcprev;m− Δmax;m; cmin;mÞ; hup; hdown
¯umðjÞ ¼ ~f½minðcprev;mþ Δmax;m; cmax;mÞ; hup; hdown
else

¯umðjÞ ¼ − ~f½maxðcprev;m− Δmax;m; cmin;mÞ; hdown; hup
umðjÞ ¼ − ~f½minðcprev;mþ Δmax;m; cmax;mÞ; hdown; hup
end if

c_prev;m¼ l½u_pred;mðjÞ; hup; hdown
end for

d¼ Qinðtkþj−1Þ ~x ¼ AðjÞ ~x þ BðjÞupredðkÞ þ DðjÞd þ βðjÞ end for

Buffer Capacity Recovery

After a period of heavy rainfall, the controller should recover the buffer capacity to handle future rainfall without flooding. This can be achieved by working with two different weighting matrices W and reference signals r. During normal operation (set-point con-trol), the diagonal elements of the matrix W corresponding to the most important water levels will have the highest values com-pared with those corresponding to the other water levels, and r is set equal to the desired set points. However, when some of the buffer capacity is used for flood prevention, the set points of the water levels of the channels downstream of the reservoir (Channels 3, 4, and 7) are set lower than the set point of the reservoir, and the diagonal elements of W corresponding to these water levels are increased. By decreasing these water levels, Gate 3 can be used to remove the excess water in the reservoir. This change in W and r is performed when one of the water levels in the reservoir is 20 cm above its set point. This is checked every time a new state estimate of the process is computed. When the buffer capacity is recovered, W and r are set equal to their initial values (set-point control). This is done when all the water levels of the reservoir are at the most 10 cm above their set point for at least 2.5 h.

Linear MPC

Because of the prediction step, the NMPC requires high computa-tion times to find optimal control accomputa-tions. This step can be avoided by working with the same linear model at every time step. To this end, the model constraint Eq. (29) in the NMPC formulation is replaced by

xðk þ 1Þ ¼ ¯AxðkÞ þ ¯BuðkÞ þ ¯DdðkÞ þ ¯β ð35Þ where ¯A, ¯B, ¯D, and ¯β are found by linearizing the discretized model equations around the nominal steady-state operating point of the river. The algorithm can still be used to find uðkÞ and ¯uðkÞ if the time-varying state-space matrices AðkÞ, BðkÞ, DðkÞ, and the vectorβðkÞ are replaced with ¯A, ¯B, ¯D, and ¯β, respectively. This controller will be called LMPC (linear MPC).

Because there is still a dependence between uðkÞ and ¯uðkÞ, the chosen control actions, and the water levels surrounding every gate, in theory several iterations involving the algorithm and solving the QP are needed. However, the simulation results show that good results can already be achieved by performing these operations only once.

State Estimator

The controllers defined in the previous section require the knowl-edge of the current state of the process x₀. However, only a very limited number of water levels are measured and none of the dis-charges. Therefore, all the other states need to be estimated. The state estimator used in this paper is the Kalman filter (Kalman 1960). This is explained here for the LMPC. The estimation of the state vector xðkÞ will be written as ˆxðkÞ. The Kalman filter uses the linear state-space model in combination with the error between

the predicted and the measured water levels to correct the state estimation through a feedback gain matrix L:

Δˆxðk þ 1Þ ¼ L½ΔyðkÞ − ΔˆyðkÞ
þ ¯AΔˆxðkÞ þ ¯BΔuðkÞ

þ ¯DΔdðkÞ ð36Þ

ΔˆyðkÞ ¼ CΔˆxðkÞ ð37Þ

where ΔxðkÞ ¼ xðkÞ − xlin; Δ ˆxðkÞ ¼ ˆxðkÞ − xlin; ΔuðkÞ ¼ uðkÞ − ulin; ΔyðkÞ ¼ yðkÞ − ylin; ΔˆyðkÞ ¼ ˆyðkÞ − ylin; yðkÞ = measured water levels; ˆyðkÞ = estimated water levels; and C = matrix that selects the measured water levels from the states. The Kalman gain L is found by minimizing the covariance of the estimation error xðkÞ − ˆxðkÞ, taking process and measurement noise into account. Both noises are assumed to be independent white noise with a normal distribution. The exact derivation for finding L can be found in Franklin et al. (1997) and Kwakernaak and Sivan (1972).

Simulation Results

This section describes the simulations results when NMPC, LMPC, and LMPC with Kalman filter are applied to the river system in Fig.1. For the first two controllers, it is assumed that all the states are measured at every time step. This allows the comparison of the control performance of these two controllers separately. Only the effect of the Kalman filter on LMPC is shown, as this controller has the lowest computation time and a similar performance as NMPC. The sampling time of every controller is 15 min, whereas the full nonlinear model of the river system generates data every minute. A spatial discretizationΔz of 50 m is used. The maximum number of iterations for the NMPC is limited to five or when the Euclidean norm of the difference of the control actions found between two consecutive iterations is smaller than 10−3. Only one iteration is used for LMPC and LMPC with the Kalman filter.

Initially, the discharges in the main river part are all3 m3=s, whereas there is no flow in the reservoir and the two side channels. This indicates that the second and the third gates are closed. The most downstream water level of the second and the fourth channels is 5 m, whereas the most downstream water level of the reservoir is 4 m. The set points for all the water levels are equal to the initial water levels. The discharge at the end of the fourth channel should remain between 0 and5 m3=s. All the gate openings should be be-tween 0 and 5 m, and they can only be changed over a maximum distance of 20 cm every time step. The differences between the ini-tial water levels and the safety levels are 0.6, 0.6, and 1.2 m for Channels 1, 2, and 5, Channels 3, 4, and 7, and the reservoir, and the differences with respect to the flood levels are 1, 0.8, and 2 m, respectively. The disturbance signal is depicted in Fig.5. Initially, the discharge entering the river system is3 m3=s, which is

Fig. 5. Disturbance signal QoutðtÞ

increased after 4 h to 12 m3=s and is decreased back to 3 m3=s 6 h later.

Table2contains the diagonal elements for the weight matrices W, S, V, and R and the elements for the weight vectors s and v. As explained in the“Problem Formulation” section, during normal op-eration the water levels of the Channels 1 and 4 and of the reservoir are the most important ones. Therefore, the weights of W corre-sponding to these water levels are set to 10, 1, and 1,000, respec-tively. All the other weights for the other water levels are fixed to 0.001. The weights for the discharges for all the channels are set to 0.001. If the buffer capacity is used to prevent the water levels from violating their safety limits, the set points for the most downstream water levels (Channel 4) are set below the set point of the reservoir, and the corresponding weights in W are increased to 1,000 such that the buffer capacity can be recovered. All the diagonal elements of R are set to 50, the weights in S and s for the slack variable ξ associated to the reservoir are set to106, the weights associated to the other channels are set to105, the weights in V and v for the slack variable ζ are assumed equal to 108, rc is set to 1.2, and the size of the prediction window NPis 15.

Fig.6shows the evolution of the maximal difference between the water levels and their set points for the three controllers. To keep the number of plots as low as possible, the results for the up-stream part of the river (Channels 1, 2 and 5) are plotted together (the top plot) as are the results for the downstream part of the river (Channels 3, 4, and 7, the middle plot). The bottom plot shows the results for the reservoir (Channel 6). Fig.7shows the applied con-trol actions for the three gates and QoutðtÞ by each of the controllers. The first conclusion that can be made from Figs.6and7is that the differences in the results between the NMPC and the LMPC are small. The control actions and the water levels have a similar gen-eral trend for both controllers. Also, the addition of the estimator in the control loop has only a limited effect on the control perfor-mance. There is only a noticeable difference in the control actions for the three controllers during the period of heavy rainfall and right after this period. However, because the trend of the control actions for all the three controllers is similar, the water levels show the same behavior. Initially, all the three controllers try to avoid using the water reservoir and react to the future rainfall by decreasing the upstream water levels. This is done by increasing the gate opening of the first gate and increasing the discharge at the end of the river to its maximal value. Once the water levels risk violating the safety levels, the controllers start using the water reservoir by opening the second and the third gates. Because of the magnitude of the

upstream disturbance, the controllers cannot avoid the violation of the safety levels; however, they can prevent the violation of the flood levels or at worst minimize their violation. Table3shows the maximal difference between the water levels of each channel and their flood levels for the three controllers. A positive value corre-sponds to the maximal violation of the flood levels, whereas a neg-ative value corresponds to the minimal margin before a flood level

(a)

(b)

(c)

Fig. 6. Evolution of maximal difference between water levels and their set points for three controllers; plots show maximal deviation of water levels at every time step of Channels (a) 1, 2, and 5; (b) 3, 4, and 7; (c) 6; flood and safety levels correspond to horizontal lines

Table 2. Diagonal Elements for Weight Matrices W, S, V, and R and Elements for Weight Vectors s and v Elements Channel 1 2 3 4 5 6 7 W ∈ R333×333 Water levels 10 · 1₁₇ 0.001 · 1₂₁ 0.001 · 1₂₁ 1 · 1₂₁ 0.001 · 1₁₁ 1,000 · 1₆₁ 0.001 · 1₁₁ 1,000 · 121a Discharges 0.001 · 1₁₈ 0.001 · 1₂₂ 0.001 · 1₂₂ 0.001 · 1₂₂ 0.001 · 1₁₂ 0.001 · 1₆₂ 0.001 · 1₁₂ S ∈ R7×7 Safety levels 105 105 105 105 105 106 105 s ∈ R7×1 Safety levels 105 105 105 105 105 106 105 V ∈ R7×7 Flood levels 108 108 108 108 108 108 108 v ∈ R7×1 Flood levels 108 108 108 108 108 108 108

R ∈ R4×4 Qgate1 Qgate2 Qgate3 Qout

Control actions 50 50 50 50

a_{Value used to recover the buffer capacity of the reservoir. The weighting factor for each of the 21 water levels changes from 1 to 1,000.}

is violated. There is almost no difference between NMPC and LMPC. There is no flooding for any of the channels with LMPC. Three channels are flooded with NMPC. However, these floodings can be neglected: the flood levels are violated with only 3 mm or less, and each flooding lasts for less than 1 min. The remaining buffer capacity in the water reservoir (Channel 6) is larger with NMPC than with LMPC. Adding a state estimator has only a minor negative influence on these results: the maximal flooding is in-creased to 2 cm, and each flooding lasts for less than 5 min. Once the disturbance is decreased again, the controllers first steer the water levels below the safety limits. Then, they keep the upstream water levels close to the safety limits (the opening of the first gate is decreased, and the second gate is completely closed), whereas the downstream water levels are further decreased below their initial set point because of the change in the matrix W and the set point r. To decrease these water levels as fast as possible, the controllers keep the downstream discharge at its maximal value while the opening of the third gate is as large as the surrounding water levels. This allows the buffer capacity of the water reservoir to be recovered quickly. Once the buffer capacity is recovered, the matrix W and the set point r are set to their initial value. The controllers close the second and the third gates to prevent the water levels in the reser-voir to increase. The controllers immediately bring the upstream and the downstream water levels to their original set point by open-ing the first gate and changopen-ing QoutðtÞ. The most significant neg-ative influence of the state estimator can be found in the water level

of the water reservoir at the end of the simulation: some part of the buffer capacity is not recovered. This can be seen in Table4, in which the average difference between the water levels and their set point for each channel at the end of the simulation is shown. This difference for NMPC and LMPC is negligible for every chan-nel. This is also the case when the Kalman filter is used except for the reservoir where there is a small difference of 2 cm.

To verify the generality of the results, a second simulation has been performed for a higher flood with a maximal discharge of 24 m3_{=s in combination with a larger reservoir. All the other} parameters of the controller, the state estimator, and the river sys-tem remained unchanged. The simulation results for this second test case are very similar as the results discussed previously. This indicates that the proposed control scheme can handle different sizes of floodings and reservoirs.

Table5contains the mean, the minimal, and the maximal com-putation times needed for computing the control actions for each step (prediction, linearization, optimization, conversion) separately and the total time by the three controllers. Because NMPC performs the predictions by using a nonlinear model, the time needed to per-form this step is much larger than that required for LMPC and LMPC with the Kalman filter. The time needed to solve the opti-mization problem with the NMPC controller is much larger than that required by the other controllers because of the need to perform several iterations to achieve convergence. The relevance of the conversion step on the total computation time needed by any of the controllers is practically negligible.

(a) (b)

(c) (d)

Fig. 7. Control actions of three controllers; horizontal lines are upper limits on control variables: (a) Gate 1; (b) Gate 2; (c) Gate 3; (d) discharge Q_outðtÞ at end of river system

Table 3. Maximal Flooding of Each Channel for Three Controllers

Channel NMPC (m) LMPC (m) LMPC + Kalman (m) 1 0.0029 −0.0263 0.0160 2 0.0023 −0.0258 0.0090 3 −0.0283 −0.0360 −0.0292 4 −0.0265 −0.0277 −0.0288 5 0.0024 −0.0284 0.0051 6 (reservoir) −0.4054 −0.3670 −0.4037 7 −0.0290 −0.0293 −0.0249

Note: Negative value indicates that there is no flooding and some buffer capacity is left.

Table 4. Average Deviation in Meters of Each Channel for Three Controllers from Their Set Points at End of Simulation (t¼ 3,000 min)

Channel NMPC (m) LMPC (m) LMPC + Kalman (m) 1 −0.0012 −0.0012 −0.0012 2 −0.0012 −0.0012 −0.0012 3 −0.0081 −0.0084 −0.0078 4 −0.0081 −0.0083 −0.0078 5 −0.0012 −0.0012 −0.0012 6 (reservoir) 0.0015 0.0011 0.0200 7 −0.0081 −0.0083 −0.0078

Table6shows the quantification of the total amount of control actions for the three controllers. For each gate m, the total gate movement is calculated with the formula

γgatem ¼

X NT−1

k¼0 ju

mðk þ 1Þ − umðkÞj ð38Þ

where NT= total number of time instants during the simulation. For Qout, the total amount of control action is quantified with the in-tegral of the absolute value of the changes of Qoutwith respect to its initial value:

γQout ¼

Z _T

0 jQoutðtÞ − Qoutð0Þjdt ð39Þ where T= end time of the simulation. This indicates by how much the controller let Qoutdeviate from its normal operation point. The higher these numbers, the more demanding the controller is for the hydraulic structures. The third gate has been used the most by the three controllers. This is for releasing the excess water in the water reservoir. Overall, the NMPC changes the gate opening more frequently than the other two controllers, especially the first gate. The difference in total gate movement for LMPC without or with the Kalman filter is mainly explained by the third gate. This gate shows some oscillations when Kalman filtering is used, which are not present when the state estimator is absent. Because all the con-trollers use the Qoutin a very similar fashion, the difference forγQout

is negligible. One can always tune the controllers further to elimi-nate the oscillations or to reduce the total movement of the gates by decreasing the upper limit on the rate of change constraint (Δmax) or by increasing the value of the diagonal elements in the weight matrix R. This has the advantage that the lifespan of the equipment can be increased, but it can at the same time increase the risk of flooding. One needs to make a trade-off between the life cycle of the structures and the risk of flooding.

Conclusions

This paper has shown that model predictive control can be used for set-point and flood control of river systems with multiple channels,

gates, junctions, and a water reservoir. The optimization problem solved by MPC at every time step can be formulated such that the controller only uses the water reservoir when there is a risk of flooding. By modifying some of the parameters of the optimization problem, the same controller can be used to empty the water res-ervoir quickly when there is no risk of flooding. When the buffer capacity is restored, the controller automatically focuses on set-point control of the other water levels in the river system.

Given that the gate equations are the main sources of nonlinear-ities, the proposed controllers work with the gate discharges as control variables in the optimization problem instead of the gate positions. Then, a conversion step is performed to calculate the cor-responding gate positions. Because the remaining equations of the nonlinear model (the Saint-Venant equations) can accurately be ap-proximated with their linear counterpart (for the problem addressed in this paper) for designing the controllers, the LMPC achieves a similar control performance as the NMPC. The big advantage with LMPC is the reduction in computation time.

Adding a Kalman filter to the control loop to estimate the unmeasured states has only a minor influence on the control per-formance. Only a small increase in the offset of the water levels in the water reservoir was observed.

Future work will be focused on applying the techniques pre-sented in this paper to a bigger river system, the Demer in Belgium, which has more water reservoirs, channels, gates and junctions. Because of the big number of states needed to model the behavior of this system, model reduction techniques will studied, e.g., proper orthogonal decomposition and Galerkin projection, and distributed control strategies.

Acknowledgments

The authors would like to thank the reviewers of this manuscript for their suggestions and remarks that undoubtedly contributed to the improvement of our initial submission.

Research was supported by Research Council KUL: GOA/10/09 MaNet, PFV/10/002 (OPTEC), FWO: Ph.D./postdoctoral grants and projects: G.0320.08 (convex MPC), G.0558.08 (Robust MHE), Belgian Federal Science Policy Office: IUAP P7/ (DYSCO, Dynamical systems, control and optimization, 2012-2017).

References

Barjas Blanco, T., Willems, P., Chiang, P. P.-K., Haverbeke, N., Berlamont, J., and De Moor, B. (2010).“Flood regulation using nonlinear model predictive control.”Control Eng. Pract., 18(10), 1147–1157. Breckpot, M., Agudelo, O. M., De Moor, B. (2012a).“Control of a single

reach with model predictive control.” River Flow 2012, Proc., Int. Conf. on Fluvial Hydraulics, R. E. Murillo Muñoz, ed., CRC Press, Boca Raton, FL, 1021–1028.

Table 5. Average, Minimal, and Maximal Computation Time (in Seconds) Needed for Each Step Separately and All Together by Three Controllers at Every Time Step

Steps

NMPC LMPC LMPC + Kalman

Mean Min. Max. Mean Min. Max. Mean Min. Max.

Prediction (s) 17.60 7.06 21.81 0.09 0.05 0.12 0.09 0.06 0.14

Linearization (s) 0.06 0.06 0.08 — — — — — —

Optimization (s) 6.26 4.53 10.17 4.98 3.50 8.41 5.07 3.31 7.76

Conversion (s) <0.01 <0.01 0.01 <0.01 <0.01 0.01 <0.01 <0.01 0.01

Total (s) 23.93 12.23 30.59 5.07 3.60 8.50 5.16 3.54 7.69

Note: Simulation was performed on a personal computer with a 2.8 GHz Intel Core i7 CPU and a random-access memory of 8 GB; Min.= minumum; Max.= maximum.

Table 6. Quantification of Total Amount of Control Actions for Three Controllers

Control variable NMPC LMPC LMPC + Kalman

Gate 1 (m) 8.04 4.39 5.47

Gate 2 (m) 3.39 2.29 1.72

Gate 3 (m) 14.16 12.39 15.08

Total gate movement (m) 25.59 19.06 22.26

γQout ðm

3_{Þ 257,936 258,504 255,667}

Breckpot, M., Agudelo, O. M., and De Moor, B. (2012b).“Model predic-tive control of a river system with two reaches.” Proc., 51st IEEE Conf. on Decision and Control, IEEE, New York.

Breckpot, M., Barjas Blanco, T., and De Moor, B. (2010).“Flood control of rivers with nonlinear model predictive control and moving horizon estimation.” Proc., 49th IEEE Conf. on Decision and Control, IEEE, New York, 6107–6112.

Burt, C. M., Mills, R. S., Khalsa, R. D., and Ruiz, V. (1998).“Improved proportional-integral (PI) logic for canal automation.”J. Irrig. Drain. Eng., 124(1), 53–57.

Chaudry, M. H. (2008). Open-channel flow, 2nd Ed., Springer, New York. Chow, V. T. (1959). Open-channel hydraulics, McGraw-Hill, New York. Clemmens, A. J., Bautista, E., Wahlin, B. T., and Strand, R. J. (2005). “Simulation of automatic canal control systems.”J. Irrig. Drain. Eng., 131(4), 324–335.

Cunge, J. A., Holly, F. M., and Verwey, A. (1980). Practical aspects of computational river hydraulics, Pitman, London.

Franklin, G. F., Powell, D. J., and Workman, M. L. (1997). Digital control of dynamic systems, 3rd Ed., Addison-Wesley, Boston, MA. Henderson, F. M. (1966). Open channel flow, Macmillan, New York. Hovd, M., and Braatz, R. D. (2001).“Handling state and output constraints

in MPC using time-dependent weights.” Proc., American Control Conf. (ACC), IEEE, New York.

Kalman, R. E. (1960).“A new approach to linear filtering and prediction problems.”J. Basic Eng., 82, 35–45.

Kwakernaak, H., and Sivan, R. (1972). Linear optimal control systems, Wiley-Interscience, New York.

Liggett, J. A., and Cunge, J. A. (1975). Numerical methods of solution of the unsteady flow equations, unsteady flow in open channels, Chapter 4, Water Resources Publications, Fort Collins, CO.

Lin, C. H., Yen, J. F., and Tsai, C. T. (2002).“Influence of sluice gate con-traction coefficient on distinguishing condition.”J. Irrig. Drain. Eng., 128(4), 249–252.

Litrico, X., and Fromion, V. (2009). Modeling and control of hydrosystems, Springer-Verlag, New York.

Litrico, X., Fromion, V., and Baume, J. P. (2006).“Tuning of robust distant downstream PI controllers for an irrigation canal pool—II. Implemen-tation issues.”J. Irrig. Drain. Eng., 132(4), 369–379.

Malaterre, P. O., and Baume, J. P. (1999).“Optimum choice of control action variables and linked algorithms. Comparison of different alternatives.” Proc., Workshop on Modernization of Irrigation Water Delivery Systems, USCID, Phoenix, 387–406.

Malaterre, P. O., Rogers, D. C., and Schuurmans, J. (1998).“Classification of canal control algorithms.”J. Irrig. Drain. Eng., 124(1), 3–10. MATLAB 3 (2011a). [Computer software]. MathWorks, Natick, MA. Mayne, Q., Rawlings, J. B., Rao, C. V., and Scokaert, P. O. M. (2000).

“Constrained model predictive control: Stability and optimality.” Automatica, 36(6), 789–814.

Negenborn, R. R., van Overloop, P.-J., Keviczky, T., and De Schutter, B. (2009). “Distributed model predictive control of irrigation canals.” Networks Heterog. Media, 4(2), 359–380.

Nocedal, J., and Wright, S. J. (2000). Numerical optimization, Springer, New York.

Puig, V., et al. (2009). “Optimal predictive control of water transport systems: Arrêt-Darré/Arros case study.” Water Sci. Technol., 60(8), 2125–2133.

Qin, S., and Badgwell, T. (2003).“A survey of industrial model predictive control technology.”Control Eng. Pract., 11(7), 733–764.

Romera, J., Ocampo-Martnez, C., Puig, V., Quevedo, J., Garca, P., and Prez, G. (2011).“Flooding management using hybrid model predictive control at the Ebro River in Spain.” Proc., 8th IWA Symp. on System Analysis and Integrated Assessment, IWA, San Sebastian, Spain, 1–8. Rossiter, J. (2003). Model-based predictive control, CRC, Boca Raton, FL. Schuurmans, J. (1997).“Control of water levels in open channels.” Ph.D.

thesis, Delft Univ. of Technology, Delft, Netherlands.

Scokaert, P. O. M., and Rawlings, J. B. (1999).“Feasibility issues in linear model predictive control.”AIChE J., 45(8), 1649–1659.

Sepúlveda, C., Gómez, M., and Rodellar, J. (2009).“Benchmark of dis-charge calibration methods for submerged sluice gates.”J. Irrig. Drain. Eng., 135(5), 676–682.

Strelkoff, T. S., and Falvey, H. T. (1993).“Numerical methods used to model unsteady canal flow.”J. Irrig. Drain. Eng., 119(4), 637–655. Sturm, T. W. (2001). Open channel hydraulics. McGraw-Hill, New York. Thai, T. (2005).“Numerical methods for parameter estimation and optimal control for the Red River Network.” Ph.D. thesis, Universität Heidelberg, Heidelberg, Germany.

Van Overloop, P. J. (2006).“Model predictive control of open water sys-tems.” Ph.D. thesis, Technische Universiteit Delft, Delft, Netherlands. Van Overloop, P. J., Clemmens, A. J., Strand, R. J., Wagemaker, R. M. J., and Bautista, E. (2010).“Real-time implementation of model predictive control on Maricopa-Stanfield irrigation and drainage district’s WM canal.”J. Irrig. Drain. Eng., 136(11), 747–756.

Van Overloop, P.-J., Schuurmans, J., Brouwer, R., and Burt, C. M. (2005). “Multiple-model optimization of proportional integral controllers on canals.”J. Irrig. Drain. Eng., 131(2), 190–196.

Van Overloop, P.-J., Weijs, S., and Dijkstra, S. (2008).“Multiple model predictive control on a drainage canal system.”Control Eng. Pract., 16(5), 531–540.

Wahlin, B., and Clemmens, A. J. (2006).“Automatic downstream water-level feedback control of branching canal networks: Theory.”J. Irrig. Drain. Eng., 132(3), 208–219.

Xu, M., van Overloop, P. J., and van de Giesen, N. C. (2011).“On the study of control effectiveness and computational efficiency of reduced Saint-Venant model in model predictive control of open channel flow.”Adv. Water Resour., 34(2), 282–290.