Agriculture plays an important part in the food chain and water resources for agriculture are essential. A problem is that the water transport systems present low efficiencies in practice. Crop yields must be optimized, and the goal of an operational water manager is to deliver water to irrigation sites accurately and efficiently. In order to fulfill this objective, we propose a centralized overall control diagram to optimize the management of the canal. Our control diagram in real time is mainly composed of two algorithms, CSE and GoRoSoBo. The first one is a powerful tool in canal management, and is able to estimate the real extracted flow in the canal and the hydrodynamic canal state from measured level data at selected points. The second one is an essential tool in the management of the canal, a feedback control algorithm operating in real time. The GoRoSoBo algorithm (Gómez, Rodellar, Soler, Bonet) is able to calculate the optimum gates trajectories for a predictive horizon taking into account the current canal state and the real extracted flow, both obtained by CSE.

## INTRODUCTION

The world population is predicted to grow from 6.9 billion in 2010 to 8.3 billion in 2030 and 9.1 billion in 2050 (UNDESA 2009) and food demand is predicted to increase by 50% in 2030 and by 70% in 2050 (Bruinsma 2009).

The most recent estimate for irrigated agriculture is an increase from 2,743 km^{3} in 2008 to 3,858 km^{3} in 2050 (FAO 2011a, 2011b). Much of the increase in irrigation water will be in regions already suffering from water scarcity. To study the impact of water scarcity, there are several simulation models available, for instance, MCG (general circulation model) or MCR (regional climate model). Most specifically in Spain, some of these models (AR4 (Rossby Centre regional atmospheric model) and HIRLAM (high resolution limited area model)) predict an increase of 1 °C in temperature and a decrease of 5% in precipitation in 2020, and as a result, a decline in water resources from these areas of 10% according to AEMET (2009) (Spanish National Weather Service). These values could increase significantly according to other studies which calculate a decrease of 20% in precipitation and an increase of 3 °C in temperature in 2050 according to the PNACC (2011).

Taking into account these estimations, the main challenge facing the agricultural sector is making 70% more food available with fewer water resources. More water resources or a more efficient use of the resources will be needed to fit this objective. We cannot waste water resources and therefore we are forced to improve the current irrigation techniques and invest in automating the operation of irrigation canals, for instance, developing new control algorithms.

Many algorithms, in particular open loop algorithms, have difficulties and failures in implementation in real canals (Rogers & Goussard 1998). These difficulties are due to deviations between the predictive/control model and the reality.

The disturbances introduced into the canal lead to important deviations between the model and reality. The disturbances caused by infiltration can be controlled by a good coating; these disturbances are not usually important and they have constant little value. Due to the difficulty of distinguishing between small disturbances and small perturbations caused by conceptual errors of the model or water infiltration losses, they are all considered negligible in canal control. The disturbances introduced by changes in demand deliveries are more dangerous because the canal could overflow or could dry depending on the magnitude of the disturbance. These disturbances are usually unscheduled demand deliveries which are caused by flow rate extractions during the irrigation cycle.

In some cases, the farmers can adjust the supplies obtained by the irrigation community before the irrigation cycle begins; these changes in demand deliveries are known. Thus, these disturbances are known in advance and can be mitigated by an open loop controller (feedforward controller) such as GoRoSo (Bautista *et al.* 2003; Wahlin & Clemmens 2006; Soler *et al.* 2013).

The main problems in a canal are the disturbances caused by climatic variations (rainfalls and associated runoff) or unscheduled demands by farmers that extract more or less flow demanded by them, because these disturbances are more difficult to mitigate by a controller. In such cases, the CSE algorithm (Bonet 2015; Soler *et al.* 2015; Bonet *et al.* 2017) is an excellent tool to approach the real extractions in the canal in real time.

One way to protect the canal from these disturbances could be with ponds built by the irrigation community or the farmers themselves. A reservoir is able to store water according to the crop requirements and regulates the volume of water provided by the canal, but it is not always possible to build large reservoirs to regulate each pool of the canal. Another option to control the disturbances in canals, which operate in steady state, would be to use the wedge storage of every pool, but this only works with very low disturbances. In all other cases, a closed loop controller (feedback controller) is needed to modify the sluice gate trajectory to return to the desired canal state in real time, which is an easy way for our overall control diagram.

Beyond the disturbances, there are factors that significantly affect the canal, such as changes in Manning's coefficient which the on-line predictive controller needs to know to recalculate the optimal sluice gate trajectory. For this reason, we define a set of algorithms that operate all together in parallel but not all of them necessarily operate in real time. All the information obtained by these algorithms is used to supply the on-line predictive controller (GoRoSoBo) to meet its objective.

In this paper, our overall control diagram has been tested in several numerical examples. The control diagram was tested in a canal with a single pool; the geometry of the canal is based on Bautista's work (Baustista & Clemmens 1997). This kind of canal has been proposed by different authors such as Wylie (1969), Liu *et al.* (1992), Chevereau (1991), and Soler (2003) to test their control algorithms due to the ease of checking the results. In a second example, the overall control diagram was tested with the test cases (Clemmens *et al.* 1998) introduced by the ASCE Task Committee on Canal Automation Algorithms according to specific guidelines established to test canal controllers.

## DIAGRAM

Crop needs and desired hydrographs for canal outlets: The hydrographs at the lateral diversion points of the main canal are calculated on the basis of the water demands. They are fixed considering the farmers' requirements and other demands accepted by the watermaster. The behavior of the canal supplying these hydrographs determines the ‘desired behavior’ (

*Y**) at several cross sections.Off-line computation of the reference trajectories: The desired behavior (

*Y**) must be transmitted to the ‘reference trajectories calculation’ algorithm that determines the positions of each gate. This algorithm calculates the optimum behavior (*Y*(reference water level)), which is the one most similar to the desired behavior that is physically possible. We call ‘_{R}*U*’ the optimum gate trajectories calculated to obtain the optimum behavior (_{R}*Y*). They must be calculated off-line (e.g., with an anticipated irrigation cycle). There is an extensive bibliography of feedforward control algorithms._{R}Off-line parameter identification: This in case of variations of empirical parameters such as Manning's coefficient or gate discharge coefficients. The water level measurements obtained during a previous period of the irrigation cycle in steady state can be transmitted to the ‘off-line parameter identification’ algorithm for estimating these empirical parameters (Figure 2). The values of these parameters can be sent to the ‘off-line reference trajectories computation’ and the ‘on-line current state computation’ to be more accurate in their respective processes.

