P Segovia
email: pablo.segovia@upc.edu

L Rajaoarisoa

F Nejjari

E Duviella

V Puig

Model predictive control and moving horizon estimation for water level 1 regulation in inland waterways 2

Keywords: Inland waterways, control-oriented modeling, model predictive control, moving 9 horizon estimation, time-delay systems

This work regards the design of optimization techniques for the purposes of state estimation and control in the framework of inland waterways, often characterized by negligible bottom slopes and large time delays. The derived control-oriented model allows these issues to be handled in a suitable manner. Then, the analogous moving horizon estimation and model predictive control techniques are applied in a centralized manner to estimate the unmeasurable states and fulfill the operational goals, respectively. Finally, the performance of the methodology is tested in simulation by means of a realistic case study based on part of the inland waterways in the north of France. The results show that the proposed methodology is able to guarantee the navigability condition, as well as the other operational goals.

Introduction 11

Inland waterways are large-scale systems, composed of natural rivers and artificial canals, 12 used mainly for transportation of passengers and freight. Indeed, fluvial transport constitutes 13 an environmentally-friendly alternative to the traditional rail and road transport modes [START_REF] Mihic | Policy and promotion of sustainable inland waterway transport in Europe-Danube River[END_REF]. In 14 order to facilitate the study of such systems, they are usually decomposed into reaches, which are 15 parts of a stream of water between two hydraulic structures such as gates or locks. Reaches are 16 usually characterized by negligible bottom slopes, and therefore the backwater effect becomes of increasing importance. This effect takes place at the downstream hydraulic structure of a canal: when the water waves impact upon the structure, the water can flow back to the upstream end, resulting in a back-and-forth mass transport known as the resonance phenomena. Other applications in the framework of free-surface water systems (e.g., sewage systems, irrigation and drainage canals) are generally not affected by the backwater effect. On the other hand, the dynamics of inland waterways are rather slow, which translates into large time delays in the network. These two features, typical of inland navigation networks, complicate its management.

Inland waterways management aims at allocating the available water resources to meet the desired objectives. The most important goal nonetheless consists in guaranteeing the navigability R3.6 condition, i.e., ensuring that the water levels are such that vessels can travel safely. Indeed, a setpoint is defined for each reach, known as the Normal Navigation Level (NNL). The objective is to keep the water levels as close as possible to the setpoints. In addition, upper and lower level bounds known as the Lower and Higher Navigation Levels (LNL and HNL, respectively) are determined around the setpoint, thus defining the navigation rectangle. If the water level of a reach is outside of the rectangle, the navigability condition can no longer be guaranteed. Other common objectives regard, for instance, minimizing the operational cost and ensuring a long lifespan of the R1.1 actuators, e.g., gates, weirs, pumps and valves.

Inland waterways are large-scale, complex systems affected by some phenomena that are not R3.17 easy to account for, e.g., demands, uncontrolled inputs, rainfall and seepage, and thus it is not easy to assume that the water levels remain within the bounds in the absence of a control strategy. This fact motivates the design of monitoring and control approaches for inland waterways. Therefore, a suitable control policy is needed in order to fulfill the objectives. Optimal control techniques have been investigated for a long time to this end. In particular, Model Predictive Control (MPC) has received a lot of attention due to its adequacy to deal with these kinds of problems. Its underlying principle consists in using a dynamic model of the process to predict the effect of manipulable inputs, subject to operational restrictions, so that the performance of the plant is optimal regarding the chosen criteria. Its ease of understanding and application has fostered its use in many different domains. In the framework of water systems, MPC has also been widely used. For instance, in [START_REF] Van Overloop | Predictive control for national water flow optimization in The Netherlands[END_REF] it was applied at regional and national scales to protect against high river flows and sea tides, to ensure navigation and to supply water during dry periods. This technique was also employed in [START_REF] Fele | Coalitional model predictive control of an irrigation canal[END_REF], using a coalitional approach to find the best compromise between communication costs and control performance for irrigation canals. The combined water supply and navigability of river systems was tackled using MPC in [START_REF] Puig | Model predictive control for combined water supply and navigability/sustainability in river systems[END_REF][START_REF] Horváth | Enhancing inland navigation by model predictive control of water levels: the Cuinchy-Fontinettes case[END_REF]. A nonlinear economic MPC was designed in [START_REF] Wang | Non-linear economic model predictive control of water distribution networks[END_REF] for water distribution networks, aiming at minimizing the economic costs associated to water treatment and pumping. A comparison of non-centralized MPC strategies for irrigation canals was carried out in [START_REF] Nguyen | Distributed Model Predictive Control of Irrigation Systems using Cooperative Controllers[END_REF], validating the benefits of cooperative control.

The use of MPC requires the vector of states to be known at current time. This information is used as the starting point to compute the set of control actions that must be applied during the future horizon. In general, the measurements of all states are not available, and thus estimates of unmeasured states must be provided to the MPC using observers. Although there exist many possibilities to estimate the states, this work employs the Moving Horizon Estimation (MHE), which is considered as the dual problem of MPC. Their combination is especially attractive since the MHE formulation corresponds also to an online optimization problem that can explicitly handle constraints [START_REF] Copp | Simultaneous nonlinear model predictive control and state estimation[END_REF]. Unlike MPC, this estimation technique started receiving wider attention only in the recent years [START_REF] Tenny | Efficient moving horizon estimation and nonlinear model predictive control[END_REF]. Indeed, the combination of MPC and MHE has been applied in diverse fields such as autonomous agricultural vehicles [START_REF] Kraus | Moving horizon estimation and nonlinear model predictive control for autonomous agricultural vehicles[END_REF], unmanned aerial vehicles [START_REF] Quintero | Robust UAV coordination for target tracking using outputfeedback model predictive control with moving horizon estimation[END_REF], preventive sensor maintenance [START_REF] Lao | Real-time preventive sensor maintenance using robust moving horizon estimation and economic model predictive control[END_REF], airborne wind energy systems [START_REF] Vukov | Real-time nonlinear MPC and MHE for a large-scale mechatronic application[END_REF] and blood glucose regulation [START_REF] Copp | Simultaneous model predictive control and moving horizon estimation for blood glucose regulation in Type 1 diabetes[END_REF]. Concerning water systems, the combination of these techniques is not so common, although it has been used for flood prevention in rivers [START_REF] Breckpot | Flood control of rivers with nonlinear model predictive control and moving horizon estimation[END_REF] and pollution mitigation for combined sewer networks [START_REF] Joseph-Duran | Output-feedback control of combined sewer networks through receding horizon control with moving horizon estimation[END_REF].

Summary of the paper and contribution

This work regards the development of MPC and MHE for inland waterways, aiming not only at guaranteeing the navigability condition of the network, but also at ensuring the rest of operational goals, linked to optimizing the control effort and reducing wear and tear of the equipment. The first steps, although focusing only on the controller design and leaving aside the topic of state estimation, were carried out in [START_REF] Segovia | Input-delay model predictive control of inland waterways considering the backwater effect[END_REF]. Since a model of the process is at the core of MPC and MHE, a control-oriented methodology based on the classical Integrator Delay Zero (IDZ) model [START_REF] Litrico | Modeling and Control of Hydrosystems[END_REF] is also developed. This new modeling formulation takes into account the challenging features of these kinds of systems, and provides flexibility in terms of adding reaches to or removing them from the case study.

The contributions of this work to the current state of the art are listed below:

R3.1

