Linear Quadratic Optimal Controller Design for Constant Downstream Water-Level PI Feedback Control of Open-Canal Systems

The key point of PI feedback control is how to design appropriate controller parameters. This paper combined the linear quadratic optimizing control theory to design an online controller for constant downstream water-level operation which could respond to different working condition properly and rapidly avoiding complex controller parameters adjusting. Based on the Integrator-Delay (ID) model simplified from Saint-Venant Equations, we established the discretized linear time invariant system of canals. Then transferred it into the state-space equations and obtained the state-feedback equation. Water level deviations and flow rate increments were chosen to form the objective function and this paper recommended values of Q and R weight matrices among it should be set according to the optimum quadratic form indicators correspondingly. This controller was applied to two practical canal systems which had diverse scales. Results showed the system under control quickly regained stability; optimizing the objective function with recommended weight matrices could well balance demands on water level deviations and flow rate changes; dynamic performances of water movements and gate movements were acceptable. Through simulations, we preliminarily proved the practicability of this online PI controller implementing LQR. This work proposed an available solutions for the design and operation of water conveyance systems around the world.


Introduction
Inter-basin water diversion projects are aiming at easing the contradictions of uneven distribution of water resources in space. Automation canal control exerts a gradually important part in water deliver projects for its high efficiency, low labour cost and satisfactory performances. The Central Arizona Project (CAP) in America develops an operation mode of control volume to avoid delay in long-distance channel and easily realize centralized control [1,2]. Control systems of Broken River in Australia are designed based on a system identification model, being effective for water saving and environmental benefits [3]. The Middle Route of Southto-North Water Diversion Project in China has successfully transported several million cubic meters of water to northern China under constant downstream water-level control [4][5][6][7].
Proportional-integral (PI) feedback algorithm is an important issue in canal control algorithms due to the extensive applicability, high robustness and good performance, but tuning of PI parameters is challenging. Over the past few decades, researchers have put forward many methods to get suitable controller parameters K p and T i [8][9][10][11][12][13][14][15][16][17]. Typically, for practical engineering, the mainly used methods are divided into two kinds: engineering tuning method and trial method. While the latter depends much on trial and error, the former have clearer principle and easier to achieve. Ziegler and Nichols [18] had completed a mature heuristic tuning called Ziegler-Nichols method though the step response curve. Litrico et al. [19] advised the Auto Tune Variation (ATV) method utilizing gain margin and phase margin of the system itself to acquire controller parameters. However when confronting with canal systems with variable cross sections and complicated through-water buildings, they might not fit so well. Nowadays with the rapid development of control theories, optimal theory, selfadaptation, fuzzy control, model predictive control, neural network and other methods are also employed to tune controller parameters.
Tuning parameters of the PI controller with optimal theory is becoming more and more intensively favored in practical applications. Baume et al. [9] proposed a global approach based on hydraulic model of full Saint-Venant equations, however the non-linear optimization problem was so complex that it needed extra independent variable numbers. Meanwhile some scholars suggested that instead of non-linear form, linear quadratic optimal design should be applied to the canal hydraulic model. Malaterre et al. [15] developed a linear model discretized Saint-Venant equations through intricate finite difference approximation. Besides, they chose upstream discharge and check gate opening as control action variables, the linearization error of the relationship between them would be fairly small only within certain realms. Wahlin et al. [20,21] put forward a fully centralized PI controller using linear quadratic regulator on the basis of linear model -ID model, which is much more concise and easier to linearization, results showed the algorithm could deal with the planned demands for large quantities of water. Their work provides a lot of valuable reference for us, however the weight matrices Q and R in LQR was obtained through trial, which may cost considerable work and rely much on workers' experience.
Thus adopted the logic of Wahlin et al. [20,21], this paper try to make a small step forward: choosing ID model as a linear hydraulic model to simulate canal flows for its conciseness and certain precision, and forming the state space expressions, then designing the PI controller combined with linear quadratic optimal theory, in which making full use of canal response characteristics to determine suitable weight matrices, so that we are capable of getting acceptable controller parameters avoiding repeated tuning.
The main content of this paper is as follows. Firstly we establish the discrete canal model based on the integratordelay (ID) model in Section 2. Then transfer the model into a set of state space equations using state space methods and obtain a state feedback equation in Section 3. Thirdly in Section 4 we implant the LQR optimal theory and based on analysis on simulated canal response characteristics, we propose a method determining the weight matrices of the objective function in LQR. Lastly on the canal control simulation platform [22], we apply the controller to two projects as test cases demonstrating its reliability.
This study is concentrating on constant downstream water-level feedback control. In other words, canals attempt to maintain a constant water level at the downstream position and it matches requirement decided by downstream. In addition, the canals mentioned are all flat with subcritical flow, and working conditions are mild. Compared with steep slopes, flat gradients have wider application in practical cases, for its being more inclined towards mild flow and preventing instability from some unreasonable managements and unknown disturbances.

