Gray-Box Model of Inland Navigation Channel: Application to the Cuinchy–Fontinettes Reach

2.1 Adaptive and Predictive Control Architecture

The main management objective of the inland navigation networks is to guarantee the seaworthiness requirement of each reach, i.e., the NNL (see Figure 1). These levels are principally disturbed by the navigation and the lock operations. During lock operations, large volume of water is withdrawn from the upstream reach and supplied to the downstream reach, causing a wave to travel to both directions: upstream to downstream, and downstream to upstream after reflection. Generally, these waves are not significant. However, sometimes, due to the dimension of the locks, the amplitude of the wave exceeds several tens of centimeters. This wave is reflected at the downstream end of the channel and travels several times back and forth. To reduce the effect of wave and to maintain the NNL, it is necessary to control the gates that are generally located beside the locks. Another possibility is to control the discharges from natural rivers. The water levels are controlled by gates and measured by tele-operating sensors. The management rules have to consider other strong constraints during drought and wet periods. During drought periods, the water for industry, agriculture, ecological aims, and drinking water has higher priority than the water for navigation. Some restrictions on navigation can be imposed during these periods. The navigation is managed by the planning of ship conveyances. However, the locks and gates have to be operated to maintain a minimal water flow (ecological flow) in order to decrease the water temperature, increase the oxygen rate in the water, and limit the modification of the ecosystem, such as the development of algae. During wet periods, some management protocols have to be applied. According to the estimation of the company responsible for the management of the French waterways, the navigation reaches are supplied by natural rivers even if the locks and gates are closed. Thus, before the occurrence of wet events, the level of each navigation channel is decreased, and their gates and locks are closed. The gates or the locks are operated only when the level of the reaches attained a critical limit beyond which the banks can be damaged. The amplitude and the frequency of extreme scenarios of drought and flood should be increased in the future as a consequence of global change. It is thus necessary to anticipate the effects of these events on the navigation reaches, by a study of their resilience and by the proposal of improved management strategies. To achieve these aims, an adaptive and predictive control architecture has been proposed in ref. [11], as depicted in Figure 2. This architecture is based on a SCADA system allowing the tele-control of the navigation network. A human–machine interface is dedicated to the supervision of the inland navigation network by a supervisor. The management constraints and rules are gathered in the management objectives and constraints generation module. To perform the management of the inland navigation network, a hybrid control accommodation module allows the determination of set points (management strategies block) according to the current state (supervision block) and the forecasting of the future state (prognosis block) of the network. These strategies can be adapted or improved according to the decision support (DS) module. This module is dedicated to test the resilience of the network against extreme scenarios and to design adaptive management strategies. In a first step, this module is used offline. Realistic scenarios of extreme events are built according to past events and to the knowledge on global change consequences. They are tested on a simulator of inland navigation networks. Thus, the DS module requires models of the dynamics of the inland navigation networks.

Figure 1.

Normal Navigation Level Inside the Navigation Rectangle.

Figure 2.

Adaptive and Predictive Control Architecture.

2.2 Inland Navigation Reach Modeling

The modeling of inland navigation reaches is necessary to study their resilience, to design adaptive and predictive management strategies, control algorithms, and FDI techniques such as proposed in Refs. [12, 18]. A navigation reach belongs to the class of open-channel systems with the particularity of very mild slope. By considering the operating points corresponding to the NNL, the dynamics of a reach can be modeled by a linear model. The structure of the model is selected in order to be suitable for the design of control strategies, FDI algorithms, etc., and particularly for MIMO systems. The proposed modeling approach assumes a first-order plus time-delay structure for every input–output pair. The output variables are chosen as the n_y measured water levels at time instant k L_i(k), i∈1, n_y in several selected points of the system (the variables and scalars are denoted in italic, the matrices in bold and italic). The input variables of the model correspond to the n_u input–output discharges Q_l(k), l∈1, n_u (see Figure 3). τ^y represents the time delay between one measurement point, L_i, and another measurement point, L_j, for all measurement points; and τ^u represents the time delay between one input discharge point, Q_i, and another input discharge point, Q_j, for all input discharge points.