• Inland waterways are large-scale systems with slow dynamics, which result in large time delays. A common approach to model a system with delays consists in augmenting the system with more states, with as many new states as delayed samples the system exhibits [START_REF] Santos | Explicit input-delay compensation for robustness improvement in MPC[END_REF]. Thus, the time delays have a direct influence on how large the resulting augmented system is. By contrast, this paper proposes a delayed representation by means of additional matrices where this augmentation procedure is not necessary.

• In order for distribution networks such as inland waterways to be completely described, it is not enough to represent their dynamic behavior. Indeed, these systems are only fully described when their static behavior is characterized. In the case of inland waterways, the mass balances at the junctions must be satisfied. This leads to a model representation that falls under a particular family of systems known as descriptor or differential-algebraic systems. The classical control theory was originally formulated for dynamic systems (without the static part), and these results are usually more involved for descriptor systems.

• Therefore, this paper proposes a new model formulation for inland waterways, for which the use of first principles yields a delayed descriptor formulation. These two features cause that the resulting model cannot be represented using the standard state-space formulation. To the best of the authors' knowledge, the formulation derived in this work constitutes a novel result.

• In addition, although this work employs only one linear model (this fact is duly justified in Section 3.4), the extension to a linear-parameter varying (LPV) formulation is rather direct using the results in [START_REF] Bolea | Gain-scheduling multivariable LPV control of an irrigation canal system[END_REF], which allows to retain the nonlinear behavior of the system in the necessary cases.

• Standard tools for control and state estimation such as the Linear Quadratic Regulator (LQR) and the Kalman filter need to be extended to deal with delayed descriptor systems. Moreover, they cannot deal with input and state/output constraints. This work proposes instead the combined use of MPC and MHE, which can be easily adapted for this model formulation. Furthermore, constraints on the inputs and the states/outputs are dealt with by these techniques in a natural manner. The review of the literature shows that the combination of MPC and MHE for delayed descriptor systems, especially in the field of open-flow water systems, is a novel approach.

The rest of the paper is organized as follows: Section 2 introduces the general waterway management problem. The control-oriented model is derived in Section 3. The MPC and MHE optimization problems are designed in Section 4. The proposed methodology is tested in Section 5 by considering a realistic case study, based on the inland navigation network in the north of France.

Finally, Section 6 draws conclusions and outlines future steps.

Notation

Throughout this paper, let R γ denote the set of column real vectors of length γ, and let R ≥0 denote the set of real non-negative scalars. Scalars are denoted by either lowercase or uppercase letters (e.g., α, a and A); vectors, by bold lowercase letters (e.g., a and b); and matrices, by bold uppercase letters (e.g., A and B). Furthermore, all vectors are column vectors unless otherwise stated, 0 denotes a zero column vector of suitable dimensions and I δ denotes the identity matrix of dimension δ. Transposition is denoted with the superscript , and the operators <, ≤, =, ≥ and > denote element-wise relations of vectors.

The waterways management problem

As mentioned before, the management of inland waterways aims at ensuring that the trans-R3.3 portation of passengers and freight is carried out safely. To guarantee seamless transport chains, the water levels must be kept inside the navigation rectangle. Furthermore, it is important that the water resources are dispatched in an optimal manner, i.e., minimizing their losses.

To this end, cross structures are operated in the waterways to regulate the levels of the reaches.

In particular, this work considers two kinds of structures: undershot gates and sharp crested weirs.

An undershot gate is a bottom opening in a wall, whose height can be regulated. Conversely, in the case of a weir, the water flows over its crest, whose elevation is also adjustable. A schematic representation of an undershot gate and a weir is given in Fig. 1, and their equations are given in Appendix B1. Note that q is the flow through the cross structure, u is the opening/elevation, and y 1 and y 2 are the upstream and downstream water levels, respectively. A two-level control architecture is typically considered in this environment [START_REF] Puig | Model predictive control for combined water supply and navigability/sustainability in river systems[END_REF]: the global control level (represented by the MPC) determines the setpoints and sends them to the local controllers available at each control structure. In turn, these local controllers are in charge of ensuring that the adequate flows are supplied through the actuators. However, this work deals only with the global control level, thus assuming that the local controllers are able to perform as desired.

The control strategy must also reject the disturbances that affect the system and interfere with the control objectives. In the framework of inland waterways, these disturbances refer to the request of lock operations by the boat masters. Indeed, boats navigate along the network until their final destination, probably along several reaches in their way. The access from one reach to the adjacent one is granted by means of locks, which typically consist in enclosures that enable boats to overcome the difference in elevation between the reaches. Lock operations require large water volumes to be withdrawn from the origin reach, which are then discharged into the destination reach. The reason for considering lock operations as disturbances is that they cannot be postponed for a long time from the moment a boat reaches a lock, and it is therefore not possible to schedule them in an optimal manner.

Control-oriented modeling

An inland waterways model can be regarded as composed of a set of elements, which are introduced and described below. Note that the physical nature of the variables, e.g., water levels, flows, openings and elevations, as well as other elements in the waterways, constrain the performance of the system.

Actuators

Gates and weirs are used to regulate the water levels in the reaches. Discharges and openings/elevations can be used in the automated control of canals, as shown in Fig. 1. The global control level must compute the optimal action and send it to the slave controller that operates the gate or weir. If the discharge is used as the control variable, the slave controller must convert the given discharge into an equivalent opening or elevation, which is not as straightforward as inverting the discharge equation [START_REF] Litrico | Conversion from discharge to gate opening for the control of irrigation canals[END_REF]. Furthermore, choosing the openings and elevations allows to link them with the local discharges and the upstream and downstream water levels at the structure, thus taking into account such complex dynamics [START_REF] Malaterre | Regulation of irrigation canals[END_REF]. For these reasons, the openings and elevations are chosen as control variables in this work. The conversion is carried out following the methodology described in [START_REF] Bolea | Linear parameter varying modeling and identification for real-time control of open-flow irrigation canals[END_REF], which basically consists in using linearized equations that describe the relationship between openings and discharges.

These elements have lower and upper operating limits

u m ≤ u m k ≤ u m , m = 1, ..., N m , (1)
where u m and u m are the lower and upper opening/elevation limits of the m-th actuator, and N m is the total number of actuators in the system.

The type of flow at the structure determines the general linearized equation to be used:

• The free-flow case is characterized by critical or super-critical flow at the structure, which overrides the effect of the downstream water level on the gate discharge. The linearized expression reads as

q (1) 2 (s) ≈ q (2) 1 (s) = k (1) y y (1) 2 (s) + k u u(s). (2)
• The submerged flow case is characterized by sub-critical flow at the structure. In this case, the discharge is affected by the downstream water level:

q (1) 2 (s) ≈ q (2) 1 (s) = k (1) y y (1)
2 (s) + k (2) y y

(2)

1 (s) + k u u(s). (3)
In both cases, q

2 is the inflow of the structure at the downstream end of reach 1, q

1 is the outflow of the structure at the upstream end of reach 2, y

2 is the water level upstream of the structure, y

1 is the water level downstream of the structure, u is the opening/elevation and k

(1)

y , k (2)
y and k u are the coefficients obtained in the linearization of the nonlinear equations of the gates R3.4

and weirs, given by (32).

Disturbances

Systems and processes are usually affected by disturbances, denoted in this work by d k . Therefore, the control strategy must minimize their effect on the system. As mentioned in Section 2, these disturbances correspond to lock operations, which makes it more difficult to stay close to the setpoints. Although lock operations are rather unpredictable and cannot be postponed for a long time, they can be somewhat anticipated. Indeed, when a boat passes through a lock, its manager informs the rest of the managers. In this way, the arrival time of the boat to the adjacent locks can be predicted, taking into account the distance and the average speed of the boat, which yields a close approximation, with an error of only several minutes. This allows the lock managers to elaborate lock operation time-series profiles ahead of time.