On-line current state computation: This module uses all values stored in the database such as ‘

*u*’ (gate position) and*y*(measured water level) which have been measured at the checkpoints in the canal during a past time horizon. If the water level measured at checkpoints is different from the desired water level, it is probably due to a disturbance introduced into the system or erroneous roughness coefficients used in the calculation. That is why it is very important to update Manning's coefficients periodically. The algorithm developed in this module is the CSE algorithm (Bonet_{M}*et al.*2017) which is able to estimate the disturbances that have modified the measured water level at checkpoints from the desired water levels. The algorithm can also obtain the current hydrodynamic state of the canal, which is very useful in any case. The hydrodynamic state of a canal is defined as the velocity and water level at each cross section of the canal.On-line predictive control: The reference trajectories plus the disturbances and the hydrodynamic canal state in real time are transmitted to the next algorithm, called ‘on-line predictive control’, which must react on-line in case of deviations between the observed data and the desired behavior at the checkpoints during each regulation period (e.g., every 5 min), due to unknown perturbations. The predictive control recalculates new gate positions (

*U*) to return to the reference behavior. There is an extensive bibliography of feedback control algorithms which recalculate the gate trajectories such as CARA (Marzouki 1989; Clemmens & Wahlin 2004; Wahlin & Clemmens 2006).

In this paper, we focus on the on-line computation task of our overall control diagram, because we are interested in improving the efficiency in water transport in ‘real time’.

The CSE solves an inverse problem implemented as a nonlinear optimization problem using the Levenberg–Marquardt method, in which the solution is the flow disturbances obtained from the water level variations measured at several points in the canal, usually next to the canal offtakes.

The GoRoSoBo solves an inverse problem (Equation (1)) implemented as a constrained nonlinear optimization problem using the Lagrange–Newton method, where the solution is the gate trajectories obtained from the output data partly from the CSE algorithm. In that sense, from the real extracted flow and the current hydrodynamic state, it is easy to obtain the water level disturbances at several points of the canal during a future horizon. Any deviation from the reference is fed back into the on-line predictive control so that this reduces the deviation of the controlled quantity from the reference. There are several canal control algorithms in the literature which use optimization methods to get the optimal gate positions to reach a particular water level target, for example, Sanders & Katopodes (1999).

*et al.*(2017). Instead, the development of GoRoSoBo algorithm is shown below: where Δ

*Y*represents the changes in water level at selected points of the canal, Δ

*U*represents a change in gate position, HIM’(

*U*) is the simplified hydraulic influence matrix that represents the influence of a gate position on the water level at different points of the canal, and HIM(

*U*) is the hydraulic influence matrix that represents the influence of a gate position on the water level and velocity along the canal.

The algorithm establishes a constrained nonlinear optimization problem to solve the last expression using the Lagrange–Newton method, as previously mentioned.

### The HIM matrix

The HIM matrix defines the influence of any gate opening over the hydraulic behavior of canal cross sections, usually limited to checkpoint sections. It is based on the full Saint-Venant equations, which describe the free surface flow in canals. In partial derivatives this system is of the hyperbolic, nonlinear, and second-order type. The two equations are based on the mass and momentum conservation.

*v*is the weighted average velocity of all the particles in a canal cross section,

*t*is the time,

*S*is the canal bottom slope,

_{0}*S*(

_{f}*y, v*) is the friction slope and

*c*is the celerity of a gravity wave, where

*A*(

*y*) is the area of the wetted surface or a cross section of the flow and

*T*(

*y*) is the top width of the free surface.

*S*=

_{fR}*Sf*(

*y*)

_{R},v_{R}*, S*=

_{fP}*Sf*(

*y*)

_{P},v_{P}*, S*=

_{fQ}*Sf*(

*y*)

_{Q},v_{Q}*,*and 0 ≤

*θ*≤ 1 is the coefficient of average time that indicates the type of numerical scheme used. When

*θ*= 1 the numerical scheme is implicit, if

*θ*= 0 is explicit, and when

*θ*= 1/2 the numerical scheme is in central differences or ‘method of the characteristic curves’.

If the flow conditions at points P and Q are known, *x _{P}, t_{P}, y_{P}, v_{P}* and

*x*are also known, so

_{Q}, t_{Q}, y_{Q}, v_{Q}*x*remain as unknown variables, which can be found by calculating the last four equations (Equation (3)) by using any of the methods solving non-linear equations, such as the Newton–Raphson method.

_{R}, t_{R}, y_{R}, v_{R}The way of calculating the influences shown in this section is closely linked to the numerical scheme of characteristic curves. However, usually this scheme is not exactly used because it gives the solution at a point R whose coordinates (*x _{R},t_{R}*) are unknown a priori. These coordinates are part of the solution and normally it is more important to know the solution of the flow conditions at specific points and at specific time instants. To solve this problem, first interpolate and then solve (Figure 3).

A structured grid like this one (Figure 2) creates a new nomenclature. Indeed, every variable will have a double index, where k refers to time and i to space. Thus, *y _{ik}* and

*v*represent the values for water level and average velocity at the co-ordinates

_{ik}*x*=

_{i}*i*Δ

*x*and

*t*=

_{k}*k*Δ

*t*where Δ

*x*and Δ

*t*are selected by the user.

*x*, and

_{P}, y_{R}, v_{R}*x*, where the values of

_{Q}*y*and

_{p}*y*are dependent on the value of

_{Q}*x*and

_{P}*x*evaluated using an interpolation function of second order too (to be coherent with the numerical procedure used) we have used the Lagrange factors (a way of representing quadratic splines). For a dummy variable

_{Q}*z*the result is: In this way the variables