Figure 3.

Time Delays τi, jy between Each Measurement Point.

Finally, time-delay matrices τ^y and τ^u are considered for each output and input variables of the model to take into account the transfer delays of the reach, leading to the following discrete time model:

where matrices A∈ℝny × ny.ny and B∈ℝny × nu.ny are the input and output matrices, respectively, and k is the current discrete time. At each time k, the vector u¯k|τu∈ℝnu.ny represents the input variable defined according to the delay matrix τ^u, and the vector y¯k|τy∈ℝny.ny represents the delayed output variables according to the delay matrix τ^y. The matrices τy∈ℕny×ny and τu∈ℕnu×ny gather the time delays between each measurement point (see Figure 3), and between each measurement point and each input and output of the system, respectively. For example, the value of τi, jy∈ℕ is the time delay between the measurement points L_i and L_j:

The time delays τi,  jy are equal to 0 for i=j.

Similarly, the time-delay matrix between the inputs can be written as

Just as in case of the measurement delay matrix, the time delays τi, ju are equal to 0 for i= j. The time-delay matrices τ^y and τ^u are supposed to be known and constant.

The elements of the matrices τ^y and τ^u can be obtained a priori by correlation methods (data-based procedure) or from the physical knowledge of the system. Consider two points along the canal, separated by a distance D: one located upstream and the other downstream. According to ref. [19], the theoretical value of the upstream time delay between these two points is evaluated by computing the integral:

with c(l) and v(l) representing the celerity and the velocity, respectively. This corresponds to the minimum time required for a perturbation to travel from the upstream point to the downstream one. Analogously, the downstream time delay τ_{d, u} can be evaluated by

and corresponds to the maximum time required for a perturbation to travel from the downstream point to the upstream one. In both cases, we recover the classic value in the uniform case when v and c are constant: τu, d = Dc + v and τd,  u = Dc − v, with c = gSb and v = QS. The cross-sectional area S and the bottom width b varies from upstream to downstream in the canal reach. Here, constant values are used. These constant values are obtained as the mean values of the different bottom widths and cross-sectional areas, respectively. g is the gravity and Q is the average flow of the canal. The vector y_k in eq. (1) is

Thus, the vector y¯k|τy is built according to the matrix τ^y such as

The vector u¯k|τu is built according to the matrix τ^u similarly as relation (7):

The input and output matrices are defined as

After the structure of the model has been defined, the objective is to identify the matrices A and B. Model (1) can be rewritten as

with M=[A B] and Φk = [y¯k|τy u¯k|τu]T. Then, the matrix M has to be determined according to an identification approach using available measured data. The data correspond to N samples of the discharges Q_i and levels L_i measured on a time interval. On the basis of relation (11), the matrix M is estimated in the following way:

with Y = [y_x+1 … y_N], Φ¯ = [Φχ ⋯ ΦN−1], and χ = max(τ^y) + 1, where max(τ^y) is the maximum entry of the matrix τ^y [see relation (3)].

Considering the characteristics of the matrices A and B, and the size of Y and Φ¯, it is easier to identify separately each line of M, before rebuilding the global matrices A and B. The zeros of matrices A and B are not considered during the identification step. It is based on the principle that each output variable can be expressed such as the first output variable:

Using this particular structure of the model, the dynamics of a real navigation reach can be identified. The CFR is presented in the next section.

3.1 Presentation

