Fractional order operators

Fractional calculus involves the generalization of standard integration and differentiation to non-integer (fractional) order fundamental operators. These are represented as _aD^α_t where a and t are the limits and α(α ∈ ℜ) is the order of the operator. Several definitions of this operator have been proposed (Podlubny, 1999), all of which generalize the standard differential/integral operator in two principal respects: a) they become the standard differential/integral operator of any order when α is an integer, b) the Laplace transform of the operator _aD_t^α is s^α (provided there are zero initial conditions), and hence the frequency characteristics of this operator are (jω)^α. This last feature is very appealing as regards its use to identify linear models from their frequency responses, as it allows us to describe systems with new asymptotic behaviors in the frequency domain (both in magnitude and phase).

One of the most frequently used definitions of the general fractional differential operator is that of Riemann-Liouville (the RL definition) (Podlubny, 1999):

where n – 1 < α < n , Γ(.) is an integer positive number, is the Euler’s gamma function and t > a. Another definition is that of Caputo which is often preferred, as the initial conditions that appear in its Laplace transform are integer derivatives of f(t) while the initial conditions of Eq. [1] are given in terms of fractional derivatives of that function.

The Laplace transform of any of these definitions under zero initial conditions for order α, (0< α < 1) is given by:

A fractional order system is that system described by the following fractional order differential equation:

where α_n > α_n–1 > ... > α₁ > α₀ and β_m > β_m–1 > ... > β₁ > β₀ are arbitrary real numbers, and , a_k (k = 0, 1, ..., n), b_k (k = 0, 1, ..., m) are constant coefficients. The fractional order transfer function is therefore given by the following expression:

where P(s) and Q(s) have no common zeros and the normalization a_n = 1, ..., α₀ = 0 has been carried out.

If a fractional system has to be simulated, then it is necessary to approximate expression [1] by means of discrete realizations. These can be obtained in two ways (Vinagre et al., 2000): i) by approximating the fractional operator using a standard transfer function in the frequency range of interest and then applying any habitual discretization technique (e.g., Tustin operator), ii) by numerically approximating the fractional operator. In the second way, the numerical approximations of the fractional derivative/integral operator are often implemented by using the following numerical generalization:

where and T is the period of discretization. This operator is formally expressed as:

and can be implemented by using:

where [.] represents the integer part, and the combinatorial function has been generalized in the following respect:

This numerical approximation is called the Grundwald-Letnikov (GL) definition of the discretized fractional operator (Podlbny, 1999), and provides sufficiently accurate results in most cases, provided that T is sufficiently small. This discrete operator may be approximated using FIR or IIR discrete filters. Expression [8] can be truncated to a fixed number of N + 1 terms of this sum (0 ≤ j ≤ N) in oder to obtain a FIR filter. This can be done because and it is called the short memory approximation. This approximation usually leads to FIR filters of a very large order (often N > 100).

The prototype hydraulic canal

The prototype hydraulic canal used as a case study in the present research is a closed-loop water variable slope rectangular canal with glass walls and a methacrylate bottom, manufactured by the Spanish company MCD-2000, and located in the Fluids Mechanics Laboratory of the Castilla-La Mancha University (Spain). The canal is 5 m long and 8 cm wide, while its walls are 25 cm high. The estimated manning coefficient is 0.001. It consists of a closed-loop water canal with an instrumental platform that integrates electro-mechanical sensors and actuators, a PLC (programmable logic controller) and a SCADA (data acquisition and supervisory system). The canal has motorized and manual adjustable slide gates (with a discharge coefficient of 0.4) that allow it to be divided into pools of different lengths, but considering its small dimensions the canal is fundamentally operated as a single main canal pool of approximately 4.7 m in length and with a downstream end operation method. The upstream gate of this main pool is a motorized undershot gate and the downstream gate is a manually adjustable overshot gate. A view of this prototype hydraulic canal is shown in Fig. 2 and its schematic representation is depicted in Fig. 3. The geometry of the downstream overshot gate is shown on the right-hand side of Fig. 3. The angular position of this gate can be manually adjusted to three values which yield gate top heights of 13, 23 and 33 mm above the canal bottom, respectively.

Figure 2. Prototype hydraulic canal in laboratory.

Figure 3. Scheme of the overall control systems of the prototype hydraulic canal. DAQ, data acquisition card; GPS, gate position sensor; US, ultrasonic level sensor; FPI/PI, fractional proportional integral/proportional integral controllers; PID, proportional integral derivative controller; DC motor, direct current motor.