*y*,

_{P}*v*,

_{P}*y*, and

_{Q}*v*become functions of

_{Q}*x*and

_{P}*x*, as follows:

_{Q}*et al.*(2017)). The interaction of this control structure with the flow can be described according to the mass and energy conservation equations (Equation (6)). where

*S*(

*y*) is the horizontal surface of the reception area in the checkpoint.

_{i}*A*(

*y*)

_{i}**v*is the incoming flow to checkpoint, defined in terms of water level and velocity.

_{i}*A*(

*y*)

_{o}**v*is the outflow from the checkpoint which continues along the canal, described in terms of water level and velocity.

_{o}*C*is the discharge coefficient of the sluice gate and

_{d}*a*is the sluice gate width. d is the checkpoint drop, and u is the gate opening.

_{c}*q*is the pumping offtake.

_{b}*q*(

_{s}*y*) is the outgoing lateral flow through the weir where

_{i}*C*is the discharge coefficient, as is the weir width and

_{s}*y*is the weir height measured from the bottom, called weir equation.

_{0}*Q*

_{offtake}(

*y*) is the outflow orifice flow where

_{i}*C*is the discharge coefficient,

_{0}*A*is the area of the offtake orifice, called orifice offtake equation. *The reverse flow is not taken into consideration for the GoRoSoBo algorithm.

_{0}*k*+ 1, that is, the incoming water level

*y*. In the same way, is defined as the existing water level at the first node of the downstream pool from the checkpoint at the same time

_{i}*k*+ 1, and

*y*the outgoing water level at the control structure. The same can be said for the velocities and where

_{o}On the other hand, , and *x _{Q}* are the unknown variables. In order to continue with the calculation of the influences of a general parameter, in our case, this parameter defines the gate position

*U*. Thus, applying the first derivative of the functions (Equation (7)) with respect to the gate position, we obtain the system of Equation (8).

*U*appears explicitly in the description. Despite the fact that the specific form of this function is still unknown, Equation (8) shows that the influence of the parameter

*U*on flow conditions at time

*k*+ 1 is the sum of the indirect influence of the conditions at instant k and the direct influence at instant

*k*+ 1 through the term ‘

*L*’, which represents the variation in the extraction flow. where:

As a summary, the method of characteristics is applied to the Saint-Venant equations in order to obtain algebraic equations to establish a relation between the influence parameter *U* and the hydrodynamic canal state; all the influences are lumped together in a global matrix, which is referred to as HIM(*U*). Based on this system of equations, and using the first derivative on an analytical process, the changes in flow behavior (water level and velocity) can be established due to a change of gate position at a point at a certain time instant.

### The optimization problem

The inverse problem (Equation (1)) is formulated as a constrained optimization problem. It is the classical non-linear problem with constraints and the method used to solve it is the Lagrange–Newton method.

*k*(Equation (9)) whose dimension is

_{F}*n*, where

_{y}*n*=

_{y}*k*

_{F}*×*

*n*is the final instant of the future time horizon and n

_{c}, k_{F}_{c}is the number of checkpoints. We define this vector as:

*x*(

*k*) which is defined as the vector containing the water level and velocity predictions established from the output data of the CSE algorithm at the time instant

*k*of all the discretization points: where

*y*(

_{i}*k*) and

*v*(

_{i}*k*) = water depth and mean velocity at point i; and

*n*= number of cross sections in which the canal is discretized. In this way, the vector

_{S}*x*(0) is the known initial condition.

The state vector at the current time defines the current hydrodynamic state; we show the state vector in a computational grid in Figure 5 (triangles).

*n*= (2 ×

_{x}*n*) ×

_{s}*k*, where ns is the number of cross sections of the canal: Because we are only interested in the water level at target points, we define a new vector called ‘prediction output vector’ that contains the water depth values given at a prescribed number of points (

_{F}*n*) from the time instant 1 to

_{c}*k*: The dimension of the prediction output vector is

_{F}*n*=

_{Y}*k*

_{F}*×*

*n*. The vector (Equation (12)) contains all water depth values exactly as the dots (x/t-dots) shown in Figure 5. If you look closely at this figure, the position of the elements of the vector (12) in the grid domain coincides with the elements of the desired water level vector.

_{c}*n*=

_{U}*n*

_{u}*×*

*K*is the number of gates, and

_{F}, nu*K*is the final operation period of the future horizon. Thus, the GoRoSoBo algorithm calculates the gate trajectories (Δ

_{F}*U*) to return from the predicted water level to the desired water level (Δ

*Y*) (1).

*et al.*1981; Fletcher 1987). In mathematical terms, the objective is to obtain the gate trajectories vector that minimizes the following performance criterion: where

*J*(

*U*) is the objective function, is the prediction vector from the regulation time step

*K*to

_{i}*K*is the desired water level vector,

_{F}, Y**Q*is a weighting matrix,

*C*is the discrete observer matrix (Malaterre 1994) and

*r*(U) is the ‘kth’ constraint function,

_{k}*I*(

*U*) is a set of equalities constraints, and

*NI*(

*U*) is a set of inequalities constraints.

*U*contains the gate trajectories (13).

## RESULTS AND DISCUSSION

### Numerical example: a canal introducing a single disturbance

The geometry of the canal proposed is based on Bautista *et al.* (1997). This canal was used by different authors including Wylie (1969), Liu *et al.* (1992), Chevereau (1991), and Soler (2003). The canal geometry adopted in our examples is based on Liu's example as well as the scheduled demand which was introduced in some scenarios.

The canal, with a trapezoidal section, is represented in Figure 7, and the general data are shown in Table 1. The characteristics of the checkpoint, sluice gate, pump station, and orifice offtake are shown in Tables 2 and 3.

Pool number | Pool length (km) | Bottom slope (%) | Side slopes (H:V) | Manning's coefficient (n) | Bottom width (m) | Canal depth (m) |
---|---|---|---|---|---|---|

I | 2.5 | 0.1 | 1.5:1 | 0.025 | 5 | 2.5 |