In the north of France, the inland navigation network allows the navigation of broad gauge ships from the regions in the north of Paris to the port of Dunkerque and to Belgium. This network covers three watersheds: the Aa, Lys, and Scarpe watersheds (see Figure 4). Owing to the topography, the canalized rivers and channels are separated by locks that allow the navigation. The CFR is located in the center of the network, between the upstream lock of Cuinchy at the east of the town Bethune and the downstream lock of Fontinettes at the southwest of the town Saint-Omer. The CFR has a major importance for the management of this inland navigation network. The first part of the CFR, i.e., 28.7 km from Cuinchy to Aire-sur-la-Lys, is called “canal d’Aire” and has been built in 1820 (see Figure 5). The second part of the channel, i.e., 13.6 km from Aire-sur-la-Lys to Saint-Omer, is called “canal de Neuffossé” and has been built in the 11th century. The direction of the flow is from Cuinchy to Fontinettes. The CFR is entirely artificial without a significant slope and with 614 different transversal cross sections that have to be considered. The CFR is managed by Voies Navigables de France, whose role is to maintain the level of the channel at NNL=19.52 NGF (Nivellement Général de la France, i.e., altitude landmark in France). To reach this aim, three points of the CFR must be controlled:

• At Cuinchy: a lock and a gate that are located side by side

• At Aire-sur-la-Lys: the gate called “Porte de Garde”

• At Fontinettes: a lock

Figure 4.

Inland Navigation Network in the North of France.

Figure 5.

Scheme of the Cuinchy–Fontinettes Navigation Reach.

The Cuinchy lock is composed of a 2-m-high chamber with a volume of 3700 m³. The maximum discharge that supplies or empties the chamber is equal to 11 m³/s. The chamber of the Fontinettes lock is 13 m high with a volume of 25,000 m³, and a maximum discharge equal to 30 m³/s. The control of the Cuinchy and Fontinettes locks is constrained by the navigation demand. The valves that supply and empty the lock chambers are controlled by the lock keeper according to several conditions: number of ships, size of ships, etc. Thus, the time of a lock operation, i.e., the discharge, is never the same. In any case, the operation of the locks causes a wave phenomenon that influences the CFR. In particular, the operation of the Fontinettes lock can cause a wave with an amplitude that exceeds 13 cm. This wave is attenuated only after 2 h.

3.2 CFR Modeling on SIC Software

The model proposed in ref. [3] is based on the Saint-Venant equations. It is a simplified model of the CFR that is dedicated for controller design: to develop a controller using that model to limit the impact of this wave. In ref. [10], the proposed gray-box model is identified with data from the SIC-1D hydrodynamic model [22] by considering the CFR such as a channel with a rectangular profile. In this article, the 614 transversal cross sections are considered. To implement the geometry of the CFR, measurement data are needed. The cross-section data are stored as x and z coordinates: x containing the horizontal distance and z is the elevation. An example of the cross section to be implemented using SIC is shown in Figure 6. These cross sections contain several data points; however, to implement them using SIC, the limit is 24. Thus, the number of points in each cross section has been reduced to 24 points. Altogether, 614 cross-section files have been implemented using SIC. Then, the locks and gates are considered as the inputs and outputs of the SIC model. The operation of the locks corresponds to a discharge during a time, with a trapezoidal pattern as shown in Figures 7 and 8 for Cuinchy and Fontinettes, respectively. Finally, SIC is used to reproduce the dynamics of the CFR, particularly the wave phenomenon. However, SIC is a deterministic model. It requires that all the data are known. In the real case, all the data are, in general, not available due to unknown inputs or outputs, noise, etc. The gray-box modeling presents the advantage to deal with some uncertainties. In this article, the gray-box model is valuated using data from SIC.

Figure 6.

Example Cross Section of the CFR. The Locations of the Measurement Are Shown with Dots, and the Cross Section Is Interpolated (Straight Line) between Them.

Figure 7.

Cuinchy Lock Operations for Different Discharge Schedules. With Black, an Operation with Maximum of 11 m³/s Discharge, with Red 9 m³/s, with Blue 7 m³/s, and with Magenta 5 m³/s.

Figure 8.