The water flows in a closed circuit from an upstream reservoir to a downstream storage reservoir in order to economize water. The return of the water to the upstream reservoir is guaranteed by an electric pump. This pump operates with a speed variator, thus allowing adjustments in frequency from 0 to 50 Hz. The total canal inflow is adjustable from 0 to 9 m3/h (≈ 2.5 L/s). The canal does not therefore have any water losses except evaporation effects, which are not significant in this case. Moreover, the laboratory canal has a water extraction system with which to simulate users’ water offtakes. Taking into account the maximum inflow permitted and the minimum downstream water level measured in our experiments (see Fig. 4), the maximum Froude number attained in the downstream of the main pool is 0.97. The regime of the flow at the end of the main pool is therefore always subcritical.

Figure 4. Different step responses of the prototype hydraulic canal.

Two ultrasonic sensors (US), located outside the top of the canal, are used to monitor and control the upstream (y_up (t)) and downstream end (y_dwe (t)) water levels. A third ultrasonic sensor is shown in Fig. 3, and this measures the water level immediately downstream of the undershot gate. When the canal has a variable slope, measurements of the second (y_dw (t)) and third (y_dwe (t)) ultrasonic sensors are different. The motorized slide undershot gate is equipped with a DC motor and a gate position sensor (GPS). The canal is also equipped with two electromagnetic flowmeters (EMF). One flowmeter is installed in the water return tube which measures the canal water inflow (Q_in (t)) pumped by the electric pump from the storage downstream reservoir toward the upstream reservoir, while the other is installed in the water extraction system tube which measures the canal water outflow (Q_out (t)) pumped by the electric pump from the canal laboratory. This is not shown in Fig. 3 because neither the water extraction system nor its outflow pump or its flowmeter are used in our experiments. These two flowmeters facilitate the supervision of the pumping operations.

Our hydraulic canal prototype uses a PC (personal computer) as a canal control station. A SCADA application is installed in this PC to ensure the automatic control and supervision of the canal. Signals from the various installed sensors are measured and recorded for the: upstream (x_up (t)) gate position, upstream (y_up (t)), downstream (y_dw (t)), and downstream end (y_dwe(t)) canal water levels, canal water inflow (Q_in (t)) and canal water outflow (Q_out (t)).

Control of this hydraulic canal is very complicated because the maneuvers needed to open and close the upstream gate produce very large changes in the upstream water levels, owing to the small dimensions of the upstream pool (see Fig. 2). In order to solve this problem, a secondary control loop that controls the upstream water level and maintains it at a fixed reference was implemented. The upstream water level is controlled by a PID controller, which acts on a speed variator (frequency converter) that varies the water flow of a pump. The primary control loop carries out the control of the downstream end water level in the main pool of the hydraulic canal. The aforementioned upstream water level control then decouples the main pool dynamics from the upstream pool dynamics, and allows us to focus on identifying only the main canal dynamics without having to pay attention to secondary dynamics caused by any interaction with the upstream pool.

All the hydraulic canal prototype controllers are installed inside the PLC, which is supervised by the SCADA application. Fig. 3 shows the overall control system for this canal. The primary control loop can be seen in the upper part of this figure, and the secondary standard control loop can be seen in the lower part. The SCADA application developed facilitates the implementation of different control strategies, such as fractional order control, standard, robust and nonlinear controls, along with various canal operation methods (upstream, downstream, downstream end, mixed, Bival, etc.), and set-point changes (y_{dwe_sp} (t)) and (y_{up_sp} (t)) of the primary and secondary control loops, respectively. This SCADA application also provides other facilities such as: 1) storage of all the input and output signals in a data-base, allowing their exportation to other programs, such as MS-Excel, Matlab, etc.; 2) a display of the state of all the input and output signals; 3) alarm generation making it possible to obtain information about verified damages, the suggestion of the corresponding actions and, in extreme conditions, even the making of automatic decisions. The main alarms are: water levels inside the hydraulic canal (in order to prevent the canal from emptying or overtopping), and the operational state of the control system devices (canal control station, DC motor, gate position and ultrasonic sensors, flowmeters, pumps, speed variators, etc.).

Dynamic models of irrigation canals and identification process

The dynamics of water flowing in irrigation canal pools is modeled by using the Saint-Venant equations. These equations are derived from mass and momentum balances and are given by (Chaudhry, 2008):

where A (x,t) is the canal cross section area, Q (x,t) is the water flow (discharge) through section A, q (x,t) is the lateral discharge, y (x,t) is the water depth, S_f (x,t) is the friction slope, S_b is the bed slope, g is the gravitational acceleration, t is the time variable, and x is the longitudinal abscissa in the direction of the flow.