Pool number | Pool length (km) | Bottom slope (%) | Side slopes (H:V) | Manning's coefficient (n) | Bottom width (m) | Canal depth (m) |
---|---|---|---|---|---|---|

I | 2.5 | 0.1 | 1.5:1 | 0.025 | 5 | 2.5 |

Number of control structure or checkpoint | Gate discharge coefficient | Gate width (m) | Gate height (m) | Step (m) |
---|---|---|---|---|

0* | 0.61 | 5.0 | 2.0 | 0.0 |

Number of control structure or checkpoint | Gate discharge coefficient | Gate width (m) | Gate height (m) | Step (m) |
---|---|---|---|---|

0* | 0.61 | 5.0 | 2.0 | 0.0 |

*Reservoir gate.

Number of control structure or checkpoint | Discharge coef./diameter orifice offtake (m) | Orifice offtake height (m) | Lateral spillway height (m) | Lateral spillway width (m)/discharge coefficient |
---|---|---|---|---|

0* | – | – | – | – |

1 | 2/0.85 | 0.8 | 2.0 | 5/1.99 |

Number of control structure or checkpoint | Discharge coef./diameter orifice offtake (m) | Orifice offtake height (m) | Lateral spillway height (m) | Lateral spillway width (m)/discharge coefficient |
---|---|---|---|---|

0* | – | – | – | – |

1 | 2/0.85 | 0.8 | 2.0 | 5/1.99 |

*Reservoir gate.

In these examples an upstream large reservoir is considered, whose water level *H*_{reservoir} is constantly 3 m throughout the test. At the end of the last pool there is a control structure with an orifice offtake and a pump station. The flow through the orifice depends on the level over the orifice and the disturbance is introduced by the pump station. This is the downstream boundary condition. This example starts from an initial steady state (Tables 4 and 5) with a specific demand delivery constant at the end of the pool (5 m^{3}/s through the orifice offtake), and the disturbance is not introduced initially.

Control structure | Initial flow rate (m^{3}/s) | Control structure | Initial water level (m) |
---|---|---|---|

Gate 1 | 5.0 | Checkpoint 1 | 1.6 |

Control structure | Initial flow rate (m^{3}/s) | Control structure | Initial water level (m) |
---|---|---|---|

Gate 1 | 5.0 | Checkpoint 1 | 1.6 |

Control structure | Flow delivered by an orifice offtake (m^{3}/s) |
---|---|

Gravity outlet 1 | 5.0 |

Control structure | Flow delivered by an orifice offtake (m^{3}/s) |
---|---|

Gravity outlet 1 | 5.0 |

### The disturbance

^{3}/s. This disturbance is introduced to the computer model as a pump flow change and we obtain the water level measurements from this mode. Once these water level measurements are introduced in the CSE algorithm, it will propose the pump flow trajectories that describe with best accuracy the variation of water level at the checkpoints during the past time horizon and the hydrodynamic state at the current time step. All these data are sent to GoRoSoBo which sets the predicted water level vector for a future horizon and recalculates the optimum gate trajectories to ensure the target level over the orifices for this future horizon. As we showed in Figure 2, all these actions are repeated every regulation period.

The constraints of the optimization problem are imposed on the gate position, so that the sluice gate opening must not be greater or smaller than *U*_{max} or *U*_{min}, respectively, and the gate movements between successive regulation periods (*dU*_{max}) should be physically acceptable (Table 6).

U (%) _{min} | U (%) _{max} | dU (%) _{max} |
---|---|---|

0.05 | 90 | 2.5 |

U (%) _{min} | U (%) _{max} | dU (%) _{max} |
---|---|---|

0.05 | 90 | 2.5 |

### Results

*et al.*2013) and with a continuous line the results obtained by the overall control diagram at the target point in Figure 8. We can compare the sluice gate trajectories obtained by GoRoSoBo vs. the sluice gate trajectories obtained by GoRoSo in Figure 9.

At the beginning of the test, the flow in the canal is steady and the scheduled delivery is constant for all demand periods, and the sluice gate position remains fixed. After the first 30 minutes (1,800 s), a disturbance is introduced into the system, although the algorithm has no notice until the next regulation period (2,100 s) once the water level is measured at the checkpoint, and the water level has already been reduced by more than 10 cm, from 1.6 cm to 1.50 cm (Figure 8).

Once the water level measurements are sent to CSE, it calculates the extracted flow vector for a past time horizon. CSE establishes the disturbance introduced in the system. This information is sent to GoRoSoBo which modifies the sluice gate trajectory to keep the water level constant at the checkpoint which is not going to increase until three regulation periods later (2,700 s). This is because once the sluice gate generates a wave for increasing the water level at the checkpoint, the wave has to reach the checkpoint (time delay). Once the wave arrives at the checkpoint, the water level increases quickly, recovering the target level of 1.6 m at 3,700 s (Figure 8).

The maximum deviation between the water level measured and desired is around 27.5 cm, from 1.6 cm to 1.325 cm, so the sluice gate movement must be quite significant in order to quickly reduce this deviation (Figure 9). The gate position changes from 0.125 to 0.34 m during a regulation period (1,800–2,100 s).

The feedforward algorithm (GoRoSo) does not modify the gate trajectory during the irrigation cycle because the scheduled demand is constant and the algorithm does not know that someone has introduced a disturbance in the canal and the water level decreases at checkpoint. The water level decreases to 1.1 m so the flow extracted by the offtake is close to 3 m^{3}/s (Figure 8). In that sense, the total extracted flow is 5 m^{3}/s, but 3 m^{3}/s are extracted by the offtake and 2 m^{3}/s are extracted by someone who introduced the disturbance. The difference in gate opening between the feedforward and the feedback algorithms at 9,000 s is due to the disturbance. A difference of 0.06 m in gate opening represents a flow of 2 m^{3}/s in this canal with this flow condition.

### Numerical example: a canal with multiple disturbances at the same time: ASCE test cases

In this numerical example, we want to evaluate our overall control diagram in a canal with several pools with multiple flow extractions at the same time. In that case, we introduce the test cases proposed by the ASCE committee to evaluate control algorithms.

