A Model Predictive Water-Level Difference Control Method for Automatic Control of Irrigation Canals

Abstract

In this paper, automatic control of the water level in an irrigation canal by automatic regulation of intermediate gates was studied. Previous scholars have proposed a water level difference control strategy that works to keep relative deviations in all pools the same for a particular situation where the operator does not have full control over the canal inflow, with the centralized linear quadratic regulator (LQR) control method used. While in practice, the deviation tolerance of pools may differ in some canals which limits the applicability of the control strategy. In this work, a weight coefficient was added to the deviation and the algorithm was improved to keep the relative deviations to certain proportions. The model predictive control (MPC) method was then used with this improved control strategy and was compared to the LQR control method using the same control strategy. The results showed that the improved strategy can keep the water level deviations in all pools to certain proportions, as is our objective. Also, under this difference control strategy, the MPC method greatly improved the control performance compared to the LQR control method.

1. Introduction

Irrigation systems are built to deliver large amounts of water from a place with sufficient water to a place where water is a scarce resource. However, a considerable amount of water is wasted due to evaporation, leakage, and lack of control [1]. Innovative and adaptive modernization of irrigation systems is a key factor in improving water use efficiencies.

Improving the operations of irrigation water delivery systems has been an important topic for several decades [2]. Canal automation has evolved to the point where most new canal designs and canal modernization projects include some level of automation [3,4,5]. Different automatic control methods have been designed, implemented, and developed for canal operation. These automatic control methods can be roughly divided into single-input, single-output (SISO) and multiple-input, and multiple-output (MIMO) controllers [6]. In the SISO control method, a single check gate is controlled according to a single water level input, such as in the proportional-integral (PI) control method [7] and in improved forms of the PI control method [8,9]. In the MIMO control method, all check gates are simultaneously controlled according to water level inputs at all monitoring points, such as in the linear quadratic regulator (LQR) [10] and model predictive control (MPC) methods [11].

Although the form of these controllers differs, they all aim to maintain the downstream water level of canal pools at a certain level. Meanwhile, high performance on-field practices require that the water be delivered with sufficient reliability, equity, and flexibility [12]. In some conditions, some control structures may not be fully controlled or cannot be controlled flexibly by an automated system, so these control structures separate the whole canal system into canal segments, with the outflow gate and inflow check gate of each segment fixed. Clemmens et al. [13] proposed another operational strategy to control water levels in long main canals with considerable transmission time or when there is no control on canal inflow and outflow. In this proposal, the water level differences between adjacent pools are the controller inputs rather than the traditional pool water level deviations. The goal of this control method is to make the water levels in all pools change at the same rate so that the main canal consequently behaves as a storage reservoir when canal inflow and outflow do not match. Guan et al. [14] applied this control strategy with a centralized LQR to a model of the Central Arizona Project (CAP) main canal, with results showing that the method is a promising way to accommodate mismatches in supply and demand through in-canal storage when there is no full control over the canal inflow. Hashemy et al. [15] used this kind of control system with a simple PI control plus filter (PIF) method to deal with delivery disturbance with a fixed head gate. Kong et al. [16] applied this difference control strategy with a PI control method to a condition in which the inflow is to be changed significantly.

As the system controls the deviation differences in adjacent pools at the downstream side of the pools, the trend of the control result is that all pool deviations tend to be the same. While under some conditions, for example, if the offtake of the pool delivered significantly more water to more important users than to other pools, the water level disturbance of the pool should be less than that of other pools and the allowable water level deviations of each canal pool could be different. This requires the water levels to change at different rates, however not chaotically. In this paper, the water level difference control strategy is improved to meet this demand. Weight coefficients are added to water level deviations and the control target is not the actual water level difference between adjacent pools, however it is the difference between water level deviations amplified via the weight coefficient. Also, as a centralized control method will be used here, it is not just a simple input change—the related parameters also change. A model predictive control (MPC) method can in general take into account future water delivery disturbances and take action before those disturbances occur when there is a scheduled delivery change, which is in contrast to the LQR method that is used in traditional water deviation control [17]. Therefore, the MPC control is used in the simulation discussed in this paper together with the water level difference control. The result of this method is also later compared to LQR method results.

2. Materials and Methods

2.1. Test Canal and Scenarios