Several methods currently exist for the solution of the Saint-Venant equations, but all of them have considerable mathematical complexities (Litrico & Fromion, 2004). These are thus often solved numerically by using specific software (e.g., Valipour, 2012; Khasraghi et al., 2015). Moreover, the direct use of these equations in the analysis of flow dynamics and the design of controllers is very difficult (Malaterre, 1995). The Saint-Venant equations are often linearized around a set-point. In this case, an equivalent first-order plus a time delay system is used to model the canal pool dynamic behavior (Rivas-Perez, 1990; Weyer, 2001). Experiments developed by some authors have demonstrated that the dynamical parameters of these linearized models may vary considerably when the discharges through the upstream gates of the canal pools vary in the operation range [Q_min, Q_max] (Litrico & Fomion, 2009; Rivas-Perez et al., 2011).

A linear parameter variable (LPV) non-integer order model of a canal pool was obtained using fractional-order LPV identification (Martinez et al., 2009). Rational local fractional-order models were obtained at different operation points, and a global LPV model was then obtained from polynomial interpolations of the parameters of these local models.

Experiments based on the response to a step like input were carried out in our hydraulic canal prototype in order to obtain a mathematical model that would accurately describe its dynamic behavior. In these experiments, the downstream gate was fixed in one of the three possible positions (13, 23 and 33 mm) above the horizontal level of the canal bottom. Different step movements were performed by the upstream gate in order to excite the canal dynamics and the downstream end water level was measured with an ultrasonic sensor. Furthermore, experiments were carried out with different upstream water levels (60, 65 and 70 mm), thus allowing us to obtain suitable and exhaustive experimental data that would permit an adequate characterization of the dynamic behavior of our canal. During the experiments, the upstream gate was in free or submerged flow, which depends on the upstream water level and the position of the downstream gate. The downstream gate was in free flow. The water levels and upstream gate positions were sampled with a period of 0.15 s. The positions of the upstream gate and the water level are given in mm.

Fig. 4 shows different step responses of our hydraulic canal prototype. This figure shows that the upstream reservoir water level is maintained approximately constant by the pump, whose speed variator is controlled by the PID. This control system compensates for the variations in the upstream water level caused by the different openings and outflows of the upstream gate (note that only fast transients that are quickly removed are produced on the upstream water level by the upstream gate opening changes).

According to the responses obtained after carrying out the experiments developed, we have proposed three simple linear models with which to characterize the dynamic behavior of our hydraulic canal prototype.

Standard first order plus time-delay model

This is the simplest and probably most frequently used canal pool model (Malaterre, 1995). It is based on simple physical considerations. In order to obtain an insight into the hydraulic process involved, let us consider a very simplified dynamic model of the main pool, which is labeled as Pool 2. Consider the following equation, which represents the flow balance just upstream of Gate 3:

where, according to Fig. 1, Q₂ (t) and Q₃ (t) are the input and output flows to the main pool, q₂ (t) is the offtake (assumed to be close to Gate 3), A is a constant, and τ is the time-delay that represents the time required by an increment of mass at the upstream gate to move to the position of the downstream gate. As Fig. 3 shows, our canal uses a submerged undershot gate like Gate 2, and an overshot gate like Gate 3 (this gate has a constant opening). These flows are therefore those given by e.g. Chaudhry (2008):

where K₂ and K₃ are constants, u₂ is the Gate 2 opening, and h is the constant height of the top of Gate 3 with regard to the canal bottom; (y_dwe (t) – h) is thus the upstream head over the gate.

In the steady state, the flow balance is where the overline signifies steady state values. The linearization of Eqs. [11] and [12] around this steady state, assuming that q₂ = 0 and taking only the input u₂ and the output y_dwe (t) as the variables representative of the process, leads to the following incremental model:

where the overline again represents initial steady state values of the gates’ opening and water levels. This linearized model is a first order plus time delay system whose parameters change dramatically depending on the operation regime. Its transfer function (taking the Laplace transforms in Eq. [13]) is:

The transfer function [14] can be expressed as:

where K is the static gain, T is the time constant, and τ is the time-delay. Expression [14] shows that these parameters change according to the canal operation regime.

An algorithm is designed in order to estimate parameters K, T and τ of Eq. [15] in each canal operation regime. The quality of the estimated model is measured using the value of the Integral Square Error (ISE), in which the error is defined between the step input response of the estimated mathematical model and the response experimentally obtained from our canal. The ISE is a performance index that is very often used in system identification (Ljung, 1999) and which is defined as:

Even though the ISE is calculated from 0 to infinity in Eq. [15], we calculate the ISE in the range from 0 to 360 s in order to save computation time, as all the experiments reached their steady state in less than 360 s (and the error between the model output and the real response would therefore become approximately zero in less than this time, if the gain K had been properly tuned). All the experiments were executed during the same time interval of 360 s.