Fontinettes Lock Operations for Different Discharge Schedules. With Black, an Operation with Maximum of 30 m³/s Discharge, with Red 25 m³/s, with Blue 20 m³/s, and with Magenta 15 m³/s.

3.3 Gray-Box Modeling

The CFR is a MIMO system with three inputs, Q₁ for Cuinchy, Q₂ for Aire-sur-la-Lys, and Q₃ for Fontinettes, i.e., n_u = 3. The levels of the CFR are measured at three points: L₁ at Cuinchy, L₂ at Aire-sur-la-lys, and L₃ at Fontinettes, i.e., n_y = 3. The time delays are determined according to relations (4) and (5), with the mean cross-sectional area S and the maximum discharge Q_max. The mean cross-sectional area S is computed according to the mean bottom width equal to 52 m and the water depth that corresponds to the NNL, i.e., 4.26 m. The maximum discharge is evaluated according to the Fontinettes lock operation and is equal to Q_max = 7.6 m³/s. The matrices τ^y = τ^u (in other words, the measurement points are located in the same place as the inputs) are given in minutes:

The vector y_k is

and the vector y¯k|τy is given by

The vector u¯k|τu is

In case of this model, τ1, 1u, τ2, 2u, and τ3, 3u are equal to zero. The modeling approach consists in identifying the coefficients of the output and input matrices defined as

3.4 Identification of the CFR Gray-Box Model

The identification task consists in estimating the coefficients of the matrices A and B according to Eq. (12). To achieve this aim, the model of the CFR that is implemented using SIC is used to generate data. The proposed scenarios consist in the simulation of the lock-and-gate operations. The first scenario, entitled “scenario 1,” is used for the identification of the matrices A and B. It corresponds to 2 days of navigation with Cuinchy and Fontinettes lock operations. During these 2 days, the CFR was crossed by 10 ships. The time of the lock operations is different (see Figure 9A). The CFR is also supplied by the gate at Aire-sur-la-Lys (see Figure 9A, with red line). This discharge allows to maintain the NNL despite the Fontinettes operations. The Cuinchy lock operations cause a wave phenomenon with a maximum amplitude of 8 cm (see Figure 9B). The operations of the Fontinettes lock cause a wave phenomenon along the CFR with an amplitude maximum equal to 13 cm (see Figure 9D). This wave phenomenon is propagated along the CFR (see Figure 9C). These signals obtained according to the SIC software are used to identify the model of the CFR. Owing to the properties of the signals, the matrix (Φ¯ Φ¯T) of eq. (12) can be singular. Thus, it is often necessary to use the Moore–Penrose pseudoinverse of this matrix. Finally, the coefficients of the matrices A and B are given by Tables 1 and 2, respectively. After the identification step, the outputs of the model y^k and the measured levels y_k = [L₁(k)L₂(k)L₃(k)]^T are compared as depicted in Figure 10. The estimated outputs in red dashed line are very close to the measured levels in blue continuous line. To estimate the effectiveness of the model, a fitting indicator (FIT) is defined as FIT = (1 − 1N − χ∑||y^ − y||2||y − ym||2) × 100, where y_m is a column vector composed of the mean value of y. For this scenario, the FIT is given for each level in Table 3. The effectiveness of the CFR model has to be evaluated in other scenarios. These tests are presented in the next section.

Figure 9.

Scenario 1. Discharges at (A) Cuinchy Q₁ (Blue Line), Aire-sur-la-Lys Q₂ (Red Line), and Fontinettes Q₃ (Black Line). Levels at (B) Cuinchy L₁, (C) Aire-sur-la-Lys L₂, and (D) Fontinettes L₃.

Figure 10.

Scenario 1. Estimated (Red Dashed Line) and Measured (Blue Continuous Line) Levels at (A) Cuinchy L₁, (B) Aire-sur-la-Lys L₂, and (C) Fontinettes L₃.