Nodes

Inland waterways are characterized by distributaries, i.e., streams that branch off from the main stream and flow away. When water streams flow into larger streams or lakes, they are referred to as tributaries. The locations in which these splittings and mergings take place are called nodes.

They are regarded as mass balance relations modeled as equality constraints given by:

0 = E u u k + E un u k-n + E d d k + E dn d k-n . (4)
Matrices E u and E un have as many rows as nodes are in the studied system, and as many R2.5

columns as controlled inputs are available. Therefore, each equation in ( 4) establishes a link among the variables involved (mass balance at the node), and thus reduces one degree of freedom.

Note that both the controlled inputs and the disturbances have an immediate and a delayed effect on the system. The delayed effect must be taken into account at the controller and estimator design stages. This issue is conveniently addressed in Section 4.

Reaches

An accurate mathematical representation of the dynamics of inland waterways is required in order to apply the MPC and MHE techniques. Indeed, a model of the system is needed in the control design stage to compute the predicted output at future time instants. Likewise, it is used to align measured and predicted values of the process, which results in the optimal state estimates.

Therefore, a special effort has to be put in the computation of a precise model.

The Saint-Venant nonlinear partial differential equations allow the accurate representation of the dynamics of open-flow water systems [START_REF] Chow | Open-channel hydraulics[END_REF]. However, the lack of a general analytical solution, as well as their extreme sensitivity to geometry errors and unmodeled dynamics, render them inadequate for control purposes. The nonlinear behavior can be retained, for instance, by considering R1.2

LPV models, which describe a class of nonlinear systems that can be modeled as parametrized linear systems, each of them designed at a different operating point (average flow along the canal). A simpler solution is to use only one linear model (obtained by linearizing the original Saint-Venant equations around an operating point and considering simplifying assumptions), provided that a single operating point is enough to describe the system dynamics. Examples of linearized models are the Integrator Delay (ID) model [START_REF] Schuurmans | Modeling of irrigation and drainage canals for controller design[END_REF], the Integrator Delay Zero (IDZ) model [START_REF] Litrico | Modeling and Control of Hydrosystems[END_REF], the Integrator Resonance (IR) model [START_REF] Van Overloop | Identification of resonance waves in open water channels[END_REF], and grey-box [START_REF] Horváth | Gray-box model of inland navigation channel: application to the Cuinchy-Fontinettes reach[END_REF] and black-box [START_REF] Weyer | System identification of an open water channel[END_REF] models. Given that the average R1.4-R1.7

flows in this work do not deviate too much from the operating point, it is not necessary to resort to nonlinear modeling strategies such as LPV models. Instead, it is enough to consider one linear model. Although the IR model might appear as the most suitable option since it explicitly takes into account the resonance phenomena, the IDZ model is chosen over the IR model for a number of reasons. First, the reflecting waves in the resonance phenomena are especially predominant in short and deep channels. Instead, the navigation reaches considered in this work are rather long.

Moreover, the IDZ model has proved to perform well when used for control purposes for these kinds of systems [START_REF] Segovia | Decentralized control of inland navigation networks with distributaries: application to navigation canals in the north of France[END_REF].

The IDZ input-output model links the discharges and the water levels at the boundaries of a reach and is given by:

  y 1 (s) y 2 (s)   =   p 11 (s) p 12 (s) p 21 (s) p 22 (s)   P   q 1 (s) q 2 (s)   , (5)
where the subscripts 1 and 2 indicate the initial (upstream) and final (downstream) ends of the R2.3

reach, y 1 (s) and y 2 (s) are the upstream and downstream water levels, q 1 (s) and q 2 (s) are the upstream inflow and downstream outflow, and p ij (s) are the different IDZ terms

p ij (s) = z ij • s + 1 A ij • s e -τ ij •s . (6)
Remark: In this work, q 1 (s) ∈ R ≥0 and q 2 (s) ∈ R ≥0 . However, the latter is an outflow, which means that it causes the water levels to diminish. Therefore, p 12 (s) and p 22 (s) are negative.

R2.1

As it can be seen in ( 6), the IDZ model contains an integrator whose gain is given by 1/A ij , a R2.1

time delay τ ij and a zero given by -1/z ij , for i, j = {1, 2}. The system can be characterized by two different behaviors in the frequency domain. In low frequencies, the behavior of the system is similar to a tank that is being filled and/or emptied. In this situation, the integrator gain and the time delay have a predominant role. The former reflects how the volume changes according to the water level variation, whereas the latter expresses the minimum time that a perturbation requires to travel from one end of the canal to the other one. Two different time delays are defined:

τ 12 = L C w -V , τ 21 = L C w + V . (7)
Equation ( 7) corresponds to the case in which both the wave celerity C w and the wave velocity V are constant. More precisely, the celerity is defined as the relative velocity of a wave with R1.3 respect to the fluid in which it travels, whereas the velocity measures the variation of the particles' position of a fluid with respect to time. In particular, τ 21 is measured from the upstream end to the downstream end, while τ 12 is measured in the inverse direction. Note also that τ 11 = τ 22 = 0, since the discharges are assumed to have an immediate effect at the locations where they take place.

On the other hand, the high frequency phenomena is approximated by the zero of p ij (s). More specifically, its constant gain approximates the oscillating modes caused by the gravity waves, which are predominant in the high frequencies.

The parameters of the first equation in ( 5) are linked to the upstream water level, while those in the second equation are linked to the downstream water level. The notation of the parameters is modified based on this fact, and is as follows:

A 11 = A 12 = A u , A 21 = A 22 = A d , τ 12 = τ u and τ 21 = τ d .
Note also that the complete model is taken into consideration. Indeed, it is common practice to design only downstream water level controllers [START_REF] Litrico | Modeling and Control of Hydrosystems[END_REF]. Instead, the full model allows to take into account the backwater effect in the upstream water level, which is of relevance due to the negligible bottom slope of the reaches.

The navigability condition restricts the water levels in the reaches. This constraint might be relaxed for a short period of time, depending on factors such as the weather condition. Thus, a relaxation parameter α α α k is considered in the constraint, and a quadratic penalty on this parameter is included in the objective function.

The navigability condition is formulated as

y r -α k ≤ y k ≤ y r + α k , (8)
with y r and y r the lower (LNL) and upper (HNL) bounds of the NNL values, respectively. These relaxation parameters α k must satisfy

α k ≥ 0. (9)
Figure 3 depicts a waterway composed of several reaches for a better understanding of the R2.2

variables introduced to formulate the problem, and how they are linked to one another. Note that the locks are not depicted in this figure, but their operations are indicated using the variable d, as defined in Section 3.2.

Remark: The subscripts u and d mean at the upstream end and at the downstream end, respectively.

gate i-1 gate i u u q (i-1) q (i) q (i+1) (i-1) (i+1) u (i) y u (i-1) y d (i-1) y u (i) d y (i) y u (i+1) (i-1) (i+1) d d d (i) y d (i-2) reach i-1 reach i reach i+1 weir i+1

Final control-oriented model