Linear Canal Control Model
Integrator-Delay (ID) model makes it possible for us to discuss linearized the hydraulic process. ID model which serves as a linearized approximation of canal flow is used to define discretized state-transition equations in order to match LQR.
The model divides the canal into two parts: one part with uniform flow and another part called backwater. Linearization and Laplace transform of Saint-Venant equations was done by Schuurmans et al. [23,24], and the expression of ID model can be expressed as below: where h = water level in the backwater part ; A s = surface of the backwater area in canals; q in = upstream inflow; q out = downstream outflow; τ = delay time of canals. The all variables denote that they are around the steady-state.
For canals under constant downstream water level, we set the target water level for the downstream point, and ordinarily it also treats as initial water level. So we get another expression: where e(t) = water level deviation of backwater part, Main parameters in ID model are the area of backwater (A s ) and delay time (τ ), which could be gotten by system identification or calculation of physical properties, or mutual authentication of these two methods furthermore. In this paper the working conditions we focused on is that discharge only changes on a small scale, therefore A s of each canal can be treated as constants. The backwater area is integrated within the region of backwater and is approximately equal to the slope of the relationship curve between h(t) ~ t . Here we mainly introduce an identified method -we plan to give a square wave to upstream flow rate (Figure 1), then record the downstream water level process and fit it (Figure 2). The slope of fitting lines is considered as 1/A s and there is an intersection of the fitting line and original water level, so τ gets, that is the difference between the intersection and planed change point in time series.  The flow in canals we discussed appears as the subcritical flow -one kind of the most common flows in canals. So we assume waves only propagate downstream. According to kinematic theory, velocity could be calculated as gravitational waves. It's a just theoretical value but it could validate the identified value d where min = minimum delay time; L n = length of backwater area; g = the gravitational acceleration, 9.8 m 2 /s, d = water depth.

State Space Expression
Through state space method, a linear system can clearly describe the causal link between state variables, control variables and output variables. And it is not only available for the single-input-single-output system, but also for the multiple-input-multiple-output system. With application of matrix and vector space theory, accompanied time domain method studying dynamic characteristics of canal systems, we transferred the ID model into discrete state space equations.

State Variables
The first issue to accomplish the transformation is to choose suitable state variables. For canal systems, we basically control the water level to achieve our management goals. Besides the PI controller makes response to the water level deviation and the change rate of it. So water level information, including water level deviation and the change relative to last time, was selected to represent the canal features. The state vector of a canal system x(k) and among it, concrete contents xj(k) characterizing Canal j, could be shown as below: where k = the number of time step; e(k) = water level deviation of backwater part for a discrete system, e(k)=h(k)-h target ; e(k) = the change of water level deviation against last moment; n j+1 = the integer portion of delay time to times step. We use subscript to show canal number.

PI Control Algorithm
The open-canal system in this paper is under the operation of constant downstream water-level, and the feedback controller is PI controller. Here we employ the incremental PI controller. Control performances of PI controller are conditional on its parameters, thus we combined the optimal theory to adjust its parameters, that is K p and T i .
where ∆u(k) = the change of flow rate against last time step as controller instructs, ∆u(k)=Q(k)-Q(k-1); K p = the proportion gain; T i = the integration time constant.

