GoRoSoBo: an overall control diagram to improve the efficiency of water transport systems in real time

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³ in 2008 to 3,858 km³ 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.

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).

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

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.

If the flow conditions at points P and Q are known, x_P, t_P, y_P, v_P and x_Q, t_Q, y_Q, v_Q are also known, so x_R, t_R, y_R, v_R 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.

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_ik represent the values for water level and average velocity at the co-ordinates x_i = iΔx and t_k = kΔt where Δx and Δt are selected by the user.

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).

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 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.

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).

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.

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³/s through the orifice offtake), and the disturbance is not introduced initially.

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).

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³/s (Figure 8). In that sense, the total extracted flow is 5 m³/s, but 3 m³/s are extracted by the offtake and 2 m³/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³/s in this canal with this flow condition.

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.

The gate discharge coefficient for the gates in Corning canal is equal to the 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³/s to 0.6 m³/s (33%).

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³/s to 1 m³/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.

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

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.

Taking into account that all water delivery changes are relevant, such as these which are 0.2 m³/s at targets 1 to 5, 0.1 m³/s at target 6, and 0.3 m³/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.

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³/s at 7,200 s because the total flow rate at pool 8 is 5 m³/s at the initial time step and the flow rate at pool 8 is 1 m³/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.

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.

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.

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.