To demonstrate the performance of the proposed water level difference controller, a simulation model of the last six pools of the Middle Route Project (MRP) for the South-to-North Water Transfer Project was built for the simulation study. The MRP is the largest water transfer project in China and it delivers water from Danjiangkou Reservoir in south central China to water-scarce areas in north China, such as Beijing and Tianjin. The main canal of the MRP is 1273 km long and its design flow rate is 350 m³/s at the upstream side and 50 m³/s at the downstream side. There are 63 check structures and 88 offtakes and there is no online reservoir. The studied canal system of the MRP is located at the junction of Hebei Province and Beijing City. At the end of the canal system, water is delivered to Huinanzhuang Pumping Station in Beijing. The total length of this canal system is about 112 km. The six pools together were treated as an independent canal in this simulation, with a constant upstream water level boundary and an uncontrolled head gate, which is not the actual condition, however is suitable for the application of the proposed control method. The initial flow of the head gate was 94.5 m³/s and the outflow of the last gate was 35 m³/s. The layout is shown in Figure 1. The six pools are denoted as Pool 1, Pool 2, Pool 3, Pool 4, Pool 5, and Pool 6. The offtake flow in pool i is denoted as d(i), the inflow gate of pool i is denoted as Gate i-1, and the outflow gate is Gate i. The water level in pool i refers to the water level at the downstream end of pool i, which is also the water level immediately upstream from Gate i as it is often the water level at the downstream end that is of concern, where the pool is at its maximum water depth.

As the last pool delivers much more flow to a larger city, the disturbance of the last pool should be less than that of other pools. Thus, the weight coefficient of the last pool should be larger than that of the other pools. The focus of this paper is to obtain a proper distribution of flow imbalance as people want this when there is limited water supply at the head source. Therefore, different water level difference demand conditions were set. In test scenario No. 1, the offtake flow in Pool 6 increased by 5 m³/s, scheduled at time 10 h; the water level deviation tolerance of all pools was the same and the deviation weight coefficients were the same. In test scenario No. 2, there was the same offtake change condition as in scenario No. 1, however the water level deviation tolerance of the pools was different. Specifically, the water level deviation tolerance of pools Nos. 1–5 was similar, while the tolerance of pool No. 6 was smaller, so the deviation weights coefficient of pool Nos. 1–5 were set to be 1, with the last pool set at 2. Also, to check the universality of MPC control under this system, another offtake change condition was also tested. In test scenario No. 3, the weight coefficients were the same as in test scenario No. 2, however the offtake flow in Pool 1 was changed.

2.2. Simulation Model

The simulation model of the channel can be described by Saint-Venant equations