State Feedback Equation
Δu and ΔG (the change of gate opening against last time step) are two obvious choices for control variable. Downstream gates of each canal receive the control instruction and then adjust their opening. Gate opening is the final variable to realize controlling the flow rate. However directly taking it as the control variable couldn't fit well. For one reason, linear quadratic optimal design requires all relationships between the variables have linear expressions. Eq. 5 reveals clearly that output results of a PI controller is the flow rate changes and as we all know, the relationship between flow rate and gate opening can't be simply described as linear, or the error would exceed expectations. So we choose the intervening variable -the flow rate as control variable, then the gate opening was calculated according the sluice flow computing formula modular outside the frame of the controller, which would not be limited by the linear demand.
where ∆u(k) = the vector of combination of ∆u(k); K = the control input, which is an aggregate of each canal's feedback controller parameters K p and T i .

Output Variables
Output variables are the water level deviation at last step. Because it is the main control target under the canal operation of constant downstream water level and could clearly show system control performance.
[ ] where y(k) = the output vector at k.

State Space Equation of ID Model
Thus combined the chosen variables with discrete ID model, state space representation shows: , 0 where; A, B and C are matrices, A = the system matrix representing the link between internal state variables of a system, B = the input matrix indicating how the input variables control state variables. C = the output matrix indicating how the output variable reacts with the state variable.
In discrete model, the sampling time, Δt, matters a lot, influcing the accuracy and stability of discretizing the continuous system. Here we try to eliminate this influence, thus Δt is chosen to be small enough.

Optimization Indicators of LQR
We choose water level deviations and flow rate increments as two evaluation ingredients of controller performances. They are most fundamental to assess but competing to some extent. The quadratic objective function is aimed at selecting the appropriate control input K by minimizing it under the constraint of linear system. Based on the state-transition equations given above, the objective function can be expressed as where Q and R = the diagonal arrays.
There is little doubt that weight matrices Q and R play an important role in the function. A set of weight matrices Q and R provides a trade-off between water deviations and flow rate changes and is the decisive component of K.
Since the relative value of Q and R influences the controller choosing, researchers usually set Q as an identity matrix and adjust R to balance their respective weights. This way makes matrices choosing convenient and the focus was transferred into choosing of R. Here we suggest an engineering tuning method of R choosing based on response characteristics of canals.
Two quadratic form indicators were introduced -nondimensional integrated square of discharge change (NISQ) and non-dimensional integrated square of error (NISE) [25], assessing fluctuations of discharges and water levels respectively. They could be gotten through simulation modelling on computer physical model, no field trials was necessary. Purpose of the objective function is to balance demands on them. So we take optimal values respectively as the basis. Besides optimization target function appears as a non-dimensional quantity, we also do some conversions, NISQ and NISE need to be multiplied by their dimensional criterion, for a canal, that is the target water level and designed flow rate. The weight matrices Q and R show as below: where min (NISE) = the minimum NISE of canal systems; min (NISQ) = the minimum NISQ of canal systems; Δt = a time step; T = simulation time; R 1 = the weighting value for the first canal. The weight value of followed canals, R i , could be set according to relative size of designed flow rate simplistically.
To summarize, the structure of this designed canal system controller shows as Figure 3. MATLAB offers a convenient mathematical tool for LQR. We build the canal control simulation platform on MATLAB and it effectively copes with different canal systems and working conditions.

Physical parameters of target canals
The controller was adapted to two gradually sloping canal systems, which had diverse scales. One case is the simplify model picking up from a large-scale project in China and the other is the ASCE test canal in California, the former's size (taking the flow rate as a typical parameter to evaluate their size) was more than ten times as bigger as the latter. And the working conditions concentrate on flow rate changing at a certain range, the largest change no more than 35% of designed flow rates. Water intakes in odd number canals located downstream start pumping water from canals simultaneously and schedule in time to reach Q out . Sketch maps show as Figure 4, specific physical parameters of each canal can be found in Table 1 ~ Table 2, also working conditions in Table 3 ~ Table 4.