*et al.*1998) in several scenarios. Each canal has eight pools separated by undershot sluice gates. All the canal pools have been discretized and numbered in the direction of flow from upstream to downstream. Geometric characteristics of canal 1 (Maricopa Stanfield canal) are shown in Tables 7 and 8 and Figure 10, and for canal 2 (Corning canal) in Tables 9 and 10 and Figure 11. In both canals, there are gravity outlet orifices at the downstream end of each canal pool but only in canal 2 is there direct pumping at the end of the last pool. The ASCE Committee proposes four cases to test feedback controllers in real time, two in the Corning canal and two in the Maricopa Stanfield. We only test one case in each canal with our overall control diagram taking into account the case with more unscheduled flow changes. The initial conditions for all tests are shown in Table 11 (Maricopa Standfield) and Table 12 (Corning canal). All scenarios are evaluated in Bonet (2015).

Pool number | Pool length (km) | Bottom slope | Side slopes (H:V) | Manning's coefficient (n) | Bottom width (m) | Canal depth (m) |
---|---|---|---|---|---|---|

I | 0.1 | 2 × 10^{−3} | 1.5:1 | 0.014 | 1 | 1.1 |

II | 1.2 | 2 × 10^{−3} | 1.5:1 | 0.014 | 1 | 1.1 |

III | 0.4 | 2 × 10^{−3} | 1.5:1 | 0.014 | 1 | 1.0 |

IV | 0.8 | 2 × 10^{−3} | 1.5:1 | 0.014 | 0.8 | 1.1 |

V | 2 | 2 × 10^{−3} | 1.5:1 | 0.014 | 0.8 | 1.1 |

VI | 1.7 | 2 × 10^{−3} | 1.5:1 | 0.014 | 0.8 | 1.0 |

VII | 1.6 | 2 × 10^{−3} | 1.5:1 | 0.014 | 0.6 | 1.0 |

VIII | 1.7 | 2 × 10^{−3} | 1.5:1 | 0.014 | 0.6 | 1.0 |

Pool number | Pool length (km) | Bottom slope | Side slopes (H:V) | Manning's coefficient (n) | Bottom width (m) | Canal depth (m) |
---|---|---|---|---|---|---|

I | 0.1 | 2 × 10^{−3} | 1.5:1 | 0.014 | 1 | 1.1 |

II | 1.2 | 2 × 10^{−3} | 1.5:1 | 0.014 | 1 | 1.1 |

III | 0.4 | 2 × 10^{−3} | 1.5:1 | 0.014 | 1 | 1.0 |

IV | 0.8 | 2 × 10^{−3} | 1.5:1 | 0.014 | 0.8 | 1.1 |

V | 2 | 2 × 10^{−3} | 1.5:1 | 0.014 | 0.8 | 1.1 |

VI | 1.7 | 2 × 10^{−3} | 1.5:1 | 0.014 | 0.8 | 1.0 |

VII | 1.6 | 2 × 10^{−3} | 1.5:1 | 0.014 | 0.6 | 1.0 |

VIII | 1.7 | 2 × 10^{−3} | 1.5:1 | 0.014 | 0.6 | 1.0 |

Target points | Gate discharge coefficient | Gate width (m) | Gate height (m) | Step (m) | Length from gate 1 (km) | Orifice offtake height (m) |
---|---|---|---|---|---|---|

0 | 0.61 | 1.5 | 1.0 | 1.0 | 0 | – |

1 | 0.61 | 1.5 | 1.1 | 1.0 | 0.1 | 0.45 |

2 | 0.61 | 1.5 | 1.1 | 1.0 | 1.3 | 0.45 |

3 | 0.61 | 1.5 | 1.0 | 1.0 | 1.7 | 0.40 |

4 | 0.61 | 1.2 | 1.1 | 1.0 | 2.5 | 0.45 |

5 | 0.61 | 1.2 | 1.1 | 1.0 | 4.5 | 0.45 |

6 | 0.61 | 1.2 | 1.0 | 1.0 | 6.2 | 0.40 |

7 | 0.61 | 1.0 | 1.0 | 1.0 | 7.8 | 0.40 |

8 | – | – | – | – | 9.5 | 0.40 |

Target points | Gate discharge coefficient | Gate width (m) | Gate height (m) | Step (m) | Length from gate 1 (km) | Orifice offtake height (m) |
---|---|---|---|---|---|---|

0 | 0.61 | 1.5 | 1.0 | 1.0 | 0 | – |

1 | 0.61 | 1.5 | 1.1 | 1.0 | 0.1 | 0.45 |

2 | 0.61 | 1.5 | 1.1 | 1.0 | 1.3 | 0.45 |

3 | 0.61 | 1.5 | 1.0 | 1.0 | 1.7 | 0.40 |

4 | 0.61 | 1.2 | 1.1 | 1.0 | 2.5 | 0.45 |

5 | 0.61 | 1.2 | 1.1 | 1.0 | 4.5 | 0.45 |

6 | 0.61 | 1.2 | 1.0 | 1.0 | 6.2 | 0.40 |

7 | 0.61 | 1.0 | 1.0 | 1.0 | 7.8 | 0.40 |

8 | – | – | – | – | 9.5 | 0.40 |

Pool number | Pool length (km) | Bottom slope | Side slopes (H:V) | Manning's coefficient (n) | Bottom width (m) | Canal depth (m) |
---|---|---|---|---|---|---|

I | 7 | 1 × 10^{−4} | 1.5:1 | 0.02 | 7 | 2.5 |

II | 3 | 1 × 10^{−4} | 1.5:1 | 0.02 | 7 | 2.5 |

III | 3 | 1 × 10^{−4} | 1.5:1 | 0.02 | 7 | 2.5 |

IV | 4 | 1 × 10^{−4} | 1.5:1 | 0.02 | 6 | 2.3 |

V | 4 | 1 × 10^{−4} | 1.5:1 | 0.02 | 6 | 2.3 |

VI | 3 | 1 × 10^{−4} | 1.5:1 | 0.02 | 5 | 2.3 |