$$\{\begin{array}{l}\frac{\partial A}{\partial t}+\frac{\partial Q}{\partial x}=q\\ \frac{\partial }{\partial t}\left(\frac{Q}{A}\right)+\frac{\partial }{\partial x}\left(\frac{Q^2}{2A^2}\right)+g\frac{\partial h}{\partial x}+g\left(S_f-S_0\right)=0\end{array}$$

(1)

where x and t are the space and time coordinates; A is the wetted area (m²); Q is the flow rate (m³/s); h is the water depth (m); S₀ is the canal bottom slope; g is the acceleration of gravity (m/s²); q is the lateral flow rate of the canal for a unit length (m²/s); and S_f is the friction slope, which is defined as

$$S_f=\frac{Q^2{n}^2}{A^2{R}^{4/3}}$$

(2)

with n being the roughness coefficient (s/m^1/3) and R the hydraulic radius (m), defined by R = A/P, where P is the wetted perimeter (m). The implicit difference scheme [18] is adopted to discretize the Saint-Venant equations and the double sweep method can be used to solve the Equations [19].

As the model tuning is a general problem for all deterministic models [20] and former study has shown that automatic control methods have a good robustness [14,21,22], the model tuning is not discussed in this paper. The downstream boundary was set as a constant water level boundary, with a water depth of 3 m. The upstream boundary was similarly set as a constant water level boundary, with a water depth of 7 m. The basic parameters and initial flow condition of each of the pools are shown in

Table 1. Basic parameters and initial flow condition of each pool.

Pool	Pool Length (km)	Bottom Width (m)	Side Slope	Slope	Downstream Initial Flows (m3/s)	Offtake Initial Flowsb (m3/s)	Target Water Depth (m)
Heading					94.5
1	26.6	21	2	9.8 × 10−5	87	7.5	4.5
2	9.7	22.5	2.75	3.9 × 10−5	70	17	4.5
3	14.9	17	1	6.2 × 10−5	42	28	4.21
4	20.8	10	2	5.4 × 10−5	42	0	4.19
5	14.7	7.5	2.5	5.1 × 10−5	42	0	4.21
6	25.4	7.5	2.5	5.3 × 10−5	35	7	3.95

Note: Each pool is composed of many sections. Parameters of the bottom width and side slope are approximate numbers.

. The simulation time step for the hydraulic model of the test canal was 1 min here, and the control time step was 10 min.

2.3. Water Level Difference Control Strategies

In cases where the inflow is fully controlled, the inflow can be flexibly adjusted to maintain constant water levels in the downstream pools. For downstream water level control, a check gate can be adjusted on the basis of the water level deviation at the downstream end of the next pool downstream. The water level deviation $e_j$ is defined as

$$e_j=y_j-S{P}_j$$

(3)

where $y_j$ is the water level at the downstream end of pool j and $S{P}_j$ is the water level set point. In this kind of control strategy, if the pool deviation tolerances differ and people want to keep water level deviations in some pools small, this can be more easily done by setting the corresponding water level weight coefficients matrix Q in the MPC control method to a larger value [23].

In cases where the inflow cannot be controlled, when outflow or offtake flow changes, the downstream pools should behave as storage reservoirs. Control actions are determined on the basis of the difference in water level deviation. Deviation difference $D_j$ is defined as

$$D_j=e_j-e_{j+1}$$

(4)

In this instance, if people want to keep water level deviations in some pools to a certain proportion, it is not feasible to merely set the corresponding weight coefficients in matrix Q in MPC differently. Since it is the water level difference that is controlled here, changing the weight coefficients in matrix Q would change the water level difference, not the water level deviation. So, in order to control the water level deviation, a weight coefficient is added to it. $D_j$ is redefined as

$$D_j=m_j{e}_j-m_{j+1}{e}_{j+1}$$

(5)

where $m_j$ is the weight coefficient of water level deviations $e_j$, which reflects the relative weight of the canal pool. Therefore, the $m_j$ of most pools should preferably be set at 1; that of the important pools, with smaller allowable water deviations, can be set larger.

As it is the redefined $D_j$ that is controlled in this case, the mathematical form of the control system involving state variables, controlled variables, and control action variables should be changed in order to design the proper controller. For control purposes, it is common to use the integrator-delay (ID) model [24] for canal pools. It assumes that a canal reach is separated into a uniform flow with the property delay time and a backwater section with the property storage area.

$$h\left(k+1\right)=h\left(k\right)+\frac{T_s}{A_s}\left\{q_{in}\left(k-k_d\right)-\left[q_{out}\left(k\right)+q_{offtake}\left(k\right)\right]\right\}$$

(6)

where h is the water level at the downstream end of the pool, $q_{in}\left(k-k_d\right)$ is the inflow to the backwater section with delay time steps $k_d$, $q_{out}\left(k\right)$ is the downstream outflow, $q_{offtake}\left(k\right)$ is the off-take outflow, $A_s$ is the average storage area and $T_s$ is the control time step.

As normally the water level deviation is controlled and the flow increment can be more directly controlled by control structures [25], Equation (4) is rewritten in an incremental form with water level deviation:

$$e\left(k+1\right)=e\left(k\right)+\Delta e\left(k\right)+\frac{T_s}{A_s}\left\{\Delta q_{in}\left(k-k_d\right)-\left[\Delta q_{out}\left(k\right)+\Delta q_{offtake}\left(k\right)\right]\right\}$$

(7)

where $\Delta e\left(k\right)$ is the increment of $e\left(k\right)$ with $\Delta e\left(k\right)=e\left(k\right)-e\left(k-1\right)$; also, $\Delta q_{in}\left(k-k_d\right)$, $\Delta q_{out}\left(k\right)$, and $\Delta q_{offtake}\left(k\right)$ are the increment of $q_{in}\left(k-k_d\right)$, $q_{out}\left(k\right)$ and $q_{offtake}\left(k\right)$, respectively.

Substituting Equation (5) into Equation (3), Equation (3) is:

$$\begin{array}{l}{D}_i\left(k+1\right)=D_i\left(k\right)+\Delta D_i\left(k\right)+m_i\frac{T_s}{A_{s,i}}\Delta q_{in}\left(k-k_{d,i}\right)-m_i\frac{T_s}{A_{s,i}}\left[\Delta q_{c,i}\left(k\right)+\Delta q_{offtake,i}\left(k\right)\right]\\ -m_{i+1}\frac{T_s}{A_{s,i+1}}\Delta q_{c,i}\left(k-k_{d,i+1}\right)+m_{i+1}\frac{T_s}{A_{s,i+1}}\left[\Delta q_{c,i+1}\left(k\right)+\Delta q_{offtake,i+1}\left(k\right)\right]\end{array}$$

(8)

where $\Delta D\left(k\right)$ is the increment of $D\left(k\right)$. $\Delta q_{c,i}\left(k\right)$ is the increment of outflow $q_{c,i}\left(k\right)$ of pool i. Equation (8) expresses the relationship between gate flow and the controlled water level difference $D\left(k\right)$ of two adjacent pools.

To establish controller design, the ID model of pools should first be determined. The two characteristics of the ID model, delay time, and storage area were calculated by applying the system identification technique [26,27]. The values are shown in

Table 2. Basic parameters and initial flow condition of each pool.

	Pool	Pool 1	Pool 2	Pool 3	Pool 4	Pool 5	Pool 6
Characteristics
As (m2)		582,524	441,176	327,869	447,761	361,446	431,655
Td (min)		70	24	35	57	41	75

.

2.4. Model Predictive Control

MPC is a control strategy that explicitly uses a simplified process model of the real system to obtain control actions by minimizing an objective function. MPC has three basic components, including a process model, an objective function, and a rolling optimization strategy [28]. A process model is used to predict the system output for some time into the future. Normally a linear invariant state-space model is used as a process model in the canal control, with the form

$$x\left(k+1\right)=Ax\left(k\right)+B_{u}u\left(k\right)+B_{d}d\left(k\right)$$

(9)

$$y\left(k\right)=Cx\left(k\right)$$

(10)

where x is the state vector; y represents the output variables of the modeled water system, which means the D(k) of all pools here; u is the vector of input variables calculated by the controller, which is the $\Delta q\left(k\right)$ of the intermediate check gate here; and d is the vector of known measurable disturbances at time step k. A represents the system matrix, $B_u$ is the input to state matrix, $B_d$ is the disturbance to state matrix, and C is the state to output matrix. The control time step is 10 min, so the delay time steps of the pools are 7, 2, 4, 6, 4, and 8, respectively. Then, Equation (6) was used for controller design to obtain the process model, with $x_{34\times 1}$ the state vector, $y_{5\times 1}$ the output vector, $u_{5\times 1}$ the input vector, $d_{6\times 1}$ the disturbance vector, $A_{34\times 34}$ the system matrix, $B_u{}_{34\times 5}$ the input to state matrix, $B_d{}_{34\times 6}$ the disturbance to state matrix, and $C_{5\times 34}$ the state to output matrix. The MPC control and LQR control could then all be designed on this process model.

Then, an output prediction is made using the process model. In the prediction process, the predicted output of the system, $y\left(k+i|k\right)$, is determined from the current state vector, $x\left(k\right)$, and the future control actions, $u\left(k+i|k\right)$. The counter i indicates the number of time steps into the future. In $y\left(k+i|k\right)$, the counter i ranges from 1 to p, the prediction horizon. While in $u\left(k+i|k\right)$, the counter i ranges from 1 to m, the control horizon should be less than or equal to p. It is assumed that there are control actions from current time step k to future time step k+m and no control actions from time step k+m+1 to future time step k+p.

The predicted values for the state and output vectors one time step into the future are expressed as

$$x\left(k+1|k\right)=Ax\left(k\right)+Bu\left(k\right)+Dd\left(k\right)$$

(11)

$$y\left(k+1|k\right)=Cx\left(k+1|k\right)=C\left[Ax\left(k\right)+Bu\left(k\right)+Dd\left(k\right)\right]$$

(12)

The prediction for the state vector and output vector two time steps into the future are

$$\begin{array}{l}x\left(k+2|k\right)=Ax\left(k+1|k\right)+Bu\left(k+1\right)+Dd\left(k+1\right)\\ =A\left[Ax\left(k\right)+Bu\left(k\right)+Dd\left(k\right)\right]+Bu\left(k+1\right)+Dd\left(k+1\right)\\ =A^{2}x\left(k\right)+ABu\left(k\right)+Bu\left(k+1\right)+ADd\left(k\right)+Dd\left(k+1\right)\end{array}$$

(13)

$$\begin{array}{l}y\left(k+2|k\right)=Cx\left(k+2|k\right)\\ =C\left[A^{2}x\left(k\right)+ABu\left(k\right)+Bu\left(k+1\right)+ADd\left(k\right)+Dd\left(k+1\right)\right]\end{array}$$

(14)

This process continues until the end of the control horizon, m, is reached. The state vector and output state vector are

$$\begin{array}{ll}\hfill x\left(k+m|k\right)& =A^{m}x\left(k\right)+A^{m-1}Bu\left(k\right)+A^{m-2}Bu\left(k+1\right)+\cdots +ABu\left(k+m-2\right)+Bu\left(k+m-1\right)\\ & +A^{m-1}Dd\left(k\right)+A^{m-2}Dd\left(k+1\right)+\cdots +ADd\left(k+m-2\right)+d\left(k+m-1\right)\\ & =A^{m}x\left(k\right)+{\displaystyle \sum _{i=1}^m{A}^{m-i}Bu\left(k+i-1\right)}+{\displaystyle \sum _{j=1}^m{A}^{m-j}Dd\left(k+j-1\right)}\end{array}$$

(15)

$$y\left(k+m|k\right)=Cx\left(k+m|k\right)=C\left[A^{m}x\left(k\right)+{\displaystyle \sum _{i=1}^m{A}^{m-i}Bu\left(k+i-1\right)}+{\displaystyle \sum _{j=1}^m{A}^{m-j}Dd\left(k+j-1\right)}\right]$$

(16)

After the control horizon has passed, the remaining output predictions are based on the free response only. At the end of the prediction horizon, p, the predicted state vector and output state vector are

$$\begin{array}{ll}\hfill x\left(k+p|k\right)& =A^{p}x\left(k\right)+A^{p-1}Bu\left(k\right)+A^{p-2}Bu\left(k+1\right)+\cdots +A^{p-m}Bu\left(k+m-1\right)\\ & +A^{p-1}Dd\left(k\right)+A^{p-2}Dd\left(k+1\right)+\cdots +ADd\left(k+p-2\right)+Dd\left(k+p-1\right)\\ & =A^{p}x\left(k\right)+{\displaystyle \sum _{i=1}^m{A}^{p-i}Bu\left(k+i-1\right)}+{\displaystyle \sum _{j=1}^p{A}^{p-j}Dd\left(k+j-1\right)}\end{array}$$

(17)

$$y\left(k+p|k\right)=Cx\left(k+p|k\right)=C\left[A^{p}x\left(k\right)+{\displaystyle \sum _{i=1}^m{A}^{p-i}Bu\left(k+i-1\right)}+{\displaystyle \sum _{j=1}^p{A}^{p-j}Dd\left(k+j-1\right)}\right]$$

(18)

An objective function, which is typically a combination of errors of output variables between a given reference and control actions over the prediction horizon, is minimized by adjusting future control actions. It is possible to assume that the output reference is that the $D\left(k\right)$ equals zero. Hence, the objective function can be expressed as

$$\underset{u(k)}{\min }J={\displaystyle \sum _{j=0}^p\left(y^T\left(k+j|k\right)Qy\left(k+j|k\right)\right)}+{\displaystyle \sum _{j=0}^{m-1}\left(u^T\left(k+j|k\right)Ru\left(k+j|k\right)\right)}$$

(19)

where Q is the weighting matrix of output and R is the weighting matrix of input. The problem can be summarized as minimizing the objective function by adjusting the future control actions u(k). Once the sequence of future control actions are determined, only the first set of control actions are implemented on the irrigation system. The system is then updated and the process repeated. This is the rolling optimization strategy of an MPC controller.

As an LQR control method is compared later, the background of LQR controllers is also introduced. The objective function of the LQR method is an infinite time domain equation as

$$\min J={\displaystyle \sum _{j=0}^{\infty }\left(y^T\left(j\right)Qy\left(j\right)+u^T\left(j\right)Ru\left(j\right)\right)}$$

(20)

So, the optimization is done in an infinite time domain and the process model is a state-space model with the form

$$x\left(k+1\right)=Ax\left(k\right)+B_{u}u\left(k\right)$$

(21)

The meaning of the parameters in Equation (21) is the same as that in Equation (8), however in Equation (21), the process model does not contain the term $B_d$ and $d$, which means that the optimization process does not use the future disturbance information in the LQR control.

For MPC, many times tuning is done through trial-and-error techniques [29]. There were still four parameters to be determined for MPC control: the prediction horizon p, the control horizon m, the cost weighting matrix Q for outputs, and the cost weighting matrix R for the control actions. The prediction horizon, p, should be long enough to include all necessary dynamics of the system. As the total delay time steps of the upstream disturbance on the water level of pool 6 was 24 (sum of the delay steps of pool 2 to 6), the prediction horizon p was set to 30 and the control horizon m was 20. Values within Q and R provide a trade-off between minimizing water level differences and minimizing check flow changes, so the matrix R can be an identity matrix. If R is relatively greater than Q, the controller will focus more on the minimization of the gate control actions than on water level difference. Therefore, a larger value of Q is preferred here. Note that a much larger value of Q may lead to instability in the control system, as the controller is designed on a simplified process model and great changes in the control actions result in great pool response and sometimes great fluctuations and even resonance. The value of the Q should be chosen carefully and in this simulation, Q was set as a diagonal matrix with an elements value of 5 for the MPC control method.

The two cost weighting matrices Q and R were also needed for LQR control. For comparison purposes, the Q matrix was also a diagonal matrix with an elements value of 5 for LQR control. However, as the results of LQR control with this Q matrix did not perform well in difference control in the Results part, a diagonal matrix Q with an elements value of 30 was also used in LQR control. The LQR control with matrix Q with an elements value of 5 is referred to as LQR-I and the LQR control with matrix Q with an elements value of 30 is referred to as LQR-II. Also, as the LQR method does not use any future disturbance information, in order to get a fair comparison of the MPC method with the LQR method in water level difference control, in a condition, the future disturbance term $d\left(k+j-1\right)$ was not used in the MPC control method. In this condition, the future disturbance was set to be zero in the process model. So, in this paper, the MPC control was used in two conditions, with and without future disturbance information. The MPC control with future disturbance information used is referred to as MPC-I and the MPC control with no future disturbance information used is referred to as MPC -II in this paper.

The MPC control method, including MPC-I and MPC-II, the LQR-I control method, and the LQR-II control method are all used with the difference control strategy for the test scenarios.

2.5. Performance Indicators

Several indicators are used here to show the performance of the controllers such as maximum absolute water level deviation (MAE), average absolute water level deviation (AAE), maximum absolute water level difference (MAD), average absolute water level difference (AAD), and the time gate flow changes. MAE, AAE, MAD, and AAD are defined as

$$\mathrm{MAE}=\max \left(\left|e_j\right|\right)$$

(22)

$$\mathrm{AAE}=\frac{\Delta t}{T}{\displaystyle \sum _{t=0}^T\left(\left|e_j\right|\right)}$$

(23)

$$\mathrm{MAD}=\max \left(\left|m_j{e}_j-m_{j+1}{e}_{j+1}\right|\right)$$

(24)

$$\mathrm{AAD}=\frac{\Delta t}{T}{\displaystyle \sum _{t=0}^T\left(\left|m_j{e}_j-m_{j+1}{e}_{j+1}\right|\right)}$$

(25)

where t is the time disturbance that happens, T is the time period from the occurrence of disturbances to the end of simulation.

3. Results

The simulation results for the test scenarios are presented in Figure 2, Figure 3, Figure 4, Figure 5, Figure 6 and Figure 7. A summary of all simulations is found in

Table 3. Summary of simulation results.

Indicators		Scenario 1				Scenario 2				Scenario 3
		MPC-I	MPC-II	LQR-I	LQR-II	MPC-I	MPC-II	LQR-I	LQR-II	MPC-I	MPC-II	LQR-I	LQR-II
MAE (m)	Pool 1	0.144	0.144	0.182	0.180	0.187	0.185	0.209	0.210	0.167	0.168	0.231	0.197
	Pool 2	0.144	0.142	0.177	0.171	0.186	0.184	0.205	0.201	0.165	0.166	0.216	0.186
	Pool 3	0.133	0.130	0.148	0.150	0.173	0.169	0.174	0.177	0.155	0.155	0.170	0.159
	Pool 4	0.118	0.115	0.112	0.124	0.156	0.150	0.135	0.150	0.138	0.138	0.117	0.128
	Pool 5	0.104	0.101	0.088	0.098	0.139	0.133	0.113	0.124	0.114	0.114	0.073	0.096
	Pool 6	0.092	0.089	0.090	0.081	0.063	0.080	0.082	0.080	0.047	0.047	0.023	0.037
AAE (m)	Pool 1	0.079	0.073	0.096	0.118	0.107	0.098	0.107	0.130	0.109	0.115	0.144	0.128
	Pool 2	0.079	0.078	0.097	0.113	0.106	0.104	0.109	0.127	0.098	0.106	0.128	0.115
	Pool 3	0.074	0.077	0.085	0.095	0.098	0.102	0.098	0.109	0.074	0.083	0.088	0.085
	Pool 4	0.070	0.074	0.068	0.075	0.090	0.096	0.082	0.088	0.052	0.058	0.046	0.053
	Pool 5	0.068	0.076	0.064	0.063	0.086	0.099	0.086	0.081	0.035	0.038	0.022	0.031
	Pool 6	0.071	0.081	0.072	0.059	0.047	0.056	0.057	0.045	0.013	0.012	0.005	0.009
MAD (m)	Pool 1	0.012	0.013	0.006	0.010	0.014	0.015	0.006	0.010	0.014	0.018	0.020	0.019
	Pool 2	0.018	0.019	0.029	0.030	0.020	0.021	0.031	0.034	0.032	0.036	0.050	0.040
	Pool 3	0.016	0.016	0.036	0.031	0.020	0.021	0.039	0.032	0.034	0.038	0.059	0.042
	Pool 4	0.014	0.021	0.025	0.027	0.019	0.035	0.025	0.029	0.029	0.031	0.045	0.034
	Pool 5	0.027	0.041	0.050	0.036	0.051	0.102	0.119	0.100	0.024	0.025	0.028	0.024
	Pool 6	-	-	-	-	-	-	-	-	-	-	-	-
AAD (m)	Pool 1	0.006	0.006	0.003	0.005	0.006	0.006	0.003	0.005	0.005	0.009	0.016	0.013
	Pool 2	0.010	0.011	0.015	0.021	0.013	0.013	0.015	0.021	0.020	0.023	0.040	0.031
	Pool 3	0.010	0.011	0.019	0.022	0.013	0.014	0.020	0.023	0.021	0.025	0.043	0.031
	Pool 4	0.006	0.009	0.011	0.015	0.008	0.013	0.014	0.017	0.017	0.020	0.024	0.022
	Pool 5	0.006	0.010	0.017	0.015	0.012	0.021	0.034	0.023	0.012	0.014	0.012	0.013
	Pool 6	-	-	-	-	-	-	-	-	-	-	-	-
T (h)	Gate 1	6.5	10.2	10.3	10.3	6.5	10.2	10.3	10.3	5.7	10.2	10.2	10.2
	Gate 2	6.6	10.2	10.3	10.3	6	10.2	10.3	10.3	5.7	10.3	10.5	10.3
	Gate 3	5.7	10.2	10.3	10.3	5.3	10.2	10.3	10.3	6.2	10.3	10.5	10.5
	Gate 4	5.2	10.2	10.3	10.3	5.2	10.2	10.3	10.3	6.8	10.5	10.5	10.5
	Gate 5	5.5	10.2	10.2	10.2	5.5	10.2	10.2	10.2	7.8	10.5	10.5	10.5

Note: The font of the maximum value of each indicator in all pools is marked in red.

.

Figure 2 and Figure 3 show the results for scenario No. 1, in which the same coefficient of all pools was used. In Figure 2, the water levels in all pools tend to decrease at a similar rate with the MPC-I, MPC-II, LQR-I, and LQR-II methods. However, after 20 h, the water level deviations were smaller at the downstream pools compared with those at the upstream using both the MPC and LQR methods. That is because although the gates at both ends are not controlled, the inflow increases and outflow decreases as the water levels decrease. The decrease in outflow was greater than the increase in inflow, so the water levels decrease trend was gentle in the downstream pools. As can be seen from Figure 3, the flow of Gate 0 (the head gate) and Gate 6 (the last-gate) also changed, however the flow change of Gate 6 was obviously greater than that of Gate 0.

By comparing the result of Figure 2a–d, it can be seen that the MPC control performed better at minimizing the water level deviation difference of adjacent pools, with a maximum AAD of 0.006m with the MPC-I control method and 0.010m with the MPC-II control method. The maximum value of the LQR-I method was 0.017 m and that of the LQR-II method was 0.015 m. Consequently, there was a smaller maximum water level deviation in all pools—0.144 m with both MPC methods compared with 0.182 m with the LQR-I control method and 0.180 m with the LQR-II control method. With better control of water level differences, the imbalance in flow was better distributed so that the maximum water level deviation caused by flow change among all pools was small. In MPC-I control, the future disturbance was considered and control actions were taken before the disturbance happened. The control result was expected to be better than LQR control. However, in MPC-II control, no future disturbance information was used, however the results were similar to those with MPC-I control and better than those in LQR methods. Two reasons can account for this. One reason is that in the MPC control, the function J minimizes the control actions in control horizon m and the water level differences in prediction horizon p, greater than m, so more water level differences are considered and the resultant control actions of the MPC method are greater to better minimize the water level differences compared with LQR-I control, which has the same weighting matrix Q. The other reason is that MPC control minimize the function in a finite time domain compared with LQR control in an infinite time domain and MPC obtains the local optimal solution while LQR obtains the global optimal solution and it is supposed that MPC control performs better than LQR in the early stage after the disturbance occurred where there is a severe water level change, however the opposite is thought to occur in the later stage. The results of Figure 2b,d also support that MPC performed better in the early stage. However, in the later stage, as the water level in all the pools decreased, the inflow of the head-gate increased and the outflow of the last gate decreased, which was not considered in the process model, and the water levels in all pools tended to stabilize at another water level rather than keep decreasing and water level differences were also controlled to zero faster than the situation where the inflow and outflow are constant. A local optimization performs better than global optimization. So, the result of Figure 2b was still better than Figure 2d.

By comparing the MPC-I results with MPC-II results in Figure 2 and Figure 3, it was found that the MPC control use of the future disturbance information performs better than MPC control without future disturbance information, however this is not so obvious. The MAD of MPC-I in pool 5 was 0.027 m compared with that value of 0.041 m of MPC-II. In the MPC-I method, control actions were taken before the disturbance happened. In Figure 3a, the gate flows of Gate 1–5 began to change at the time 6.5 h, 6.6 h, 5.7 h, 5.2 h, 5.5 h. All are ahead of 10 h, when the offtake flow changes, while in Figure 3b, the gate flows all begin to change around 10 h, only after the disturbance happens. Although advance actions were taken, the feed-forward flow changes were small.

Also, by comparing Figure 2c,d, this shows that LQR-II was better than LQR-I in minimizing water level difference. That is because in LQR-II, the values of the elements of the weighting matrix Q are greater. However, there were also greater water level fluctuations and flow change fluctuations in Figure 2d and Figure 3d as compared with Figure 2c and Figure 3c. As greater values in matrix Q may result in significant flow change and water level fluctuations, they are more likely to cause water level resonance, so the values are set at 30 here instead of at a larger value. This also shows that by increasing the values in matrix Q in the LQR method to get a better control of water level difference, this is more dangerous compared with the MPC control.

Figure 4 and Figure 5 show the results for scenario No. 2 in which the weight coefficient of Pool 6 was set at 2 while the others were set at 1. Figure 3 shows that with the weight coefficients used, the water level deviation in Pool 6 was smaller than that in the other pools, almost half of the water-level deviation in Pool 5. In the MPC-I control method, the maximum water level deviation in Pool 6 was 0.063 m and was 0.139 m in Pool 5 in the MPC-II control method. The maximum water level deviation in Pool 6 was 0.080 m and 0.133 m in Pool 5. With the LQR-I method, the maximum water level deviations in these two pools were 0.082 m and 0.113 m, respectively, and with the LQR-II method, 0.082 m and 0.124 m, respectively. This means that with the proposed strategy, the water level deviations can be a certain proportion as we hope, and by using Equation (10) to build and design the controllers, the water level and flow control processes are reasonable. The water level deviations in Pools 1–5 gradually approached the same value, while the water level in Pool 6 approached another value in both the MPC and LQR methods. Compared with results for scenario No. 1, the water level deviations of Pools 1–5 were much bigger. The maximum water level deviations of Pools 1–5 were 0.187 m, 0.185 m, 0.209 m, and 0.210 m with the MPC-I control method, the MPC-II control method, the LQR-I method and the LQR-II method in scenario No. 2, and larger than those in scenario No. 1 with values of 0.145 m, 0.143 m, 0.182 m, and 0.180 m. That is because the small water level deviation in Pool 6 in scenario No. 2 requires small inflow and outflow imbalance in Pool 6, while the flow imbalance in other pools is greater, so the water levels in these pools decreased more.

Also, in scenario No. 2, both MPC control performed better than LQR control methods. This time, the MPC-I control performed better than MPC-II control did obviously. The MAD in pool 5 is 0.051 m in MPC-I control compared with that of 0.102 m in MPC-II control. As the weight coefficient of Pool 6 was set at 2, much more flow changes were required in the upper pools than in scenario No.1. Feed-forward flow changes happened in MPC-I and the results were better. It can be concluded that the more the flow changes, the more obvious the advantages of the MPC-I method will be.

In Figure 5d, the flow fluctuations in Pools 1–5 are much more obvious using the LQR-II method in scenario No. 2 than in scenario No. 1 because the greater flow increase of Gate 5 was required to reduce the decreased rate of the water level in Pool 6 as compared with scenario No. 1. The maximum flow change of Gate 5 was 4.8 m³/s and 5.2 m³/s in scenario No. 1 and scenario No. 2, respectively.

Figure 6 and Figure 7 demonstrate the control results for scenario No. 3, in which the disturbance occurred in Pool 1. Although the water disturbance happened in Pool 1, the water level deviation in Pool 1 decreased quickly when the disturbance happened at 10 h, however at 35 h, the water level deviation in Pools 2–5 was close to the water level deviation in Pool 1, and the water level deviation in Pool 6 is about half of that with the MPC control methods. Similar to previous results, both MPC control performed better than LQR control methods. The maximum MAE in MPC-I, MPC-II, LQR-I, and LQR-II control are 0.167 m, 0.168 m, 0.231 m, and 0.197 m, respectively. The maximum MAD in MPC-I, MPC-II, LQR-I, and LQR-II control are 0.034 m, 0.038 m, 0.059 m, and 0.042 m, respectively.

Forward control actions were taken with the MPC-I control method. In Figure 7a, the gate flow of Gate 1–5 began to change at time 5.7 h, 5.7 h, 6.2 h, 6.8 h, and 7.8 h, respectively, with the MPC-I control method, while in Figure 7b with the MPC-II control method, the times were 10.2 h, 10.5 h, 10.5 h, 10.5 h, and 10.5 h, respectively. Although the feed-forward controls were taken, the flow changes were all particularly small as feed-forward change may cause the water level to increase, which may increase the water level differences on the contrary. So, the MPC-I control performed almost the same as the MPC-II control.

The LQR-II control also performed better than the LQR-I control, however the flow fluctuation is also much more obvious in the LQR-II control in Figure 7d. The results of the LQR-II control were much more similar to the MPC-II control compared to the LQR-I control with the same Q with the MPC-II control.

4. Conclusions

In this paper, a water-level difference control strategy was improved to keep the water level deviations in all pools to certain proportions in a particular situation where the operator does not have full control over the canal inflow. Then, the control strategy was conducted with MPC control and LQR methods, respectively, to get an automatic control of the water level. A simulation model of a canal system was built to show the control results. There are several important conclusions that can be made based on the various simulation results. These conclusions are summarized below.

• Water level difference control allows the operator to have no full control of the head gate and the tailgate while automatically controlling all check gates in between. It reveals flow mismatches by causing the water levels to rise or fall at the same rate.
• By adding a weight coefficient to the water level deviation to construct water level difference and with several changes in controller design, the control method can make the water levels rise or fall at different rates in the proportion that people want with flow mismatches, consequently changing the water level deviations with the proportion.
• Both the LQR and MPC control methods with the proposed control strategy work to minimize the water level difference, however the MPC control performs better even with no future disturbance information taken into account as a local optimization is better than global optimization in LQR control in water level difference control.
• The MPC control method performs better when future disturbance information is taken into account and can take feed-forward control before disturbance happens. However, the more upstream the disturbance occurs, the less obvious this advantage is.