In our identification approach, therefore: i) we manually determine the time-delay value τ by means of a simple inspection of the experimental response, and ii) parameters K and T are obtained from an optimization procedure that minimizes the aforementioned ISE index. The nominal values obtained for the laboratory hydraulic canal model (values obtained in the case of the most frequent operation regime) are: K₀ = 0.372, T₀ =10.975 and τ₀ =4.05. They correspond to an upstream reservoir water level of 70 mm, a downstream end water level of 50 mm in the main pool, and a gate opening of 15 mm. However, these parameters have wide ranges of variation around their nominal values when the operation regime changes.

We repeated the identification process of model [15] using the canal step input response data obtained in all the possible operation regimes: we covered the entire flow range, we used different upstream water levels, and we opened the downstream gate in the main pool using different settings. All the possible combinations of the following values were considered:

where u_i and u_f are the initial and final values of a gate manoeuvre (the gate position undergoes a step from u_i to u_f). This yielded a total amount of 1080 different experiments that were carried out. Model [15] was fitted to the response data obtained in each of these experiments. This set of experiments yielded the following ranges of model parameter variations: 0.086 ≤ K ≤ 0.775, 1.1 ≤ T ≤ 48.94 and 2.25 ≤ t ≤ 5.7.

The accuracy of the model [15] estimated in the different operation regimes was sufficient for control purposes, but we wished to explore whether, by just slightly increasing the complexity of the model, a better fit to the experimental data could be achieved, which would obviously have a significant impact on the design of the controller.

We have three parameters (K, T and τ) in the first order plus time delay model. In the following subsections we explored the impact on the model accuracy if a fourth parameter is added to the model. We therefore fit models with four parameters to our already recorded experimental data, and we considered two new model structures, which are natural extensions of the previous one. The first is of an integer-order nature (the second order plus time delay model), while the second is of a fractional-order nature (a fractional-order plus time delay model), which implies using only a fractional order derivative operator s^α in the model [15].

Second order plus time-delay model

This is a standard integer-order model that has also often been used to model main irrigation canal pools (see e.g. Rivas-Perez et al., 2008):

where: K is the static gain, T₁, T₂ are the time constants, and τ is the time-delay.

The identification method applied in this case consisted of the following: i) parameters K and τ were taken from the previous identification process, because the time-delay and the static gain would remain the same in both models [15] and [18], as they had been obtained from the first and the last instants of the response respectively, which were not influenced by the transient response, and ii) parameters T₁ and T₂ were again obtained from an optimization procedure that minimized the ISE performance index using model [18].

We ran this identification process on the data of the 1080 experiments, and yielded the following ranges of model parameter variations: 0.086 ≤ K ≤ 0.775, 1 ≤ T₁ ≤ 48.92, 0 ≤ T₂ ≤ 13.18 and 2.25 ≤ τ ≤ 5.7.

In many of the identification experiments carried out, the second order model with a time-delay fitted converted into a first order model with a time-delay, because the optimum T₂ obtained was zero. This meant that adding a second pole to the model did not, in many cases, improve the accuracy attained with the first order plus time-delay model.

There is another integer-order model with four parameters to be tuned that could be used:

but this model is not suitable for our canal as it implies a discontinuity in the step input response (just when the response starts varying, after the time-delay has passed), and we verified that the experimental responses of our canal pool did not have such discontinuity.

Fractional-order plus time-delay model

The previous subsection allows us to conclude that increasing the complexity of the transfer function obtained from an integer-order differential equation by adding one more parameter does not produce more accurate models of our main canal pool. In this subsection we therefore model our canal pool dynamics with a fractional-order differential equation plus time delay. We propose a simple fractional-order model that also has four parameters to be tuned (a model with the same complexity of the second order plus time-delay model shown in the previous subsection and with only one more parameter to be tuned than in the first order plus time-delay model):

where K is the static gain, T is the time constant, τ is the time-delay, and τ is the fractional order.

The identification method applied in this case was similar to that used in the previous subsection. It consisted of: i) removing parameters K and τ from the first order plus time-delay identification process, for the same reason as that shown in the previous subsection, and ii) once again obtaining parameters T and α from an optimization procedure that minimized the ISE performance index using model [20].

We again ran this identification process on the data of the 1080 experiments, and yielded the following ranges of model parameter variations: 0.086 ≤ K ≤ 0.775, 1.17 ≤ T ≤ 48, 2.25 ≤ τ ≤ 5.7 and 0.59 ≤ α ≤ 1.32.

Table 1. Comparison of the three models with the error index calculated in an interval of 360 s and in a time interval that lasts until the instant at which the experimental response has reached its steady state