VII | 2 | 1 × 10^{−4} | 1.5:1 | 0.02 | 5 | 1.9 |

VIII | 2 | 1 × 10^{−4} | 1.5:1 | 0.02 | 5 | 1.9 |

Pool number | Pool length (km) | Bottom slope | Side slopes (H:V) | Manning's coefficient (n) | Bottom width (m) | Canal depth (m) |
---|---|---|---|---|---|---|

I | 7 | 1 × 10^{−4} | 1.5:1 | 0.02 | 7 | 2.5 |

II | 3 | 1 × 10^{−4} | 1.5:1 | 0.02 | 7 | 2.5 |

III | 3 | 1 × 10^{−4} | 1.5:1 | 0.02 | 7 | 2.5 |

IV | 4 | 1 × 10^{−4} | 1.5:1 | 0.02 | 6 | 2.3 |

V | 4 | 1 × 10^{−4} | 1.5:1 | 0.02 | 6 | 2.3 |

VI | 3 | 1 × 10^{−4} | 1.5:1 | 0.02 | 5 | 2.3 |

VII | 2 | 1 × 10^{−4} | 1.5:1 | 0.02 | 5 | 1.9 |

VIII | 2 | 1 × 10^{−4} | 1.5:1 | 0.02 | 5 | 1.9 |

Target points | Gate width (m) | Gate height (m) | Step (m) | Length from gate 1 (Km) | Orifice offtake height (m) | Lateral spillway height (m) |
---|---|---|---|---|---|---|

0 | 7 | 2.3 | 0.2 | 0 | – | 3 |

1 | 7 | 2.3 | 0.2 | 7 | 1.05 | 2.5 |

2 | 7 | 2.3 | 0.2 | 10 | 1.05 | 2.5 |

3 | 7 | 2.3 | 0.2 | 13 | 1.05 | 2.5 |

4 | 6 | 2.1 | 0.2 | 17 | 0.95 | 2.3 |

5 | 6 | 2.1 | 0.2 | 21 | 0.95 | 2.3 |

6 | 5 | 1.8 | 0.2 | 24 | 0.85 | 1.9 |

7 | 5 | 1.8 | 0.2 | 26 | 0.85 | 1.9 |

8 | – | – | – | 28 | 0.85 | 1.9 |

Target points | Gate width (m) | Gate height (m) | Step (m) | Length from gate 1 (Km) | Orifice offtake height (m) | Lateral spillway height (m) |
---|---|---|---|---|---|---|

0 | 7 | 2.3 | 0.2 | 0 | – | 3 |

1 | 7 | 2.3 | 0.2 | 7 | 1.05 | 2.5 |

2 | 7 | 2.3 | 0.2 | 10 | 1.05 | 2.5 |

3 | 7 | 2.3 | 0.2 | 13 | 1.05 | 2.5 |

4 | 6 | 2.1 | 0.2 | 17 | 0.95 | 2.3 |

5 | 6 | 2.1 | 0.2 | 21 | 0.95 | 2.3 |

6 | 5 | 1.8 | 0.2 | 24 | 0.85 | 1.9 |

7 | 5 | 1.8 | 0.2 | 26 | 0.85 | 1.9 |

8 | – | – | – | 28 | 0.85 | 1.9 |

Pool number | Offtake initial flow (m^{3}/s) | Check initial flow (m^{3}/s) | Unscheduled offtake changes at 2 hours (m^{3}/s) | Check final flow (m^{3}/s) |
---|---|---|---|---|

Heading | – | 2.0 | – | 2.0 |

1 | 0.2 | 1.8 | – | 1.8 |

2 | 0.0 | 1.8 | 0.2 | 1.6 |

3 | 0.4 | 1.4 | −0.2 | 1.4 |

4 | 0.0 | 1.4 | 0.2 | 1.2 |

5 | 0.0 | 1.4 | 0.2 | 1.0 |

6 | 0.3 | 1.1 | −0.1 | 0.8 |

7 | 0.2 | 0.9 | – | 0.6 |

8 | 0.9 | 0.0 | −0.3 | 0.0 |

Pool number | Offtake initial flow (m^{3}/s) | Check initial flow (m^{3}/s) | Unscheduled offtake changes at 2 hours (m^{3}/s) | Check final flow (m^{3}/s) |
---|---|---|---|---|

Heading | – | 2.0 | – | 2.0 |

1 | 0.2 | 1.8 | – | 1.8 |

2 | 0.0 | 1.8 | 0.2 | 1.6 |

3 | 0.4 | 1.4 | −0.2 | 1.4 |

4 | 0.0 | 1.4 | 0.2 | 1.2 |

5 | 0.0 | 1.4 | 0.2 | 1.0 |

6 | 0.3 | 1.1 | −0.1 | 0.8 |

7 | 0.2 | 0.9 | – | 0.6 |

8 | 0.9 | 0.0 | −0.3 | 0.0 |

Pool number | Offtake initial flow (m^{3}/s) | Check initial flow (m^{3}/s) | Unscheduled offtake changes at 2 hours (m^{3}/s) | Resulting check flow (m^{3}/s) |
---|---|---|---|---|

Heading | – | 13.7 | – | 2.7 |

1 | 1.7 | 12.0 | −1.5 | 2.5 |

2 | 1.8 | 10.2 | −1.5 | 2.2 |

3 | 2.7 | 7.5 | −2.5 | 2.0 |

4 | 0.3 | 7.2 | – | 1.7 |

5 | 0.2 | 7.0 | – | 1.5 |

6 | 0.8 | 6.2 | −0.5 | 1.2 |

7 | 1.2 | 5.0 | −1.0 | 1.0 |

8 | 0.3 + 2.0* | 2.7 | −2.0* | 0.7 |

Pool number | Offtake initial flow (m^{3}/s) | Check initial flow (m^{3}/s) | Unscheduled offtake changes at 2 hours (m^{3}/s) | Resulting check flow (m^{3}/s) |
---|---|---|---|---|

Heading | – | 13.7 | – | 2.7 |

1 | 1.7 | 12.0 | −1.5 | 2.5 |