The final, equivalent state-space representation is derived step by step. This new formulation R3.5 allows to coordinate current and delayed information in a systematic manner. Indeed, the model is described by variables with an immediate and a delayed effect. Those variables with a delayed effect are provided to the control and state estimation algorithms as parameters, ensuring that their values are taken into account adequately, which is crucial for a satisfactory performance of the algorithms. Remark: In the state-space model formulation, the notation q represents the discharges, whereas the variable u is saved for the openings and elevations, and will be used later on. y 1 (s) = p 11 (s) q 1 (s) -p 12 (s) q 2 (s), y 2 (s) = p 21 (s) q 1 (s) -p 22 (s) q 2 (s). [START_REF] Kraus | Moving horizon estimation and nonlinear model predictive control for autonomous agricultural vehicles[END_REF] Then, ( 6) is substituted in [START_REF] Kraus | Moving horizon estimation and nonlinear model predictive control for autonomous agricultural vehicles[END_REF], taking into account the parameter naming adopted in Section 3.4, leading to

y 1 (s) = z 11 s + 1 A u s q 1 (s) - z 12 s + 1 A u s e -τus q 2 (s), y 2 (s) = z 21 s + 1 A d s e -τ d s q 1 (s) - z 22 s + 1 A d s q 2 (s). (11)
In order to simplify the task, the delays are initially dropped, and then reincorporated when the state-space representation is obtained. Additionally, a convenient manipulation of (11) leads to

y 1 (s) = 1/A u s + z 11 A u q 1 (s) - 1/A u s + z 12 A u q 2 (s), y 2 (s) = 1/A d s + z 21 A d q 1 (s) - 1/A d s + z 22 A d q 2 (s). (12)
A standard transformation of (12) (see Chapter 2.5 in [START_REF] Ogata | Modern Control Engineering[END_REF] for more details) yields the statespace representation

ẋ(t) =   0 0 0 0   x(t) +   1 -1 1 -1   q(t), y(t) =   1 Au 0 0 1 A d   x(t) +   z 11 Au -z 12 Au z 21 A d -z 22 A d   q(t). (13)
Model ( 13) is discretized with a sampling time T s as follows:

x k+1 =   1 0 0 1   x k +   T s -T s T s -T s   q k , y k =   1 Au 0 0 1 A d   x k +   z 11 Au -z 12 Au z 21 A d -z 22 A d   q k . (14)
The time delays are re-incorporated into [START_REF] Copp | Simultaneous model predictive control and moving horizon estimation for blood glucose regulation in Type 1 diabetes[END_REF], which yields

x k+1 =   1 0 0 1   x k +   T s 0 0 -T s   q k +   0 -T s T s 0   q k-n , y k =   1 Au 0 0 1 A d   x k +   z 11 Au 0 0 -z 22 A d   q k +   0 -z 12 Au z 21 A d 0   q k-n , (15)
with q k-n the vector of discharges delayed n samples (n = τ /T s , with • the ceiling function).

In practice, τ d ≈ τ u , which leads to a single value of n.

Finally, the disturbances introduced in Section 3.2 are incorporated to the model. Since these R2.4

lock operations are also flows, and the locks are next to the actuators, their effect on the system is the same as the controlled discharges. Thus, the matrices for controlled discharges and disturbances are the same, leading to

x k+1 =   1 0 0 1   x k +   T s 0 0 -T s   q k +   0 -T s T s 0   q k-n +   T s 0 0 -T s   d k +   0 -T s T s 0   d k-n , y k =   1 Au 0 0 1 A d   x k +   z 11 Au 0 0 -z 22 A d   q k +   0 -z 12 Au z 21 A d 0   q k-n +   z 11 Au 0 0 -z 22 A d   d k +   0 -z 12 Au z 21 A d 0   d k-n . (16)
As mentioned before, ( 16) must be obtained for each reach in the case study. Then, q k and q k-n must be substituted in each case by either (2) or (3) accordingly. It can be anticipated that this substitution will cause delayed states to appear in the model. Indeed, q k = f y k , u k , and Although the step-by-step derivation of the final model is given in Section 5, it is convenient to present its final structure at this stage, since it will be used in Section 4 to design the controller and the estimator. Therefore, the general model formulation of a system with n x states, n u inputs and n y outputs is

thus q k-n = f y k-n , u k-n , with y k-n = g (x k-n ),
x k+1 = Ax k + A n x k-n + B u u k + B un u k-n + B d d k + B dn d k-n , (17a)
y k = Cx k + C n x k-n + D u u k + D un u k-n + D d d k + D dn d k-n , (17b)
with is the output equation. The mass balances given by ( 4) can be formulated by means of constraints, as it is now, or they can be incorporated into (17a) as shown in [START_REF] Grosso | Chance-constrained model predictive control for drinking water networks[END_REF].

x k ∈ R nx , u k ∈ R nu , y k ∈ R ny ,
Model [START_REF] Segovia | Input-delay model predictive control of inland waterways considering the backwater effect[END_REF] corresponds to the case of only one delay in the network. The general case for a system with multiple delays given by the set S = {n 1 , n 2 , ..., n p } reads as

x k+1 = Ax k + B u u k + B d d k + n i ∈S (A n i x k-n i + B un i u k-n i + B dn i d k-n i ) , (18a)
y k = Cx k + D u u k + D d d k + n i ∈S (C n i x k-n i + D un i u k-n i + D dn i d k-n i ) . (18b)
Furthermore, (4) must be modified as

0 = E u u k + E d d k + n i ∈S (E un i u k-n i + E dn i d k-n i ) ( 19
)
in order to account for multiple time delays.

Finally, note that the delayed terms in ( 17) and ( 18) cause these models to not be representable using the standard state-space formulation. While the theory of other classical state feedback control techniques might not be used for this representation, the combination of MPC and MHE can deal with these models in a suitable manner. Furthermore, this formulation also allows a flexible and more compact notation of a system with delayed variables. Indeed, a common approach to represent such systems consists in the augmentation procedure described in [START_REF] Santos | Explicit input-delay compensation for robustness improvement in MPC[END_REF], where the delay effect is incorporated as a dead-beat dynamic to obtain an undelayed representation. However, a downside of this methodology lies in the large dimensionality of the resulting description. By contrast, no augmented model needs to be derived in the case of MPC and MHE.

Control design and state estimation

This section deals with the control design and state estimation, for which the model of the system is needed. An approach is presented for each of them, namely MPC for control and MHE for state estimation, for which the complete structure of the resulting multi-objective optimization problems is given. The last part of the section provides some insight on the combination of MPC and MHE in simulation.

Control design: the MPC approach

Modern inland waterways are complex, multivariable systems whose management requires the use of advanced control methods [START_REF] Van Overloop | Model predictive control on open water systems[END_REF]. MPC is characterized by several interesting features that are very suitable for these kinds of systems [START_REF] Camacho | Model Predictive Control[END_REF][START_REF] Maciejowski | Predictive control: with constraints[END_REF]:

• The model of the system captures the dynamic and static interactions between input, output and disturbance variables.

• The physical constraints on inputs and outputs can be handled in a systematic manner.

• Multiple operational goals can be taken into account simultaneously.

• It is particularly suitable for those systems for which the disturbances can be forecasted.

The main principle of MPC resides in computing a sequence of inputs that causes the predicted response of the system to move to the desired setpoint in an optimal manner while respecting the constraints. The constraints imposed by the elements that make up the model have already been defined in Sections 3.1-3.4. On the other hand, the set of operational goals is defined below.

Operational goals and multi-objective function

One or more operational goals are expected to be achieved during the process. To this end, a certain criterion is optimized in the computation of the control signals. This criterion is usually built as the weighted sum of several terms, where each of them represents an operational goal.

Note that the set of operational goals that can be taken into account is not unique. In this work, the following are considered:

• Maintaining the water levels close to the setpoints: This is the most important objective to be fulfilled. Its mathematical formulation reads as