Simulation results of canal performance under designed controller
Following the controller designed below, R 1 is 0.002 in Canal System 1, and is 0.44 in Canal System 2, under the examples' working conditions, simulated canal responses show as Figure 5. We also conclude relevant performance indicates to evaluate performances in a variety of aspects, including the water level, flow rate, gate movement and stable time (   MAE is the maximum absolute error, StE is the steady-state error, IAE is the integral of absolute magnitude of error, NIAQ is the non-dimensional integrated absolute discharge change, NIAW is the non-dimensional integrated absolute gate movement; Stable Time means, at that time, water level deviation is no more than 2% of target and gate opening changes no more than 1 mm.
On the whole, the controller fits well. In the transition time, it maintains strong over-damping features, and the water level and flow rate in canal fluctuate slowly. The weight matrices in optimization target function provide us with a balanced control effect of the water level and flow rate in canals, which both behave gentle rolls and are inclined to be stable. There exist more fluctuates in Canal System 1, we owe it to the specificity of canal system itself: head difference around gates is much smaller and looks more vulnerable to flow rate demands. Besides, the largest overshoot of gate movements is no more than 30 cm and it avoids excessively repeated regulation. The systems regain stability in 24 hours both showing excellent celerity and accuracy. Moreover, the final error between stable state and target state is within an acceptable range. For Canal System 2 due to the light control effect, StE couldn't be exactly eliminated during our simulation, around -2 cm at the last moment

Influence of Q and R weight matrices
The optimized controller parameters under different R 1 are summarized in Table 6. R 1 determines the controller parameters and it clearly indicates that smaller R 1 produces much more aggressive controllers: K p takes on ascend trend, meanwhile T i appropriately cuts down to avert the controller out-of-balance. But here a couple of interesting results are worth mentioning. The details of K p and T i are somewhat different in different canals. For instance, firstly the change rule of K p doesn't represent as strict linear descending curve as R 1 grows. Secondly minification or magnification of K p and T i for Canal System 1 is greater, We speculate, that is because characteristics of water surface themselves play an important part, they distinguish the integral properties.  We took the Canal 1 in Canal System 2 as an example, other canals show similar responses. As Figure 6 shows, all processes feature as an over-damped system and the controller set by smaller R 1 contributes to more control actions. As a consequence, it cuts down stable time and steady-state error, but also causes severer water level oscillation and more frequent gate movements. At a certain range of R 1 , we reckon their control performances are not too bad. Case R 1 = 0.44 corresponds to the calculated value from the method mentioned above. It is noteworthy that this R 1 may not the best but it efficiently avoids deteriorations of any dynamic performance. We may draw to the conclusion that the mentioned R determination could meet our basic control needs and it could be regarded as a mean judgment and weight on water level and flow rate. If we have further needs, for example requiring better flow rate performances, we may could increase or decrease the weight of R based on this weight value.

Conclusions
Based on the ID canal model and PI feedback control model, implanting the LQR optimal theory, a canal control model for discrete systems is established. A set of suitable PI controller parameters can be obtained by adjusting the weight matrices Q and R according to systems' response characteristics. The following conclusions are drawn regarding the application of LQR optimal design: 1) The designed PI feedback controller with linear quadratic optimal theory for gradually sloping canals under constant downstream water-level operation works well. This paper proposed a method to set the R matrix according to canal designed data to avoid cumbersome trial procedure. Simulations of two different canal systems show that this method is feasible for different canals. The controller is able to superiorly balance the demand on water level, flow rate changes and gate movements compared with others. This engineering tuning method of R is practicable for open-channel controller design.
2) In this study, different set of R matrix is compared with typical response processes. LQR controllers at a certain range of R can control the water level deviations within acceptable range and bring the flow rate to the target value thus it is much robust comparing to local PI controllers. Smaller R 1 produces more sensitive controller mainly through enlarging K p meanwhile reducing T i . The diversity of canal response characteristics makes the reasonable range of R 1 different: for narrow and deep channel, it is much larger for its sensitivity to water level and inherent physical parameters.
It should be noted that this controller in this way may not be the best on the whole interval but it could serve as an optimal solution in the controller parameters choosing. Our study suggests a quite promising tuning approach for the PI feedback controller. We plan to combine LQ optimal theory with other feedback controllers in our further study.