2 | 1.8 | 10.2 | −1.5 | 2.2 |

3 | 2.7 | 7.5 | −2.5 | 2.0 |

4 | 0.3 | 7.2 | – | 1.7 |

5 | 0.2 | 7.0 | – | 1.5 |

6 | 0.8 | 6.2 | −0.5 | 1.2 |

7 | 1.2 | 5.0 | −1.0 | 1.0 |

8 | 0.3 + 2.0* | 2.7 | −2.0* | 0.7 |

*Changes in downstream pump discharge with no change in offtake flow.

The gate discharge coefficient for the gates in Corning canal is equal to the Maricopa Stanfield.

#### Scenario: Test case 1-2 (Maricopa Stanfield)

The unscheduled deliveries are more significant at target 8 in test case 1-2, where the flow rate changes from 0.9 m^{3}/s to 0.6 m^{3}/s (33%).

#### Scenario: Test case 2-2 (Corning canal)

The unscheduled water deliveries in test case 2-2 are very relevant, as these water deliveries are a huge amount of flow and these unscheduled water deliveries represent a high percentage with respect to the total flow of the canal. In this case, the unscheduled water deliveries are 81% on the total flow and these water changes are more significant at target 8, where the flow rate changes from 5 m^{3}/s to 1 m^{3}/s. For this reason, it should be noted that 15 minutes after introducing the unscheduled water deliveries, the water level at target 8 is 1.93 meters, so the lateral spillway must operate. This is the maximum flow rate change at a target in all test cases.

### Constraints

The constraints are determined as a certain percentage of the gate height, imposed at the gate positions (Table 13).

U (%) _{min} | U (%) _{max} | dU (%) _{max} | |
---|---|---|---|

Maricopa Stanfield | 2 | 90 | 5 |

Corning canal | 0.5 | 90 | 5 |

U (%) _{min} | U (%) _{max} | dU (%) _{max} | |
---|---|---|---|

Maricopa Stanfield | 2 | 90 | 5 |

Corning canal | 0.5 | 90 | 5 |

### Results

We show the results obtained in each test case dividing these into water level at checkpoints and gate trajectories which were also divided in a graph with only four checkpoints or four gates, respectively, to analyze the results properly.

#### Test case 1-2 (Maricopa Stanfield)

^{3}/s at checkpoint 8 is introduced, with a total flow rate change at target 8 from 0.6 m

^{3}/s to 0.9 m

^{3}/s, that is, a flow rate change of around 50% from the initial one.

Taking into account that all water delivery changes are relevant, such as these which are 0.2 m^{3}/s at targets 1 to 5, 0.1 m^{3}/s at target 6, and 0.3 m^{3}/s at target 8, in a canal that has no storage capacity and a steep slope, this is a hard test for a control algorithm.

We can check that the biggest change in gate position at 18,000 s is gate 8 which is logical because the more important flow change is at checkpoint 8.

The sluice gate trajectories have the same movement as the water level at the checkpoints, due to the gate trajectories and water levels fluctuating around the desired solution.

The water level decreases in all targets at 9,900 s due to changes in gate trajectories, and the water level returns to the desired water level in all checkpoints at 21,600 s.

The water levels were almost equal to the desired levels at the end of the test, with a maximum error between them of around 0.3 cm.

#### Test case 2-2 (Corning canal)

Redirecting the water level at checkpoints to the desired water level in a centralized system is not the duty of only one sluice gate. It is possible that some sluice gates have increased the water level error in some checkpoints at a regulation period (for instance, checkpoints 4 and 5) when there are no unscheduled offtake changes in these checkpoints, but the main objective is to redirect the measured water level to the desired values in all checkpoints as soon as possible.

The water level decreases at all checkpoints at 25,200 s (7 hours) and the water level returns to the desired value at all the checkpoints at 30,600 s. The unscheduled water delivery is quite important at target 8, introducing a total flow change of 4 m^{3}/s at 7,200 s because the total flow rate at pool 8 is 5 m^{3}/s at the initial time step and the flow rate at pool 8 is 1 m^{3}/s at 7,200 s, so the flow change is close to 80% of the total flow rate of the pool.

The gate movements are almost zero in the last 2 hours of the test because the water levels at checkpoints are equal to the desired values at the end of the test, with a maximum error of around 0.5 cm.

### Performance indicators

These indicators were introduced to compare different control algorithms, as it is not easy to evaluate them and judge the controller's ability to deliver water. This is the reason why the ASCE Task Committee devised these test cases.

Although not all these performance indicators evaluate the deviations between the measured and the desired water level at the checkpoints, they take into account other variables such as the changes in flow.

We compare the performance indicators obtained with GoRoSoBo in the test cases with the same indicators obtained with other control algorithms like CLIS (Liu *et al.* 1998) and PILOTE (Malaterre 1995) (see Tables 14 and 15). CLIS is based on an inverse solution method of the Saint-Venant equations, and it is designed for the automation of demand-oriented systems. PILOTE is a LQR closed-loop controller and it is obtained from the steady-state solution of the Riccati equation, a Kalman filter is used to reconstruct the state variables and the unknown perturbations from a reduced number of observed variables.

Test case 1-2 tuned unscheduled | ||||||||
---|---|---|---|---|---|---|---|---|

MAE (%) | IAE (%) | StE (%) | IAQ (m^{3}/s) | |||||

12–24 h | 12–24 h | 12–24 h | 12–24 h | |||||

Max | Average | Max | Average | Max | Average | Max | Average | |

CLIS | 34.5 | 14.2 | 5.0 | 2.0 | 3.6 | 1.1 | 0.2 | 0.1 |

PILOTE | 43.0 | 24.9 | 9.2 | 5.2 | 11.2 | 2.9 | 2.9 | 1.4 |

GoRoSoBo | 33.5 | 10.3 | 5.0 | 1.2 | 2.1 | 0.7 | 6.7 | 3.6 |

Test case 1-2 tuned unscheduled | ||||||||
---|---|---|---|---|---|---|---|---|

MAE (%) | IAE (%) | StE (%) | IAQ (m^{3}/s) | |||||