J 1 k = (y k -y r ) (y k -y r ) , (20)
with y r the vector of NNL values.

• Cost reduction: This term reflects the economic costs derived from operating the available equipment. It can be formulated as

J 2 k = γ u k u k , (21)
with γ the vector of known costs associated to the equipment operation.

• Smoothness of the control signal: In order to avoid wear and tear, and increase the lifespan of the equipment, it is a common practice to penalize the control signal variation between consecutive time instants:

J 3 k = ∆u k ∆u k , (22)
with ∆u k = u k -u k-1 .

• Penalty in the relaxation parameter: α k , which was introduced in (8), is penalized to ensure that the water levels are outside the navigation interval as little as possible:

J 4 k = α k α k . ( 23
)
The multi-objective function J that gathers the control objectives can be described by

J (u k , y k , α k ) = Hp k=1 4 j=1 β j J j k , (24)
where H p is the prediction horizon and β j are the weights of the j-th objective. Note that H p must be chosen according not only to the system dynamics (settling time), but also to take into account the system delays. Therefore, H p > t s + max S, where t s is the settling time (in samples), and S was defined for [START_REF] Litrico | Modeling and Control of Hydrosystems[END_REF]. Moreover, in order to set the weight of each objective in a multi-objective R3.8

optimization problem, the procedure described in [START_REF] Toro | Tuning of predictive controllers for drinking water networked systems[END_REF] can be used.

MPC formulation

Gathering the control-oriented model, the system constraints and the multi-objective function, the design of the MPC follows classical approaches [START_REF] Camacho | Model Predictive Control[END_REF][START_REF] Maciejowski | Predictive control: with constraints[END_REF]: an optimization problem is solved over a prediction horizon, minimizing a cost function while respecting the system constraints. The first component of the vector of control inputs is extracted from the solution and is applied to the system, and the rest are disregarded. This procedure is repeated at each time instant, following a R2.8

receding-horizon strategy.

The optimization problem is given by min

{u i|k } k+Hp-1 i=k , {y i|k } k+Hp-1 i=k , {α i|k } k+Hp-1 i=k J u i|k , y i|k , α i|k (25a)
subject to:

x i+1|k = Ax i|k + A n x i-n|k + B u u i|k + B un u i-n|k + (25b) B d d i|k + B dn d i-n|k , i ∈ {k, ..., k + H p -1},
y i|k = Cx i|k + C n x i-n|k + D u u i|k + D un u i-n|k + (25c) D d d i|k + D dn d i-n|k , i ∈ {k, ..., k + H p -1}, 0 = E u u i|k + E un u i-n|k + E d d i|k + E dn d i-n|k , i ∈ {k, ..., k + H p -1}, (25d)
u m ≤ u m i|k ≤ u m , i ∈ {k, ..., k + H p -1}, (25e)
y r -α i|k ≤ y i|k ≤ y r + α i|k , i ∈ {k, ..., k + H p -1}, (25f)
α i|k ≥ 0, i ∈ {k, ..., k + H p -1}, (25g)
x j|k = xMHE j , j ∈ {k -n, ..., k}, (25h)
u l|k = u M P C l , l ∈ {k -n, ..., k -1}, ( 25i
)
where k is the current time instant, i is the time instant along the prediction horizon and k + i|k indicates the predicted value of the variable at instant k + i using information available at instant R1.8

k. Remark: j and l are used to indicate the use of past information, for which the considered time intervals are different than the one described by i. Equations (25b) and (25c) correspond to the model described by ( 17), (25d) are the mass balances given in (4), and (25e)-(25g) are the constraints given in ( 1), ( 8) and ( 9), respectively. Equation (25h) sets the values of the delayed R2.10b

states according to the solution provided by the MHE (noted as xMHE i

) in past iterations. Note that the MHE will be introduced in Section 4.2. These delayed values are provided to the MPC as parameters. Similarly, the delayed control actions obtained by the MPC (noted as

u M P C i ) in
previous iterations are also provided as parameters by means of (25i).

The optimal solution is given by the sequences

{u i|k } k+Hp-1 i=k , {y i|k } k+Hp-1 i=k , {α i|k } k+Hp-1 i=k
1 . As R3.9

it was stated before, only u k|k is applied to the system, according to the receding philosophy

u M P C k u k|k . (26)

State estimation: the MHE approach

The control strategy presented in Section 4.1 uses the states to compute the set of optimal control actions. The system states oftentimes are not directly measurable, and therefore they need to be estimated from the available data using a state estimator.

Thus, the problem to be solved is that of designing an observer that fully reconstructs the system states. In this work, the MHE is used for this purpose. The main principle of this technique consists in formulating the estimation problem as a quadratic program using a moving estimation window of a fixed size [START_REF] Rao | Constrained linear state estimationa moving horizon approach[END_REF][START_REF] Allgöwer | Nonlinear predictive control and moving horizon estimation -an introductory overview[END_REF]. Indeed, it is assumed that only part of the available information of the system (inputs and outputs) is considered, which is shifted in time to consider the most recent information. Otherwise, the computational burden renders the full-information problem impractical to solve, as more and more data are processed with time. In this way, a truncated sequence of state estimates is computed at each time step instead of the full-state sequence to make the problem tractable [START_REF] Rao | Moving horizon strategies for the constrained monitoring and control of nonlinear discrete-time systems[END_REF].

The formulation corresponding to the optimization problem solved by the MHE reads as

1 {u i|k } k+Hp-1 i=k
{u k|k , u k+1|k , • • • , u k+Hp-1|k }; y i|k and α i|k are defined in the same manner

min {x i|k } k i=k-N xk-N|k -x k-N P -1 xk-N|k -x k-N + (27a) k-1 i=k-N w i|k Q -1 w i|k + v i|k R -1 v i|k
subject to:

w i|k = xi+1|k -Ax i|k + A n xi-n|k + B u u i|k + B un u i-n|k + (27b) B d d i|k + B dn d i-n|k , i ∈ {k -N, ..., k -1}, v i|k = y i|k -Cx i|k + C n xi-n|k + D u u i|k + D un u i-n|k + (27c) D d d k+i|k + D dn d i-n|k , i ∈ {k -N, ..., k -1}, 0 = E u u i|k + E un u i-n|k + E d d i|k + E dn d i-n|k , i ∈ {k -N, ..., k -1}, (27d)
x r ≤ xi|k ≤ x r , i ∈ {k -N, ..., k -1}, (27e) xj|k
= xMHE j , j ∈ {k -N -n, ..., k -N -1}, (27f)
u l|k = u M P C l , l ∈ {k -N -n, ..., k -1}, (27g)
y m|k = y m , m ∈ {k -N, ..., k -1}, (27h)
with (27b) accounting for the system disturbances and (27c) for the measurement noise. Additionally, (27d) describes the static part of the model, (27e) defines the valid interval of the state variables, and (27f), (27g) and (27h) set the values of the delayed states, inputs and outputs, respectively, following the same ideas as in [START_REF] Schuurmans | Modeling of irrigation and drainage canals for controller design[END_REF]. The value x k-N in (27a) corresponds to the most R3.10 likely initial state vector, and is chosen based on the available knowledge of the system, whereas xk-N|k is the first value of the optimal state sequence computed by the MHE at time instant k.

The error in this initial guess, given by xk-N|k -x k-N , is weighted by means of the matrix R3.11 P -1 , which indicates the confidence into the initial state, and its tuning allows to guarantee the boundedness of the estimation, as discussed in [START_REF] Rao | Constrained State Estimation for Nonlinear Discrete-Time Systems: Stability and Moving Horizon Approximations[END_REF]. On the other hand, Q -1 and R -1 are the weighting matrices inverses of suitable dimensions linked to the confidence in the quality of the model and the measurements, respectively. The larger these matrices are, the lesser the confidence in the associated term is, as the matrices are inverted. These inverses are directly related to the co-variance matrices only in the case of linear systems with zero-mean uncorrelated random variables for unknown disturbances [START_REF] Boegli | Real-Time Moving Horizon Estimation for Advanced Motion Control[END_REF]. In any other situation, e.g., constrained states, this connection is only an approximation.

The MHE problem ( 27) is formulated as follows: at the current time instant k, N input-output pairs [(u k-N , y k-N ) : (u k-1 , y k-1 )] shall be available. Therefore, N is the length of the moving estimation window, which bounds the size of the problem. The resulting least-squares problem is solved, yielding the optimal sequence {x i|k } k i=k-N2 .However, as is the case in the MPC problem, only one value in the sequence is considered, and the rest are discarded. In the MHE problem, this corresponds to the last value, that is, xk|k . Therefore,

xMHE k xk|k . ( 28
)
In the next iteration, for k = k + 1, the truncated data sequence is updated, and becomes