12–24 h | 12–24 h | 12–24 h | 12–24 h | |||||

Max | Average | Max | Average | Max | Average | Max | Average | |

CLIS | 34.5 | 14.2 | 5.0 | 2.0 | 3.6 | 1.1 | 0.2 | 0.1 |

PILOTE | 43.0 | 24.9 | 9.2 | 5.2 | 11.2 | 2.9 | 2.9 | 1.4 |

GoRoSoBo | 33.5 | 10.3 | 5.0 | 1.2 | 2.1 | 0.7 | 6.7 | 3.6 |

Test case 2-2 tuned unscheduled | ||||||||
---|---|---|---|---|---|---|---|---|

MAE (%) | IAE (%) | StE (%) | IAQ (m^{3}/s) | |||||

12–24 h | 12–24 h | 12–24 h | 12–24 h | |||||

Max | Average | Max | Average | Max | Average | Max | Average | |

CLIS | 21.1 | 14.9 | 7.6 | 2.8 | 0.7 | 0.4 | 9.7 | 5.5 |

PILOTE | 34.2 | 17.1 | 10.6 | 7.1 | 8.8 | 4.3 | 10.4 | 6.1 |

GoRoSoBo | 13.6 | 7.8 | 6.3 | 2.1 | 0.8 | 0.5 | 15.2 | 11.7 |

Test case 2-2 tuned unscheduled | ||||||||
---|---|---|---|---|---|---|---|---|

MAE (%) | IAE (%) | StE (%) | IAQ (m^{3}/s) | |||||

12–24 h | 12–24 h | 12–24 h | 12–24 h | |||||

Max | Average | Max | Average | Max | Average | Max | Average | |

CLIS | 21.1 | 14.9 | 7.6 | 2.8 | 0.7 | 0.4 | 9.7 | 5.5 |

PILOTE | 34.2 | 17.1 | 10.6 | 7.1 | 8.8 | 4.3 | 10.4 | 6.1 |

GoRoSoBo | 13.6 | 7.8 | 6.3 | 2.1 | 0.8 | 0.5 | 15.2 | 11.7 |

The *MAE* is the maximum deviation of the controlled water level at the checkpoint with regard to the desired water level. From the values obtained by GoRoSoBo with this performance indicator, GoRoSoBo shows the best results in all tests.

The *IAE* is the integrated deviation of controlled water level on the target water depth. From the maximum and average *IAE* values obtained by GoRoSoBo, this algorithm shows the best results in three of four cases.

The *StE* is the deviation of controlled water level at steady state over target water depth, so we only consider the last 2 hours of the irrigation cycle to calculate this performance indicator. GoRoSoBo obtained the best results in three cases and it was the second in two cases.

The *IAQ* gives us an idea of how many changes in flow rate are introduced by the gates to reach the desired flow rate at the end of the test. GoRoSoBo shows the highest values on this performance indicator. When a control algorithm has, as an only objective, to maintain a desired water level at the checkpoints, the gate trajectories are calculated for that objective and some performance indicators benefit more (i.e., MAE/IAE) than others (i.e., IAQ). In that sense, when significant flow changes are introduced in the system, quick changes on gate trajectories are necessary to recover the desired value at checkpoints as soon as possible. For this reason, increasing the constraints more in gate movements to obtain better values of this performance indicator is not our principal objective. We can conclude that the magnitude of flow changes through the gates is acceptable and the control algorithm maintains the measured water level close to the desired water level and there exists a functional constraint to restrict the gate movement.

## CONCLUSIONS

The GoRoSoBo algorithm is able to find the optimum gate trajectory during a predictive horizon taking into account demand deliveries, initial gate trajectories, desired water level vector, disturbances, and the current canal state which are obtained by CSE. All these data are introduced into GoRoSoBo and the algorithm recalculates the optimum gate trajectory to keep the desired water level constant at checkpoints.

The GoRoSoBo algorithm uses the Lagrange–Newton method to solve a constrained optimization problem. This method is considered the most efficient when the Jacobian matrix and the Hessian matrix, which are used in the computation of the gate trajectories, have been compiled.

The introduction of constraints was absolutely necessary to ensure stability in our optimized problem, due to inherent instability in the unconstrained problem, which is caused by the condition number of the Hessian and the HIM matrixes. Not all elements of the HIM have similar values; there are gates that have a significant influence in certain checkpoints and little influence in others. This disparity of influence between elements of the matrix inevitably leads to a band matrix, of course, badly conditioned. This was a reason to use the Marquardt coefficient which improves the Hessian matrix condition.

The watermaster can be more or less strict about the constraints in the optimization problem, because he must consider the main priority of the canal. For instance, if the main priority is to maintain the water level constant at any price, the constraints on the gates would be lax, so the constraints will not have an important role in computing the gate trajectories by GoRoSoBo. If the main priority is to keep the water level constant, with reduced gate movements, the constraints would be more restrictive, so they will have an important role in the optimization problem.

The *IAQ* value is much higher in GoRoSoBo than in other controllers proposed. This index is strongly linked to the main priority of the canal. In the case where the priority is to keep the water level constant at the checkpoints, the gate movements are significant and the flow rate variations through the sluice gate also, so the *IAQ* value will also be significant.

One of the main problems of our overall control diagram is the computation time. In cases where the predictive horizon is large, the regulation interval is small, the canal length and the number of checkpoints are significant, the calculation time of the algorithm is too large for operating in real time and much more so if it is necessary to update the hydraulic influence matrix in every regulation time step. A solution with the time conflict would be to parallelize the algorithm with OpenMP to reduce the CPU time, as current workstations have several processors with dozens of cores which would increase the speed-up in solving the problem. Another possibility to speed-up the calculating process would be if GoRoSoBo calculated the gate trajectories using the HIM in a previous time step, while at the same time, another algorithm calculated the HIM at the current time step to be used by GoRoSoBo in the next time step.

The results obtained by GoRoSoBo using our overall control diagram have been very satisfactory, as we have shown in the case of a simple canal and the test cases by the ASCE.