[(u k -N , y k -N : (u k -1 , y k -1 )], which is equivalent to [(u k-N +1 , y k-N +1 ) : (u k , y k )]. Then, note
that the oldest measurement pair (u k-N , y k-N ) is dropped, and the newest measurement pair (u k , y k ) is incorporated, following the moving horizon philosophy.

Simulation

Once the MPC and MHE are designed, they must be integrated in the simulation loop. The solution of the controller (the optimal control inputs), together with the measurements, are fed into the estimator. In turn, the estimator computes the optimal state estimates, which are used in the next time instant by the controller to compute the new set of control inputs.

The real system is equipped with sensors, which provide water level measurements. Unfortu-R1.9

nately, these real values are not available in this work. Since the MHE requires the measurements to estimate the states, these must be generated in simulation, using the output equation (17b).

The effect of this limitation is that the estimator cannot be used for the first time at k = 1, but at k = N + 1. Indeed, an input-measurement pair (u k , y k ) will be generated at each time instant k, using the solution of the controller. Thus, the necessary data to compute the state estimates will not be available until N samples have elapsed. In addition, the system must be manually R2.11

initialized by selecting any feasible state vector.

Finally, note that, at time k, the MPC yields the sequence {u i|k } k+Hp-1 i=k

. Thus, it is necessary to know the disturbances until, at least, the instant k + H p -1. By contrast, the MHE sets its starting point N samples in the past, and reconstructs the optimal sequence of state estimates until the current time instant k. This requirement in terms of available information is fulfilled based on the policy introduced in Section 3.2, which allows to anticipate future lock operations.

Case study

This section illustrates the performance of the MPC and the MHE by means of a realistic case study system, which is based on part of the inland waterways in the north of France. First, the system is described, emphasizing its physical features, so that the modeling step yields a model as close as possible to the real system. Then, the experimental design step is regarded, describing the considered scenarios, which try to represent faithfully an average navigation day. Finally, the results coming from the experimental design are presented and discussed.

System description

The inland waterways in the north of France is linked with the Belgian and Dutch inland waterways, and is managed by Voies Navigables de France3 (VNF). Its main objective is that of guaranteeing the navigability condition, which is achieved by keeping the water levels inside the navigation rectangle defined by the LNL and the HNL, and as close as possible to the NNL.

This inland navigation network consists of more than fifty reaches that are interconnected by locks, gates and weirs. Part of it is depicted in Fig. 4, which shows the two reaches considered in the case study. The i-th reach is labeled as NR i , and its setpoint (NNL) is specified in red. In R3. [START_REF] Vukov | Real-time nonlinear MPC and MHE for a large-scale mechatronic application[END_REF] addition, the locks that connect adjacent reaches are labeled in black.

A more schematic view of the system is depicted in Fig. 5, resulting in a four-reach case study.

The case study choice is motivated by the following reasons:

• It features a distributary, which branches off from NR 1 at an intermediate point and flows away, to the lock of Don. This topology is regarded as of special interest, since the mass balance at this natural bifurcation (not controlled) is not straightforward to model. Indeed, a possible approach for this situation is shown below. • Besides, NR 2 is an important reach in this network for two reasons [START_REF] Segovia | Modeling and fault diagnosis of flat inland navigation canals[END_REF]: its strategic location, which allows dispatching water among the three major catchments in the region; and its downstream lock in Fontinettes, which performs the largest lock operations in terms of volume, and is therefore responsible for the largest disturbances. Being able to deal with R3.14 the worst-case scenario can give a feel for the magnitude of the disturbances that the control strategy attempts to reject.

Since the bifurcation is of natural type (uncontrolled), this node can be eliminated, based on an estimation of the ratios of the total flow for each stream after the bifurcation. Indeed, it can be considered, according to VNF, that the each of the flows after the bifurcation correspond to 50% of the flow before the bifurcation. This yields the simplified, final three-reach case study scheme given in Fig. 6. Note that the reaches are renamed for convenience, and also the nodes for labeling purposes.

The physical parameters of the case study are summarized in Table 1. Note that the lengths of the reaches are approximately the same, which results in the same time delay (in samples) for which is deemed reasonable due to the slow dynamics of the system. This allows to use model [START_REF] Segovia | Input-delay model predictive control of inland waterways considering the backwater effect[END_REF], and thus there is no need to resort to the more general case [START_REF] Litrico | Modeling and Control of Hydrosystems[END_REF].

Table 1: Parameters of the reaches LNL, NNL and HNL are the relative lower, normal and higher navigation levels (with respect to the bottom of the reach), L is the length of the reaches, w r is the bottom width, m r is the side slope of the reach (m r = 0 for a rectangular cross section), s b is the bottom slope (s b = 0 for a flat reach), n r is the Manning roughness coefficient and Q s is the operating point considered when linearizing the Saint-Venant equations. Indeed, it is considered that an average flow of 1.2 m 3 /s comes from upstream of Douai, and that it is divided into two equal parts after it. The numerical state-space matrices computed from these parameters are given in Appendix A.

Reach nr. LNL [m] NNL [m] HNL [m] L [m] w r [m] m r [m/m] s b [m/m] n r [s/m 1/3 ] Q s [m 3 /s] (1,2) 3.95 4
In order to estimate the states, the measurements of the system are needed. Depending on the structure of the problem, e.g., the topology of the network or the connections between reaches, not all measurements are required. This statement can be realized by inspecting matrix C, which links the states and the measurements. In the present case, it is not necessary to consider the measurement y

C(1,2) k
, since the associated state has an effect on the downstream level y

C(4) k
, given by the off-diagonal, nonzero entry in the fourth row of C. Thus, in order to show the effectiveness of the approach, the six states will be reconstructed with only five measurements, i.e., assuming that the water level y

C(1,2) k
is not available. Thus, this value will be obtained from the state estimates.

The lock operations that take place in Douai, Don, Cuinchy and Fontinettes disturb the system.

Their average magnitudes and durations are given in Table 2. On the other hand, the same four nodes are equipped with controlled devices which allow to dispatch water to fulfill the control objectives. In particular, Douai and Fontinettes are equipped with undershot gates, whereas Don and Cuinchy are equipped with weirs. It is considered that both the gates and the weirs can deliver a maximum flow of 10 m 3 /s, which will have to be converted into maximum gate openings and sill elevations, respectively. Their nonlinear expressions, as well as their linearized equations, are given in Appendix B.

Experimental design

A 24-hour scenario, depicted in Fig. 7, is designed by considering a lock operation time-series model for a typical navigation profile. In the real system, the following existing management restrictions must be taken into account:

• A day is divided in two periods: navigation and stoppage. Boats are only allowed to navigate during the navigation period, which starts each day at 6 a.m. and finishes after fourteen hours, at 8 p.m. The navigation is interrupted until the next day at 6 a.m.

• The current policy allows a maximum of two lock operations per hour.

Besides, the scenario does not consider changes in the setpoints, thus assuming that the navigation conditions do not change during the simulation. Such modifications typically occur due to changes in the weather condition, e.g., flood and drought episodes, which might require to readjust the LNL, NNL and HNL values.

In order to test this scenario, it is necessary to compute the IDZ model of each reach as shown in [START_REF] Litrico | Modeling and Control of Hydrosystems[END_REF], and then discretize them using T s . The global model ( 17) is built by stacking the IDZ model of each reach, so that the MPC and MHE can be used together. interval [6127, 6807] s, or equivalently [5.11, 5.67] samples. Thus, a unique delay of 6 samples can be considered for all the reaches, according to the ceiling rule introduced in [START_REF] Breckpot | Flood control of rivers with nonlinear model predictive control and moving horizon estimation[END_REF]. In this way, by setting the prediction horizon values equal to 12 samples, the waves are given enough time to reflect at the downstream end and travel back.

Results

As mentioned before, any feasible initial state might be considered, as there is no available data that can be used to estimate the initial conditions. Once the system reaches the steady state, it is disturbed by considering the 24-hour navigation profile depicted in Fig. 7. The final state-space model is computed, and then the results yielded by the estimator and the controller are shown.

Control-oriented model

For the sake of convenience, the steps are outlined here, while the final system matrices as well as the nonlinear and linearized equations are detailed in the Appendices A and B, respectively.

First, the methodology described in [START_REF] Litrico | Modeling and Control of Hydrosystems[END_REF] must be applied for each reach, and then these partial models [START_REF] Joseph-Duran | Output-feedback control of combined sewer networks through receding horizon control with moving horizon estimation[END_REF] are stacked to build the global model. Next, the nonlinear equations of the gates and weirs are linearized around the NNL and Q s values given in Table 1. The linearized equations must be substituted in the model. Note that q D does not depend on its upstream water elevation.

Since this water level is outside of the scope of the control problem, it is considered that there is enough water upstream of Douai, and thus that this level remains constant.

Next, it is necessary to substitute the water depths by the states in the linearized equations.

The relationship between these variables provided by matrix C can be exploited to this effect.

Once this step has been completed, the linearized equations can be substituted in the model. As it was mentioned in Section 3.5, the delayed expressions q k-n = f x k-n , u k-n introduce the delayed states in the final formulation. Note also that, since q k depends on the states, a rearrangement and grouping of terms is required. The final model is then given by ( 17), with

x k+1 = x D(1,2) k+1 x C(1,2) k+1
x D(1,3) k+1

x Don k+1 x

C(4) k+1 x F k+1 , x k = x D(1,2) k x C(1,2) k x D(1,3) k x Don k x C(4) k x F k , x k-n = x D(1,2) k-n x C(1,2) k-n x D(1,3) k-n x Don k-n x C(4) k-n x F k-n , u k = u D k u C k u Don k u F k , u k-n = u D k-n u C k-n u Don k-n u F k-n , d k = d D k d C k d Don k d F k , d k-n = d D k-n d C k-n d Don k-n d F k-n , y k = y D(1,2) k y D(1,3) k y Don k y C(4) k y F k .
Since the bifurcation node can be eliminated as proposed in Fig. 6, (4) is no longer needed.

Estimator

The comparison between the optimal estimated states given by the MHE and the real states obtained in simulation are depicted in Figure 8. At each time instant, the MHE is fed with the corresponding sequence of N input-measurement pairs, which are used to compute the optimal sequence of states.

In general, the values provided by the estimator match the real states with no significant error.

In addition, the real states are noisier than the estimated states. This behavior is in line with the nature of the observer, which acts as a filter, smoothing the predictions. Furthermore, the constraints on the state bounds are satisfied. Such bounds are not even depicted in Fig. 8 for the benefit of a better visualization of the results, as the states are far from the bounds.

In order to ensure a quantitative comparison between the real and the estimated states, the similarity of both signals is quantified by means of the correlation coefficient. Given a pair of signals (m t , n t ) with M observations each, the correlation coefficient is defined as

ρ m,n = 1 M -1 M i=1 m i -µ m σ m n i -µ n σ n , (30)
where µ m and σ m are the mean and standard deviation of m t , respectively, and µ n and σ n are the mean and standard deviation of n t . This coefficient is bounded between 1 and -1: the closer this coefficient is to 1 (respectively -1), the stronger the positive (respectively negative) correlation between the pair of signals is. The correlation coefficients between the real and the estimated states are summarized in Table 3. It can be stated that the performance of the MHE is satisfactory, since all the correlation coefficients are very close to 1, which indicates a strong, positive correlation. Indeed, the main goal of the estimator is that of reconstructing an accurate state sequence that is not directly measurable. This procedure is used to achieve the final goal of this work, which consists in fulfilling the control objectives, so that the desired system performance is attained. Therefore, state estimation is regarded as a tool employed by the controller in pursuit of the final goal.

Controller

The estimated states are used by the controller in order to compute the sequence of future optimal inputs, applying only its first component. It must be recalled that the real system measurements are not available in this work, and thus they must be obtained using the output equation.

Based on the control objectives defined in Section 4.1.1, two main results are looked at: the water levels and the control signals, depicted in Figs. 9 and 10, respectively. The modeling simplification introduced in Section 5.1 led to considering three reaches, where their upstream and downstream water levels are arranged by columns in Fig. 9. The flow diversion is considered to take place downstream of the gate of Douai, and thus the water levels in Douai for reaches [START_REF] Mihic | Policy and promotion of sustainable inland waterway transport in Europe-Danube River[END_REF][START_REF] Van Overloop | Predictive control for national water flow optimization in The Netherlands[END_REF] and (1,3) might be different. Recall also that y C(1,2) k is not measured, but computed from the state estimates.

Regarding the control objectives linked to the water levels, it can be stated that the MPC is able to keep the levels very close to the setpoints despite of the disturbances. To quantify the

T P = 1 - 1 H p Hp k=1   y k -y r 1 2 y r -y r   2 . (31)
Equation ( 31) was introduced in [START_REF] Segovia | Input-delay model predictive control of inland waterways considering the backwater effect[END_REF] as a modification of the standard root relative squared error. These tracking performance indices are defined as the relative error between the predicted levels y k and the setpoints y r (NNL values). The denominator, given by 1 2 y r -y r , equals the semi-amplitude of the symmetric [LNL, HNL] interval, which is the maximum allowed variation from y r . The squaring emphasizes larger differences, which is of interest in this case, since it focuses on the water levels y k that are far from the setpoints y r .

The numerical values of the indices for each water level are summarized in Table 4. Note that TP is bounded between 0 and 1, where 1 corresponds to the perfect tracking performance.

Therefore, it can be stated that the MPC provides satisfactory results in terms of keeping the water levels close to the setpoints. 9 shows that the water levels are never outside of the navigation rectangle, and therefore the penalty on this behavior, represented by α α α k in (23), equals 0.

On the other hand, the control objectives ( 21) and ( 22) are linked to the control signals, depicted in Fig. 10. The undershot gates in Douai and Fontinettes work at their maximum capacity for a maximal gate opening, whereas the weirs work at their maximum capacity when the sill elevation is zero. This can be realized by the minus sign in the linearized weir equations. When the water levels reach the setpoints, the flows delivered by the actuators should be minimum, i.e., minimum gate openings and maximum sill elevations. However, rejecting disturbances results in continuously operating the actuators, as these disturbances go against the control objectives. Therefore, it can be seen how the actuators work close to their maximum capacity only during short periods of time, and always within the equipment design range. During most of the simulation, the inputs are far from the physical limits of the actuators, thus taking into account the cost reduction objective [START_REF] Litrico | Conversion from discharge to gate opening for the control of irrigation canals[END_REF].

Regarding the smoothness of the control signals given by the operational goal [START_REF] Malaterre | Regulation of irrigation canals[END_REF], the control actions present some peaks. However, there are no large variations between consecutive control actions, especially compared with the design range of the actuators. This behavior should result in a long lifespan of the equipment. Although the weight of this objective could be increased, this would probably interfere with the rest of objectives.

Conclusions

This work regarded the design of state estimation and control techniques in the framework of inland waterways, often characterized by negligible bottom slopes and large time delays. To this end, the MHE and MPC approaches were considered, respectively. This choice obeys to the appropriateness of these optimization techniques to solve the considered problem. In addition, the existing duality between both can also be exploited. Even though this work targets inland waterways as the field of application, for which the time-delay structure comes directly from the control-oriented modeling formulation, the general methodology can be applied to any kind of system characterized by the transportation of mass, energy or information.

The main objective of this work was that of fulfilling a set of operational goals, each of them linked to an aspect of the system performance. This task was taken care of by the MPC, which requires the knowledge of the states in order to determine the optimal predictive control law.

However, as the states are not directly measurable in this problem, the MHE was used for this purpose. Both techniques were synchronized, using the previously derived control-oriented model formulation. Finally, the inland waterways in the north of France were used to build a realistic example in order to demonstrate the performance of the methodology. Although the simulation results are deemed satisfactory, the centralized architecture employed in this work can encounter implementation problems due to the large dimensionality of the network. Thus, non-centralized implementation approaches might be considered in the future to deal with this limitation, possibly continuing the first steps carried out in [START_REF] Segovia | Distributed Input-Delay Model Predictive Control of Inland Waterways[END_REF].

From the proposed approach, it has been shown that the MPC is able to reject the disturbances that are caused by the lock operations. These disturbances are assumed to be known in this work, but this seldom happens in real applications. Thus, strategies to estimate the effect of unknown disturbances, such as the unknown input observer (UIO), might be considered in the future. This class of observers assumes no a priori knowledge about such inputs, which is an interesting feature to be exploited for fault diagnosis purposes. Moreover, combining fault diagnosis and non-centralized control and state estimation could lead to implementable reconfiguration strategies, aiming at guaranteeing that the system exhibits an acceptable performance even in the presence of faults.

Both MPC and MHE use a dynamic model of the process. In this work, the IDZ model was used as the starting point to build a control-oriented model. The IDZ is a physical model, and is originally formulated for a canal with two inputs and two outputs. This formulation was extended to consider larger portions of inland waterways that comprise more than one canal. To do so, the original model was manipulated, yielding a more convenient state-space representation, which also allows to handle the delays in the system in a more suitable manner. This work tackled the particular case in which the lengths of all the reaches are approximately equivalent, leading to only one delay in the network. In addition, this work considers only one operating point during the whole simulation. Although this is sufficient when the setpoints do not change during the simulation, this assumption is rather limited in the case of large operating ranges. In this regard, LPV models and Takagi-Sugeno fuzzy models could be considered. Remark: The delay of the system is equal to 6 samples, as stated in Section 5.2.

A =              

B. Equations of the gates and weirs B1. Nonlinear equations

The nonlinear equations that describe the dynamics of the gates and the weirs in the considered case study are:

Undershot gate: q = C dg L g u 2gy 1

(32a)

Weir: q = C dw 2g (y 1 -u) 3/2 (32b)

where q is the discharge through the structures; C dg and C dw , the discharge coefficients, equal to R3.15 0.6 and 0.4, respectively; L g , the gate width; u, the gate opening or elevation; g, the acceleration of gravity; and y 1 , the upstream water level.

B2. Linearized equations

Using the NNL and Q S values in Table 1, the linearized versions of ( 32

2 Figure 1 :

 21 Figure 1: Water level regulation with: (a) an undershot gate. (b) a sharp crested weir.

Figure 2

 2 Figure 2 depicts a navigation canal consisting of three reaches, separated from one another by cross structures. The red solid arrows represent uncontrolled inputs and outputs such as natural bifurcations. On the other hand, the green dashed arrows indicate the flows generated due to the lock operations. Finally, the blue dotted arrows represent controlled actions, carried out by gates and weirs in order to regulate the water levels. Notice that the locks are often built next to a control structure.

Figure 2 :

 2 Figure 2: Navigation canal and its water-resource exchanges

Figure 3 :

 3 Figure 3: Navigation canal schematic with the variables involved

  and f and g are the corresponding relationships among the variables.

Figure 4 :Figure 5 :

 45 Figure 4: Part of the inland waterways in the north of France

Figure 6 :

 6 Figure 6: Simplified schematic diagram of the case study

Figure 7 :

 7 Figure 7: The considered lock operation profile

3 ]Figure 8 :

 38 Figure 8: State estimates (blue solid lines) and computed states (dash-dot gray lines)

Figure 9 :

 9 Figure 9: Water levels (blue solid lines), NNL (black dashed lines), and LNL and HNL (red dashed lines)

Figure 10 :

 10 Figure 10: Gate openings (blue solid lines) and physical limits (dash-dot red lines)

  values for the MPC problemβ 1 = 20 ; β 2 = 2 ; β 3 = 5 ; β 4 = 1.D. Weighting values for the MHE problemP -1 = I nx = I 6 ; Q -1 = I ny-1 = I 5 ; R -1 = I nx = I 6 .

  and A, A n , B u , B un , B d , B dn , C, C n , D u , D un , D d and D dn are time-invariant matrices of suitable dimensions. The state equation is given by (17a), and (17b)

Table 2 :

 2 Lock operations

	Lock	Dispatched water volume [m 3 ] Duration [min]
	Douai	18000	20
	Cuinchy	12000	20
	Don	12000	20
	Fontinettes	30000	20

Table 3

 3

	: Correlation coefficients
	D (1,2) C (1,2) D (1,3)	Don	C (4)	F
	0.8389 0.8856 0.8841 0.9276 0.9654 0.9403

Table 4

 4

	: Tracking performances
	D (1,2) C (1,2) D (1,3)	Don	C (4)	F
	0.9930 0.9877 0.9904 0.9907 0.9693 0.9717
	Furthermore, Fig.

{x i|k } k i=k-N {x k-N |k , xk-N+1|k , • • • , xk|k }

http://www.vnf.fr

Acknowledgments

This work has been partially funded by the Spanish State Research Agency (AEI) and the European Regional Development Fund (ERFD) through the projects DEOCS (ref. MINECO DPI2016-76493) and SCAV (ref. MINECO DPI2017-88403-R). This work has also been partially funded by AGAUR of Generalitat de Catalunya through the Advanced Control Systems (SAC) group grant (2017 SGR 482).

The authors want to thank VNF for the support and information provided regarding the inland waterways management.

Appendix

All the numerical values used in this work are presented here for the purpose of reproducibility.

A. Numerical model

The final state-space matrices are: