Wiener-Hopf optimal control of a hydraulic canal prototype with fractional order dynamics

ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

Contents lists available at ScienceDirect

ISA Transactions journal homepage: www.elsevier.com/locate/isatrans

Research article

Wiener-Hopf optimal control of a hydraulic canal prototype with fractional order dynamics Vicente Feliu-Batlle a,n, Daniel Feliu-Talegón b, Andres San-Millan b, Raúl Rivas-Pérez c a

Universidad de Castilla-La Mancha, Escuela Técnica Superior de Ingenieros Industriales. Av. Camilo Jose Cela s/n, 13071 Ciudad Real, Spain Instituto de Investigaciones Energeticas y Aplicaciones Industriales (INEI), Campus Universitario de Ciudad Real, 13071 Ciudad Real, Spain c Havana Technological University (CUJAE), Department of Automatic Control and Computer Science, Calle 114, No. 11901, Marianao, 19390 La Habana, Cuba b

art ic l e i nf o

a b s t r a c t

Article history: Received 15 September 2016 Received in revised form 7 June 2017 Accepted 16 June 2017

This article addresses the control of a laboratory hydraulic canal prototype that has fractional order dynamics and a time delay. Controlling this prototype is relevant since its dynamics closely resembles the dynamics of real main irrigation canals. Moreover, the dynamics of hydraulic canals vary largely when the operation regime changes since they are strongly nonlinear systems. All this makes difficult to design adequate controllers. The controller proposed in this article looks for a good time response to step commands. The design criterium for this controller is minimizing the integral performance index ISE. Then a new methodology to control fractional order processes with a time delay, based on the Wiener-Hopf control and the Padé approximation of the time delay, is developed. Moreover, in order to improve the robustness of the control system, a gain scheduling fractional order controller is proposed. Experiments show the adequate performance of the proposed controller. & 2017 ISA. Published by Elsevier Ltd. All rights reserved.

Keywords: Wiener-Hopf filter Fractional order control Hydraulic canal control ISE optimization Gain scheduling control

1. Introduction Nowadays, the effective management of limited available water resources is a constant source of research because the water is the most valuable resource in the world [1]. The irrigated agriculture is the largest consumer of freshwater: about 70% of all freshwater withdrawals go to irrigation systems [2,3], especially in arid and semi-arid areas that are mainly characterized by low rainfall and high evaporation [4,5]. The water distribution is carried out in these systems through irrigation main canals, whose mainly objective is to convey water from its source down to its final users [6]. Huge water losses take place in these canals. The research focused on improving the management and efficiency of this class of plants has therefore a high scientific, economic and social interest [7]. The design of effective controllers of water distribution is an efficient way to solve the problems of water management, reduction of the huge losses of water in irrigation main canals, and promotion of sustainable development of irrigation areas [8,9]. However, the design of these controllers is not a simple task, since irrigation main canals are complex systems with distributed parameters over long distances, significant time-delays, strong n

Corresponding author. E-mail address: [email protected] (V. Feliu-Batlle).

nonlinearities and dynamics that change with the operation conditions [10]. Several strategies have been proposed to control irrigation main canals, which are mainly based on heuristic, classical, optimal, predictive, intelligent or robust algorithms [11,12]. Probably, the most widely used solution is the classical PID controller, due to its robustness, accuracy and easiness of implementation in the field [8,10]. Since canal dynamics exhibits time delays, control systems based on the Smith predictor and its modifications [13] have also been proposed. However, many studies have shown that the previous solutions do not perform well in the cases in which the irrigation main canals have dynamic behaviors that vary with the operation regime [14,15]. The design and implementation of controllers for irrigation main canals require mathematical models that accurately describe the dynamic behavior of these plants under realistic conditions [8]. However, the construction of these mathematical models can be a very difficult and laborious task [10]. Different linear models of irrigation main canal pools have been proposed. They are based either on Saint-Venant equations [8,16,17] or on the use of system identification tools [18–20]. A review of the main approaches used in the development of models for irrigation main canal pools can be found in [21]. Some authors proposed a linear time invariant (LTI) model of an irrigation main canal pool based on a first order with a time delay transfer function [8,18,22] of the form

http://dx.doi.org/10.1016/j.isatra.2017.06.012 0019-0578/& 2017 ISA. Published by Elsevier Ltd. All rights reserved.

Please cite this article as: Feliu-Batlle V, et al. Wiener-Hopf optimal control of a hydraulic canal prototype with fractional order dynamics. ISA Transactions (2017), http://dx.doi.org/10.1016/j.isatra.2017.06.012i

V. Feliu-Batlle et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

2

Gi(s ) =

Ki e−τi s 1 + Ts i

(1)

It is quite clear that the complete system dynamics is not modelled adequately with a linearized model and that the dynamical parameters of the model may experience variations when the canal operation regime changes [8]. The degree of adequacy of model (1) for the design of high performance control systems is often not that required owing to the nonlinear behavior of the canal, to diffusive phenomena that are not adequately described and to parameter uncertainties [10]. There are therefore significant open research problems in this field yet. Fractional calculus is a direct expansion of the traditional integer order calculus. It allows the order of the derivative of a function to be a fractional number instead of integer [23]. Researchers have found that fractional order differential equations are more adequate than integer order equations to model certain processes, thus providing an excellent tool with which to describe the dynamics of these processes [24]. In particular, fractional order calculus has shown to be very effective in modeling distributed parameter processes that involve partial differential equations, as occurs in electrochemical, thermal or hydraulic processes [24]. Regarding irrigation canals, a discrete fractional order model that included a time delay was proposed in [25]. This model involved 16 parameters to be tuned. Recently, a simpler continuous fractional order model has been proposed in [26], which involved only 4 parameters to be tuned and yielded a reasonable accuracy. This second model includes a fractional order derivative and a time delay, being of the form

Gf (s ) =

Kf 1 + Tf s α

e−τf

specifications and the resulting time specifications is often unaccurate in the case of plants with time delays. Examples of design of controllers in the frequency domain for FROPTD processes are [39], in which Ziegler–Nichols-type rules have been obtained for tuning fractional order PID controllers, or [40], in which the parameters of fractional order PID controllers have been obtained by solving a min-max optimization problem of a H∞ norm. Very few works have been devoted to design controllers that achieve time specifications like overshoot, settling time, integral performance indexes, etc., in FROPTD processes. We only mention the works carried out to design model predictive controllers for fractional order processes, e.g. [41], which, though do not intend to control processes with time delays, can be easily modified to control FROPTD processes. Fractional order controllers have already been applied to canal automation. But these controllers are designed using integer order models with time delays. Some examples are [9,42] in which the standard control loop was used and [13,22] in which the Smith predictor control scheme was used. An implementation of a fractional order controller in a real irrigation canal is described in [43]. The only work up to date that carries out the design of a fractional order controller based on a FROPTD process is described in [44]. This article proposes a gain scheduling controller with a pole placement for a laboratory prototype of irrigation canal, which is based on the complicated fractional order canal model developed in [25]. This work guarantees closed loop stability, it does not pursue any specific time specification (though it proposes a general way of improving the time response by achieving pole placement), and provides simulated results. The following issues can be stated from the previous review:

s

(2)

This model improved in 30% the accuracy attained by model (1) in describing canal pool dynamics, at the cost of having to tune only one parameter more: the fractional order of the model. This fractional model will be used in this article in the design of the control system. Some authors have shown that the use of fractional order models in the design of controllers yields control systems that are often of fractional order too, and that are more effective than the traditional integer order controllers whose design is based on integer order models [27,28]. In particular, controllers for simple fractional order processes (having transfer functions of low orders and few coefficients to tune) without time delay were developed in [29] and, more recently, in [30,31]. Several research works have also been carried out on the control of simple fractional order systems with time delays. Some relevant results are subsequently mentioned. The stability of LTI fractional order systems with time delays (FROPTD) has been studied in [32], and the case of LTI fractional order systems with a fractional order delay has been studied in [33]. An algorithm for stabilizing FROPTD systems using fractional order PID controllers was proposed in [34] and the set of all the PID controllers that stabiIize the class of FROPTD systems that have an unstable pole was obtained in [35]. Though an important research effort has been done to design stable controllers for FROPTD processes, relatively little effort has been devoted to design controllers that achieve time or frequency specifications in these processes. Controllers that achieve frequency specifications can be designed for FROPTD processes by applying methods that are commonly used to design fractional order controllers that achieve frequency specifications in integer order processes, e.g. [36–38]. Often, these methods combine the achievement of frequency specifications and robustness objectives. However, the correspondence between the achieved frequency

1) Methodologies to design controllers that achieve time specifications in FROPTD processes have to be developed. 2) Design methodologies that yield simple tuning rules for the controllers are required since canal dynamics changes with the operation regime and the controller parameters may need to be changed accordingly in real time, as for example in gain scheduling or other adaptive control techniques. 3) At present, the techniques used to design controllers for FROPTD processes require an “a priori” definition of the controller structure. This structure is commonly the fractional order PID controller. However, it can not be guaranteed that this structure is the best one in many applications. 4) Techniques to design controllers based on FROPTD models can improve the performance of the canal automation systems, since they are based on fractional order models that in some cases describe more accurately than integer order models the canal dynamics. 5) Experimental results of control systems of FROPTD processes have not been reported in the scientific literature, as far as we know. This paper addresses the above issues. It proposes the automation of a laboratory canal prototype using a controller that minimizes the integral of the square error (between the process output and the command signal) performance index (ISE). This index has been extensively used to design controllers for integer order plants (e.g. [45]). The controller that minimizes the ISE is given by the Wiener-Hopf design method of optimal controllers, see e.g. [46]. This method has the following interesting features: (a) the minimized index is a measure of the performance of the time response of the process, (b) the controller design procedure gives a closed form solution which is simple and does not require any iterative procedure, being thus suitable for adaptive control, and (c) the method yields the controller with the optimal structure, i.e, it does not require the “a priori” knowledge of the

Please cite this article as: Feliu-Batlle V, et al. Wiener-Hopf optimal control of a hydraulic canal prototype with fractional order dynamics. ISA Transactions (2017), http://dx.doi.org/10.1016/j.isatra.2017.06.012i

V. Feliu-Batlle et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

controller structure. The Wiener-Hopf design method seems therefore quite well suited to solve the first three mentioned issues. Moreover, this method was extended to control fractional order plants in [47]. However, this method has not been extended to the case of FROPTD processes yet. The main novelties of this article are therefore: (1) to extend the Wiener-Hopf design method to the case of FROPTD processes by using the Padé approximation of the time delay, (2) to design its application to the control of a hydraulic canal and 3) demonstrate its feasibility and performance by carrying out experimentation in a laboratory prototype of a hydraulic canal. This paper is organized as follows. Section 2 describes briefly the laboratory prototype of a main canal pool. Section 3 presents the identification carried out of the process dynamics, and the fractional order model is justified. Section 4 develops a WienerHopf controller for an integer order model of the process. Section 5 proposes a new fractional order controller based on the WienerHopf methodology developed in [47] and shows the advantages over the previous one. Section 6 performs a robustness analysis of the proposed control system taking into account: a) errors produced by the approximation of the real delay by the Padé expression and b) variations in the parameters produced by changes of the canal operation regime. A gain scheduling control system is developed in Section 7. Section 8 presents some experimental results and Section 9 offers some conclusions.

2. Facility description Our experimental prototype of hydraulic canal is a variable slope rectangular canal with glass walls which is located in the Mechanics of Fluids Laboratory of the Castilla-La Mancha University (UCLM), Spain. The water in this canal flows in a closed circuit such that waste of water is avoided. Fig. 1 shows a view of this canal. Fig. 2 depicts a schematic representation of it. Canal dimensions are 5 m long, 0.08 m wide, and 0.25 m walls high. Its upstream slide gate is a motorized undershot gate, while it downstream slide gate is a manually adjustable overshot gate. It allows the division of the canal into pools of different lengths. However, considering its small dimensions, the canal is operated mainly as a single main canal pool with a downstream end operation method. The canal is equipped with an instrumental platform that integrates ultrasonic level sensors (ULS), gate position sensor (GPS), DC motor, flow sensor (FS), electric pumps and a speed variator. The control system includes a SIMATIC programmable logic controller (PLC) S7-300, a control station based in a personal computer (PC), and a SCADA (data acquisition and supervisory system) (see Fig. 2). The water flows from an upstream reservoir to a downstream storage reservoir. An electric pump with a speed variator, that allow adjustments in frequency from 0 to 50 Hz, performs the water return to the upstream reservoir. The total canal inflow is

3

adjustable from 0 to 9 m3/h ( E 2.5 l/s). Three ultrasonic level sensors (ULS), located outside the top of the canal, are used to monitor and control the upstream, downstream and downstream end water levels. The SCADA application has been installed in the control station (PC) and ensures the control and supervision of the canal. Signals are measured from the installed sensors and are recorded. The most important measured variables are: the upstream gate position ( x up(t )), the upstream ( yup (t )), downstream ( ydw (t )), and downstream end ( ydwe (t )) canal water levels, and the canal water inflow ( Q in(t )), which is pumped by the electric pumps from the downstream storage reservoir toward the upstream reservoir. In a real canal with a submerged undershot gate (as is our case), the water inflow to a main pool depends mostly on the gate opening and the difference between the upstream water level and the downstream water level, according to the expression

(

)

Q (t ) = Kq⋅x up 2⋅g⋅ yup − ydw . In a real canal, the changes in yup as consequence of the water discharge of the upstream pool in the main pool are small because of the big volume of water of the upstream pool. However, in our case, the upstream pool is relatively small and the discharge of water on the main pool may produce large variations on yup . Then, the flow Q in provided by the pump aims to maintain approximately constant yup . By doing this, the water inflow to the main pool only depends on the gate opening, which is the usual case in real canals (coupling between adjacent pools may also be taken into account, but its effect is often neglected and, in any case, it can be compensated by several ways, e.g. by a feedforward term). Then, in order to simulate a real canal, the pump flow Q in must vary in order to reduce the changes of yup , which is the magnitude that must be kept approximately constant in the usual operation of a real canal driven by gravity. Then control of this small prototype of hydraulic canal is complicated. Since the dimensions of the upstream pool are small (see Fig. 1), the maneuvers performed by the upstream gate produce large changes in the upstream water level which produces an undesired coupling between this pool (upstream reservoir in Fig. 2) and the main canal pool. A secondary control loop that keeps the upstream water level in a fixed reference has therefore been implemented in order to solve this problem. A PID controller has been designed for this purpose, which acts on a speed variator (frequency converter) that varies the water flow of the pump. The primary control loop carries out the control of the downstream end water level in the main pool of the hydraulic canal. Since the mentioned upstream water level control decouples the main pool dynamics from the upstream pool dynamics, we can focus our research into the identification and control of only the main canal dynamics without having to take care of secondary dynamics caused by any interaction with the upstream pool. The controllers of water distribution in the canal are implemented in the PLC, which is supervised by the SCADA application. Different control strategies and set-point changes can be implemented through the SCADA application. The SCADA application also provides other facilities such as the storage of canal operation signals in a data-base (allowing their exportation to other programs), alarms generation making it possible to obtain information about verified damages, the suggestion of the corresponding actions and, in extreme conditions, even making decisions automatically.

3. Fractional order model of the hydraulic canal prototype

Fig. 1. Hydraulic canal prototype of the UCLM.

The dynamics of our canal are described by the Saint-Venant equations, which are nonlinear hyperbolic partial differential equations (e.g. [6]). As mentioned before, linearized models

Please cite this article as: Feliu-Batlle V, et al. Wiener-Hopf optimal control of a hydraulic canal prototype with fractional order dynamics. ISA Transactions (2017), http://dx.doi.org/10.1016/j.isatra.2017.06.012i

4

V. Feliu-Batlle et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

Fig. 2. Schematic representation of the prototype hydraulic canal.

employed to design canal controllers. The fitting error used in [26] to compare the model accuracy is based on the integral squared error (ISE) between the experimental data and the response of the fitted model. It is of the form:

Error =

Fig. 3. Upstream gate opening manoeuvres and downstream end water level responses.

around some flow regimes are often used in order to design canal controllers. The parameters of these approximated LTI models change depending on the operation regime of the canal. Consequently, experiments based on the responses to step inputs of our hydraulic canal prototype were carried out in order to obtain linear dynamic models that accurately describe its dynamic behavior around several flow regimes. In these experiments, the downstream gate was fixed and the upstream gate performed different step movements in order to excite the canal dynamics. The downstream end water level ( ydwe (t )) is the variable to be controlled. Then this variable is regarded as the output of the process, and is measured by an ultrasonic sensor. Experiments were carried out in [26] with different upstream water levels (60, 65 and 70 mm), different positions of the downstream gate (13, 23 and 33 mm above the horizontal level of the canal bottom) and different opening step manoeuvres of the upstream gate. This yielded 1080 operation regimes. Then 1080 data vectors were recorded. Given a LTI model, the carried out identification procedure yielded different parameters of the model in function of the considered operation regime. Moreover, several LTI models were proposed in [26] to fit the recorded data, and their accuracies in describing the canal pool dynamics were compared. This work demonstrated that a fractional order model of the form (2) significantly outperformed the accuracy achieved by model (1) in describing the canal response, which is the model usually

ISE tf

Δy¯dwe

, ISE =

∫0

tf

2 eiden (t ) dt .

(3)

where is the time of duration of each experiment, Δy¯dwe is the difference between the steady state downstream end water level after the upstream gate opening operation and the steady state downstream end water level before the upstream gate opening e e operation, and eiden(t ) = Δydwe (t ) − Δydwe (t ), where Δydwe (t ) is the response of the fitted model and the symbol Δ means incremental variable with respect to the initial steady state. It was shown in that article that reductions higher than 15% in the index Error were attained if model (2) were used instead of the standard model (1). Moreover, it was shown that other linear models with a complexity similar to (2), i.e. having four parameters to be identified, yielded less accurate models than this one. Fig. 3 shows a sequence of different step responses of our hydraulic canal prototype. The upstream reservoir water level is maintained approximately constant to a reference value of * = 60 mm by the secondary closed loop (pump with the speed yup variator controlled by a PID). The downstream gate was in free flow. The water levels and upstream gate positions are given in mm and were sampled with a period Ts = 0.13 s. The position of the upstream gate is shown in the upper subplot and the water level is shown in the lower subplot. It can be noticed that the responses of the canal change with the operation regime. In this work, we choose a particular operation regime, in which the downstream end water level is ydwe,0 = 54.36 mm. The step input applied to the upstream gate in order to excite the system, allowing the identification of an LTI model of the canal dynamics in this operation regime, is shown in the upper subplot of Fig. 3 surrounded by an ellipse. The obtained response is shown in the lower subplot of Fig. 3, also surrounded by an ellipse. This operation regime will be automated later by closing a control loop. Several controllers will be proposed which will be compared after. Models (1) and (2) were fitted to the data of the considered operation regime, yielding the following identified parameters:

Model ( 1): Ki = 0.273, Ti = 2.34 and τi = 4.9 Model ( 2): Kf = 0.273, Tf = 1.82, α = 0.85 and τf = 4.9.

(4) (5)

Please cite this article as: Feliu-Batlle V, et al. Wiener-Hopf optimal control of a hydraulic canal prototype with fractional order dynamics. ISA Transactions (2017), http://dx.doi.org/10.1016/j.isatra.2017.06.012i

V. Feliu-Batlle et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

5

Then the use of the fractional order models (2) in the design of the control system of our main canal pool is justified and may possibly yield controllers with a better performance than the ones that can be obtained from the integer order model (1). This will be analyzed in the following sections.

4. Controller design methodology for the integer order model As mentioned in the Introduction, we propose to design a controller for the process (1) with parameters (4) that minimizes the ISE index:

ISE =

Fig. 4. Upstream gate opening and downstream end water level responses of the main canal pool and the identified models in the chosen operation regime.

Fig. 4 shows the responses provided by these two models to the step command input produced by a sudden change of the upstream gate from 20 mm to 10 mm (see the upper plot of Fig. 3). It shows that model (2) with parameters (5) reproduces better the real behavior of the canal than model (1) with parameters (4). The fitting error (3) of the transfer function resulting with (4) is 0.0556 and the fitting error of the transfer function resulting with (5) is 0.0461, which implies and improvement in the fitting accuracy of 17% (reduction of the fitting error). Note that, in the calculation of error (3), is the time elapsed from the application of the step input to the instant at which the experimental response reaches its steady state, which is approximately 60 s. Table 1 shows the models (2) fitted to the experimental responses obtained with the twelve different steps of the sequence of Fig. 3. The first column enumerates the transfer functions identified from the step responses (steps are ordered from left to right in Fig. 3 and the step n° 4 is the one used to obtain parameters (4) and (5)). The second column shows the initial downstream water level ydwe, i and the third one the final downstream water level ydwe, f after the manoeuvre (both are steady state values). The last column of this table shows the reduction in the Error index (3) produced by fitting models (2) instead of fitting models (1)*. This table shows strong variations in the gain and the time constant in function of the operation regime, while the time delay changes only slightly. Table 1 Parameter values obtained in the experimental identification of Gf (s ) for all the main canal pool responses to the step sequence of Fig. 3.

∫0

∞

* (t ) − Δy (t ))2 dt, (Δydwe dwe

which is function of the difference between the process output * (t ). Other integral perforΔydwe (t ) and the reference signal Δydwe mance indexes like the IAE, ITSE or ITAE could also be minimized in order to obtain optimal controllers. However, the ISE is an index especially adequate to our problem because it weights more large errors and tends to reduce the influence of small errors. Fig. 4 shows that the steady state response of the main canal pool has some noise caused by the low accuracy of the sensor (measurements have a resolution larger than 0.05 mm) and some small waves on the water surface (high frequency dynamics that have not been considered in the modelling). Since these effects are usually neglected in canal automation, we are interested in minimizing an index like the ISE, that presents low sensitivity to these small steady state errors. Note that the other mentioned integral indexes are much more sensitive to steady state errors than the proposed one. Then designing controllers that minimize them may yield controllers less efficient in removing the large errors produced during the transient response. Moreover, there are not closed form solutions to obtain the controllers that minimize these other integral performance indexes, and their determination requires iterative procedures, unlike in the process of obtaining the controller that minimizes the ISE index. * (t ) is generated by We assume that the reference signal Δydwe passing a step command Δy*(t ) through an integer order prefilter of the form

F (s ) =

b⋅s + a , s+a

ydwe, i

ydwe, f

τf

Kf

Tf

α

ΔError (%)

1 2 3 4 5 6 7 8 9 10 11 12

57.4 56.86 56.08 54.36 51.64 54.18 56.03 57.05 57.66 56.08 51.67 56.00

56.86 56.08 54.36 51.64 54.18 56.03 57.05 57.66 56.08 51.67 56.00 57.65

4.8 5 5 4.9 5.2 4.8 5 5 5.1 5.2 5.2 5.2

0.055 0.078 0.172 0.273 0.255 0.185 0.102 0.061 0.079 0.221 0.216 0.083

1.3 1.95 1.82 1.82 1.43 2.86 2.21 1.43 1.56 2.08 1.69 1.56

1 0.9 0.85 0.85 0.85 0.85 0.95 1 0.9 0.9 0.8 0.9

0 3 12 17 15 9 1 0 5 17 23 7

n It is stressed that it is only a small set of the 1080 identified models obtained in different operation regimes.

0
0≤b≤1

(7)

Since b ≤ 1, this prefilter is a phase lag compensator. Moreover, this prefilter has unity gain ( F (0) = 1) and represents different behaviours in function of the coefficient b : if b = 0, a first order prefilter F (s ) = a/(s + a) is obtained, and if b = 1, a unity gain prefilter F (s ) = 1 is obtained. The Laplace transform of the resulting reference signal is therefore:

* (s ) = ΔY dwe

Number of the step

(6)

F (s ) , s

(8)

which becomes a step command if b = 1. Fig. 5 shows several reference signals obtained with this filter. Moreover, assuming that a settling time ts is required for this signal, the following relationship is easily obtained:

a=−

⎛ 0.05 ⎞ 1 ⎟ ⋅ ln⎜ ⎝1 − b⎠ ts

,

(9)

which is true for values b ≤ 0.95 since the band error associated to the settling time is ±0.05. 4.1. Wiener-Hopf controller Given a rational transfer function G(s ) and the standard control loop, the optimal controller that minimizes the ISE index (6) is

Please cite this article as: Feliu-Batlle V, et al. Wiener-Hopf optimal control of a hydraulic canal prototype with fractional order dynamics. ISA Transactions (2017), http://dx.doi.org/10.1016/j.isatra.2017.06.012i

V. Feliu-Batlle et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

6

Fig. 6. Wiener-Hopf filter scheme.

Moreover, it is verified that if b = 1:

1 + Ti⋅s ψn(τi⋅s ) ⋅ , Ki s⋅χi, n (s )

Ci, n(s ) =

being χi, n (s ) =

ψn(τ⋅s )

ψn(x) =

⎛ n ⎞ ( 2n − k ) ! k x ( 2n ) ! k=0

∑ ⎜⎝ k ⎟⎠

(10)

Denote Gi, n(s ) as the transfer function (1) whose exponential term is approximated by Pn(τ⋅s ):

Gi, n(s ) =

being χi, n (s ) =

⎛⎛ ⎞ ψ ( −τi ⋅ a) ⎞ ⎟⎟ ⋅ s + a⎟⎟ ⋅ ψn(−τi ⋅ s ) a ⋅ ψn(τi ⋅ s ) − ⎜⎜ ⎜⎜ 1 − n ψn( τi ⋅ a) ⎠ ⎝⎝ ⎠ s ⋅ ( s + a)

(16)

a polynomial.

Proof: Given the relationships in Fig. 6, the Wiener-Hopf spectral-factorization solution of the problem in closed form is [46]:

⎡ ΔY * (s )⋅ΔY * ( − s )⋅G( − s ) ⎤ 1 dwe dwe ⎢ ⎥ * (s )− ⎣ * (s )+ G(s )−⋅ΔY dwe G(s )+⋅ΔY dwe ⎦−

W (s ) =

(17)

where X (s ) groups all the stable factors of the numerator and denominator of X (s )⋅X ( − s ) and X (s )+ groups the remaining of X (s )−⋅X (s )+ = X (s )⋅X ( − s ). Moreover, X (s )⋅X ( − s ), so that ⎡ ΔYdwe ⎤ * (s ) ⋅ ΔY * (−s ) ⋅ G(−s ) dwe ⎢ ⎥ groups all the terms of the partial fraction G(s )+ ⋅ ΔY * (s )+ ⎣ ⎦− dwe expansion of

Ki Pn(τi⋅s ) 1 + Ts i

)

−

n

,

a polynomial, and if b = 0:

(

obtained by applying the Wiener-Hopf design method of optimal controllers, see e.g. [46]. Consider transfer function (1), which is irrational because it has a time delay term. The Wiener-Hopf method can not be therefore applied. One way to overcome this problem is to substitute in (1) the time delay by the nth-order Padé approximation (we consider only Padé approximations in which the degree of the numerator and denominator are equal):

ψn( − τ⋅s )

s

⎛ ⎞ ψ ( −τi ⋅ a) ⎜ 1 − n ⋅s + a⎟⋅ψn(τi⋅s ) ψn( τi ⋅ a) ⎠ 1 + Ti⋅s ⎝ Ci, n(s ) = , ⋅ Ki s⋅( s + a)⋅χi, n (s )

* (t ) . Fig. 5. Reference signals Δydwe

e−τs ≈ Pn(τ⋅s ) =

(15)

ψn(τi ⋅ s ) − ψn(−τi ⋅ s )

* ( s ) ⋅ ΔY * ( −s ) ⋅ G(−s ) ΔYdwe dwe G( s )+ ⋅ ΔY *

dwe

(11)

where the coefficients of the transfer function Pn(τ⋅s ) are given in Table 2 for some values of n. The next theorem gives the optimal controller for this process. Theorem 1. Assume a process Gi, n(s ) given by (10) and (11) and a * (t ) given by (7) and (8). Then the controller command signal Δydwe Ci, n(s ) that minimizes index (6) is given by.

(s )+

whose poles are in the complex

left-half plane, i.e., the pole in the origin and the pole −a of (7). * (s ) by (7) Consider that G(s ) is given by (10) and (11) and ΔY dwe and (8). Taking into account that ψn(τi⋅s ) is strictly Hurwitz [48], and all the roots of ψn( − τi⋅s ) lie therefore in the complex right-half K ⋅ψ τ ⋅s plane, it is obtained that G( − s) = − T ⋅ s − 1i / Tn( ⋅ iψ )−τ ⋅ s , G(s)− = s + 11 / T , i i ( i) n( i ) G(s )+ = s+a/b s ⋅ ( s + a)

−K 2 i , T 2 ⋅ ( s − 1 / Ti) i

* (s )⋅ΔYdwe * ( − s) = ΔYdwe

* (s )+ = − and ΔY dwe

b2 ⋅ ( s − a / b) s ⋅ ( s − a)

−b2 ⋅ ( s + a / b) ⋅ ( s − a / b) , s2 ⋅ ( s + a ) ⋅ ( s − a )

* (s )−= ΔYdwe

. Substituting these expressions

in (17) and operating yields

1 + Ti⋅s ( μ⋅s + a/b)⋅ψn(τi⋅s ) Ci, n(s ) = ⋅ Ki s⋅( s + a)⋅χi, n (s )

(12)

−

if 0 < b < 1, being

χi, n

s⋅( s + 1/Ti )⋅( s + a) ⎡ Ti ( s + a/b)⋅ψn(τi⋅s ) ⎤ ⎥ ⋅⎢ s + a/b ⎣⎢ Ki s⋅( s + a)⋅ψn( − τi⋅s ) ⎥⎦

Wi, n(s ) =

( s + a/b)⋅ψn(τi⋅s) − ( μ⋅s + a/b)⋅ψn( − τi⋅s) (s ) = s⋅( s + a)

(13)

a polynomial, and

⎞ ψ ( −τi⋅a) 1 ⎛1 μ= − ⎜ − 1⎟⋅ n ⎝ ⎠ ψn( τi⋅a) b b

(14)

Table 2 Coefficients of polynomials ψn(x ) of the Padé approximation. n

ψn(x )

1

1+

2

1+

3

1+

4

1+

1 ⋅x 2 1 ⋅x 2 1 ⋅x 2 1 ⋅x 2

(18)

Since all the roots of ψn(τi⋅s ) lie in the complex half-plane, the roots of ψn( − τi⋅s ) lie in the complex right-half plane, and it is obtained that

⎡ ( s + a/b)⋅ψ (τ ⋅s ) ⎤ 1/b n i ⎢ ⎥ = + ⎢⎣ s⋅( s + a)⋅ψn( − τi⋅s ) ⎥⎦ s −

( 1 − )⋅ 1 b

ψn( −τi ⋅ a) ψn( τi ⋅ a)

s+a

=

μ⋅s + a/b s⋅( s + a)

(19)

being μ of the form (14). Substituting (19) in (18) gives that

1 ( 1 + Ti⋅s )⋅( μ⋅s + a/b) ⋅ Ki s + a/b

Wi, n(s ) =

(20)

Subsequently, controller Ci, n(s ) is obtained substituting (20) in the expression: +

+

+

1 2 ⋅x 12 1 2 ⋅x 10 3 2 ⋅x 28

Ci, n(s ) = +

+

1 ⋅x 3 120 1 3 1 ⋅x + ⋅x4 84 1680

Wi, n(s ) 1 − Wi, n(s )⋅Gi, n(s ) ⋅

=

1 + Ti⋅s Ki

( μ⋅s + a/b)⋅ψn(τi⋅s) ( s + a/b)⋅ψn(τi⋅s) − ( μ⋅s + a/b)⋅ψn( − τi⋅s)

(21)

Please cite this article as: Feliu-Batlle V, et al. Wiener-Hopf optimal control of a hydraulic canal prototype with fractional order dynamics. ISA Transactions (2017), http://dx.doi.org/10.1016/j.isatra.2017.06.012i

V. Feliu-Batlle et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

It can be easily checked that s = 0 is a root of the denominator of (21). Taking into account that μ is of the form (14), it is also easily obtained that s = − a is another root of the denominator of (21). Then controller (21) can be expressed as (12) with (13). Controller (15) is easily obtained substituting b = 1 in expression (12) and (13) and simplifying the factor ( s + a). Controller (16) is obtained from (15) by multiplying its numerator and denominator by b and calculating limits if b → 0 (note that in this case lim b⋅μ = 1 − ψn( −τi⋅a) /ψn( τi⋅a)). □ b→0

Lemma 1. If 0 < b < 1, n > 0 and a > 0, then it is verified that μ > 1. Proof: Rearrange μ as

⎛ ψ ( −τi⋅a) ⎞ 1 ψ ( −τi⋅a) ⎟⎟⋅ + n μ = ⎜⎜ 1 − n ψn( τi⋅a) ⎠ b ψn( τi⋅a) ⎝ If n > 0 and a > 0, it ψn( −τi⋅a) /ψn( τi⋅a) < 1, and it

1−

ψn( −τi ⋅ a) ψn( τi ⋅ a)

(22) can be easily checked is therefore verified

that that

> 0. Taking into account that: 1) function (22) is

strictly increasing as a consequence of the previous statement, 2) we have that μ = 1 if b = 1 and 3) 1/b > 1 in the range 0 < b < 1, it is yielded that μ > 1 in the considered interval of b, and the lemma is proven. □ Remark 1. Optimal controllers provided by Theorem 1 are always improper. Denoting as nc the relative degree of controllers Ci, n(s ), it is obtained that: 1) In the case that 0 < b < 1, since μ > 1, the relative degree is nc = − 1. 2) In the case that b = 0, since it is always verified that

1−

ψn( −τi ⋅ a) ψn( τi ⋅ a)

> 0, the relative degree is nc = − 1.

3) In the case that b = 1: nc = − 1 if n is odd and nc = − 2 if n is even. Subsequently, a theorem about the stability of the closed loop systems obtained with controllers Ci, n(s ) in the case of using the process with the delay (1) is proposed. Theorem 2. Let us denote the open loop transfer function of the process (1) compensated with the controller Ci, n(s ), given by Theorem 1, as:

L i, n(s ) = Gi(s )⋅Ci, n(s ) =

( μ⋅s + a/b)⋅ψn( τi⋅s)⋅e−τi⋅ s ( s + a/b)⋅ψn(τi⋅s) − ( μ⋅s + a/b)⋅ψn( − τi⋅s)

(23)

Then it is verified that: 1) L i, n(s ) has always a pole in the origin. 2) L i, n(s ) has a pole in −a if b < 1. 3) In the case that b = 0, the best that can be achieved is a marginally stable closed loop system. 4) In the case that 0 < b ≤ 1, the closed loop system is unstable for even values of n and can only be stabilized using odd values of n. Proof: Statements 1) and 2) are trivially obtained from the form of the denominator of the controllers Ci, n(s ) given in Theorem 1. In the case that b = 0 we have that

L i, n(s ) =

⎛ ⎜ 1 − ⎝

(

a⋅ψn( τi⋅s ) −

((

ψn( −τi ⋅ a) ψn( τi ⋅ a)

1−

⎞ ⋅s + a⎟⋅ψn(τi⋅s ) ⎠

)

ψn( −τi ⋅ a) ψ n ( τi ⋅ a)

)⋅s + a)⋅ψ ( −τ ⋅s) n

i

−τi ⋅ s

e

7

( − 1)n + 1⋅e−j ⋅ ω ⋅ τi . Then this Nyquist plot crosses the point ( − 1, 0) an infinite number of times yielding unstable or marginally stable closed loop systems. In the case that 0 < b ≤ 1, the behaviour of the frequency response of (23) at high frequencies is

limω→∞L i, n(j⋅ω) =

μ 1 − μ⋅( − 1)n

(25)

If n were odd, expression (25) would be limω →∞L i, n(j⋅ω) =

μ . 1+μ

According to Lemma 1, μ > 1 if b < 1. Moreover, it is easy to check that μ = 1 if b = 1. In both cases limω→∞( j⋅ω) < 1, which means that the point ( − 1, 0) would not be surrounded an infinite number of times by the Nyquist plot. μ If n were even, expression (25) would be limω →∞L i, n(j⋅ω) = 1 − μ . Taking into account again that μ ≥ 1 if b ≤ 1, it is obtained that

limω→∞( j⋅ω) > 1

( lim

( j⋅ω)

ω→∞

= ∞ if b = 1), which means that

the point ( − 1, 0) would be surrounded an infinite number of times by the Nyquist plot, the closed loop system being then unstable. □ 4.2. Modified Wiener-Hopf controller The controllers obtained from Theorem 1 have several drawbacks that have to be solved prior to their implementation in our canal: 1) According to Remark 1, controllers Ci, n(s ) are improper and cannot then be implemented. 2) According to Theorem 2, even values of n cannot be used in the Pade approximation because they yield unstable closed loop systems when implemented with the real process (1), though they yield stable closed loop systems when implemented with the approximated process (11). 3) It is desired that the time response Δydwe (t ) have a very small overshoot in order to the water level do not surpass the canal walls and do not waste water. This subsection develops methods to overcome these drawbacks. The following lemma proposes a modification of the standard control scheme in order to overcome the first drawback. Lemma 2. Consider a process G(s ) of relative degree ng > 0, a controller C (s ) of relative degree nc (it could be negative), and a prefilter F (s ) of relative degree nf > 0, which can be decomposed in two prefilters F1(s ) and F2(s ) such that F (s ) = F1(s )⋅F2(s ), being the denominators of F1(s ) and F2(s ) strictly Hurwitz polynomials. Let the relative degree of these last two prefilters be n f > 0 and n f > 0 re2 1 spectively (it is verified that nf = n f + n f ). Then the two control 1 2 schemes shown in Fig. 7 have the same transfer function.

H (s ) = F (s )⋅

G(s )⋅C (s ) 1 + G(s )⋅C (s )

(26)

between the reference Δy*(t ) and the output Δydwe (t ). Moreover, the transfer functions involved in the scheme (b) of this figure are proper if it is verified that: (1) ng ≥ n f and (2) nc + n f ≥ 0. 1 1 Proof: It is obvious that scheme (a) has the transfer function (26). Operating in the second scheme it is obtained that

U (s ) = C (s )F1(s )⎡⎣ F2(s )⋅ΔY *(s ) − F1−1(s )⋅ΔYdwe(s )⎤⎦ ,

ΔYdwe(s ) = G(s )U (s )

(24)

whose frequency response at high frequencies is limω→∞L i, n (j⋅ω)=

(27)

and substituting the first equation into the second and operating yields (26) too. The stated conditions on the relative degrees

Please cite this article as: Feliu-Batlle V, et al. Wiener-Hopf optimal control of a hydraulic canal prototype with fractional order dynamics. ISA Transactions (2017), http://dx.doi.org/10.1016/j.isatra.2017.06.012i

V. Feliu-Batlle et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

8

(1) turns to be unstable. Then controller (29) will be subsequently used. 4.3. Controller tuning Fig. 4 shows that the settling time of the experimental response to a step input is about 35 s. Controllers used in irrigation main canals are commonly designed to exhibit settling times between 10 and 20 times the value of the time delay τi , e.g. [8], which are bigger than the canal open loop settling time. In particular, a settling time ts of 58 s (about 12 times τi ) has been chosen in this work. Moreover, a value b = 0.5 is chosen, because of robustness issues that will be analyzed in a subsequent section. Expression (9) yields that a = 0.04 . Then the prefilter for the reference signal is

F2(s ) =

0.5⋅s + 0.04 s + 0.04

(30)

The controller (29) obtained for the nominal process with parameters (4) is: Fig. 7. Control schemes: (a) the standard scheme and (b) the proposed modified scheme.

of the transfer functions involved are easily obtained from the expressions of the blocks that appear in the diagram of Fig. 7(b). □ Remark 2. The two schemes of Fig. 7 are only equivalent from the point of view of the relation between ΔY *(s ) and ΔYdwe(s ). For example, the transfer functions between ΔYdwe(s ) and a disturbance at the input of the process G(s ) are G(s ) /(1 + G(s )⋅C (s )) in scheme (a) and G(s )⋅ 1 + G(s )⋅C (s )⋅( 1 − F1(s )) /(1 + G(s )⋅C (s )) in scheme (b).

(

)

The previous lemma is applied to our problem with a prefilter F (s ) of the form:

F (s ) = F1(s )⋅F2(s )

(28)

in which F2(s ) is the prefilter (7) and F1(s ) = 1/( 1 + ε⋅s), where ε is chosen small enough to do not modify noticeably the optimized ISE, while allowing to solve the first drawback. Since the relative degree of Gi, n(s ) is ng = 1 ∀ n, Condition 1 of the lemma states that n f = 1 and, consequently, Condition 2 states 1

that nc ≥ − 1. A prefilter of the form proposed for F1(s ) is suitable for this. In this case, the feedthrough block of Fig. 7(b) is ε⋅s⋅G(s ). The second drawback can be avoided by making n odd and b > 0. The third drawback is avoided by choosing an appropriate prefilter. The closed loop time responses of process (1) with controllers Ci, n(s ) designed from Padé approximations n = 1 and n = 3 with a step reference ( F2(s ) = 1) have overshoots higher than 25 % , which is considered unacceptable. It is proposed to use prefilters (7) with b < 1, that provide command signals smoother than a step (see Fig. 5) and, therefore, responses with lower overshoots. The previous constraints suggest to use n odd and a prefilter (7) with 0 < b < 1. Application of Theorem 1 with prefilter (7) to process (11) approximated with n = 1 gives:

Ci,1(s ) = μ=

(

τ

)(

( 1 + Ti⋅s)⋅ 1 + 2i ⋅s ⋅ μ⋅s + 2 ⋅ s⋅(s + a) Ki⋅τi⋅( μ + 1) 2⋅b − a⋅b⋅τi + 2⋅a⋅τi 2⋅b + a⋅b⋅τi

a b

)

Ci,1(s ) =

4.637⋅s 3 + 4.189⋅s 2 + 1.072⋅s + 0.055 s⋅(s + 0.04)

(31)

Making F1(s)¼ 1/(1 þ 0.1  s), the transfer functions to be implemented are obtained:

Ci,1(s )⋅F1(s ) = ε⋅s⋅ Gi,1(s ) =

4.637⋅s 3 + 4.189⋅s 2 + 1.072⋅s + 0.055 , s⋅(s + 0.04)⋅( 1 + 0.1⋅s ) −0.0669⋅s 2 + 0.0273⋅s 5.773⋅s 2 + 4.79⋅s + 1

(32)

Fig. 8 shows the simulated closed loop responses to a step command Δy*(t ) passed through the prefilter F (s ) = F1(s )⋅F2(s ) where F2(s ) is given by (30) and F1(s ) = 1/( 1 + 0.1⋅s) - of the process (1) with parameters (4) using the Padé approximation of the delay with n = 1, and control system (32). This figure shows that the closed loop response of the process (1) with its Padé approximation reaches the command signal in about 2⋅τi s. Note that the simulated responses are hereafter normalized assuming that the process has zero initial state and a unity step command Δy*(t ) is applied. The ISE index using the control system (32) and a prefilter F (s ) (with a F2(s ) of the form (30) and a F1(s ) with ε = 0.1) in the case of simulating the model (11) with Padé approximation n = 1 is 1.48. The ISE index obtained using only F2(s ) (assuming F1(s ) = 1) in the closed loop simulation of the model (11) is 1.46. It means that the increase in the ISE caused by adding F1(s ) is only 0.02, which is negligible.

,

(29)

The controller obtained using Padé approximation n = 3 is unstable (it has two poles in the right half-plane). Then, though the closed loop control system using process (11) with n = 3 is stable, the closed loop control system using the delayed process

Fig. 8. Simulated closed loop downstream end water level responses of the process with a Padé approximation of n = 1 and control system (32) in the cases: (a) model (11) using parameters (4) and (b) model (33) using parameters (5).

Please cite this article as: Feliu-Batlle V, et al. Wiener-Hopf optimal control of a hydraulic canal prototype with fractional order dynamics. ISA Transactions (2017), http://dx.doi.org/10.1016/j.isatra.2017.06.012i

V. Feliu-Batlle et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

The gain crossover frequency of the open loop transfer function, Gi,1(s )⋅Ci,1(s ), is ωc = 0.26 rad/s. Since it is verified that

Gi,1(j⋅ω) = Gi(j⋅ω) , the open loop transfer function using the real process Gi(s ) instead of the approximation Gi,1(s ) would have the same gain crossover frequency. However, the phase margin of Gi,1(s )⋅Ci,1(s ) is 115∘ while the phase margin of Gi(s )⋅Ci,1(s ) is 107∘ . Then the error in the phase margin (phase margin is a good index of the relative stability and damping of the closed loop system) committed by using the first order Padé approximation instead of the model with the real delay is about 7% . This error is small enough to let us assume that the controllers designed using the Padé approximation with n = 1 would reasonably work with the real delayed process. Moreover, this assumption will be supported by the results obtained in the experimental section. It is mentioned that a comparison of the simulated responses of the closed loop system with the real delay and with the Padé approximation yield differences only in the first 2⋅τi s: basically, the response of the system with the time delay is zero in the first τi s while the response with the Padé approximation is not.

5. Controller design using the fractional order model

Kf 1 + Tf s α

P1(τf ⋅s )

(33)

Assume that the scheme of Fig. 7(b) were used for this process with the control system (32). Fig. 8 shows the closed loop response of this system to a command (8) with a prefilter F (s ) = F1(s )⋅F2(s ) that is given by (30) and ε = 0.1. This response is worse than the one obtained in the previous section using Gi,1(s ): it is slower and has an ISE index of 4.03, which is significantly higher than the one obtained in the case of the integer order model. This motivates the redesign of the controller using the FROPTD model (2). 5.1. Wiener-Hopf controller The design of controllers using the Wiener-Hopf method was extended to fractional order processes in [47]. In order to carry out a further extension of the method that allows the design of controllers for FROPTD processes, the time delay is substituted again by the Padé approximation:

Gf , n(s ) =

Kf 1 + Tf s α

Pn(τf ⋅s )

( ) ( )

⎞ ψn −τf ⋅a 1 ⎛1 − ⎜ − 1⎟⋅ ⎠ ψ τ ⋅a b ⎝b n f

(34)

Cf , n(s ) =

Cf , n(s ) =

1 + Tf ⋅s α Kf

⋅

( μ⋅s + a/b)⋅ψn(τf ⋅s) s⋅( s + a)⋅χf , n (s )

(35)

( s + a/b)⋅ψn(τf ⋅s) − ( μ⋅s + a/b)⋅ψn( − τf ⋅s)

a polynomial, and

s⋅( s + a)

ψn(τf ⋅ s ) − ψn(−τf ⋅ s ) s

(38) a polynomial, and if b = 0:

⎛⎛ ⎞ ψ ( −τf ⋅ a) ⎞ ⎜⎜1 − n ⎟⋅s + a⎟⋅ψn(τf ⋅s ) ψn( τf ⋅ a) ⎠ 1 + Tf ⋅s α ⎝ ⎝ ⎠ Cf , n(s ) = ⋅ , Kf s⋅( s + a)⋅χf , n (s )

being χf , n (s ) =

) ⎞⎟ ⋅ s + a⎞⎟ ⋅ ψ (−τ ⋅ s) ) ⎟⎠ ⎟⎟⎠ n f

⎛⎛ ψn −τf ⋅ a ⎜⎜ a ⋅ ψn(τf ⋅ s ) − ⎜ ⎜ 1 − ⎜ ψn τf ⋅ a ⎝⎝

( (

s ⋅ ( s + a)

(39)

a polynomial.

Proof: It is omitted because it is very similar to the one given for Theorem 1. □ Lemma 1 applies also to this case. Moreover, since the open loop transfer function of the compensated system is

( μ⋅s + a/b)⋅ψn( τf ⋅s)⋅e−τf ⋅ s ( s + a/b)⋅ψn(τf ⋅s) − ( μ⋅s + a/b)⋅ψn( − τf ⋅s)

(40)

which has the same form as L i, n(s ) in (23), the results of Theorem 2 also apply to this case. Remark 3. Optimal controllers provided by Theorem 3 are always improper. Denoting as nc the relative degree of controllers Cf , n(s ), it is obtained that: 1) In the case that 0 < b < 1, since μ > 1, the relative degree is nc = − α . 2) In the case that b = 0, since it is always verified that

1−

( (

) )

ψn −τf ⋅ a ψn τf ⋅ a

> 0, the relative degree is nc = − α .

3) In the case that b = 1: nc = − α if n is odd and nc = − 1 − α if n is even. Lemma 2 has to be applied again in order to implement controller (35). Its extension to fractional order processes is immediate. A prefilter F (s ) = F1(s )⋅F2(s ) is used in which F2(s ) is of the form (7). Since the relative degree of . is ng = α ≤ 1 ∀ n, Condition 1 of the lemma states that n f = α and, consequently, Condition 1

2 states that nc ≥ − α . These facts imply that prefilter F1(s ) can be of the form and F1(s ) = 1/( 1 + ε⋅s α ) being a > 0, and the feedthrough block of Fig. 7(b) is therefore ε⋅s α⋅Gf , n(s ). Application of Theorem 3 using a Padé approximation of order n = 1 yields

Cf ,1(s ) = μ=

τf

)(

(

)(

1 + Tf ⋅s α ⋅ 1 + 2 ⋅s ⋅ μ⋅s + 2 ⋅ s⋅(s + a) Kf ⋅τf ⋅( μ + 1)

a b

)

,

2⋅b − a⋅b⋅τf + 2⋅a⋅τf 2⋅b + a⋅b⋅τf

(41)

5.2. Controller tuning

if 0 < b < 1, being

χf , n (s ) =

1 + Tf ⋅s α ψn(τf ⋅s ) , ⋅ Kf s⋅χf , n (s )

being χf , n (s ) =

and the following theorem is proposed. Theorem 3. Assume a process Gf , n(s ) given by (34) and a command * (t ) given by (8) with a prefilter (7). Then the controller signal Δydwe Cf , n(s ) that minimizes index (6) is given by.

(37)

Moreover, it is verified that if b = 1:

L f , n(s ) = Gf (s )⋅Cf , n(s ) =

Assume that the previous simulation were carried out using the fractional order model (2) with parameters (5) and a Padé approximation n = 1:

Gf ,1(s ) =

μ=

9

(36)

The same reference signal as in the integer order case is used to minimize the ISE. Then the prefilter (30) is used and the controller (41) is obtained for the nominal process with parameters (5):

( 1 + 1.82⋅s )⋅( 1.982⋅s 0.85

Cf ,1(s ) =

2

+ 0.943⋅s + 0.055

s⋅(s + 0.04)

) (42)

Please cite this article as: Feliu-Batlle V, et al. Wiener-Hopf optimal control of a hydraulic canal prototype with fractional order dynamics. ISA Transactions (2017), http://dx.doi.org/10.1016/j.isatra.2017.06.012i

V. Feliu-Batlle et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

10

Using the F1(s ) proposed in the previous subsection with ε = 0.1, the transfer functions to be implemented are obtained:

( 1 + 1.82⋅s )⋅( 1.982⋅s + 0.943⋅s + 0.055) , (s )⋅F (s ) = s⋅(s + 0.038)⋅( 1 + 0.1⋅s ) 0.85

Cf ,1

2

1

condition to prevent positive roots in (47) is imposing that the coefficient of the term x 2 be positive. This condition is expressed in (44). Moreover, it is easy to check that

0.85

d Lf ,1

(43)

Since the same prefilter F2(s ) is used in the control of (1) and (2), and the same delay is considered in both models ( τi = τf ), the open loop transfer function obtained in this section, Gf ,1(s )⋅Cf ,1(s ), and the one obtained in the previous section, Gi,1(s )⋅Ci,1(s ), are equal. The only difference between the closed loop responses to step commands ( Δy*(t ) in Fig. 7) of both systems is caused by the use of different prefilters F1(s ). The ISE index is 1.48, which is equal to the one obtained in the previous section. This illustrates the little influence of having changed prefilter F1(s ).

Lemma 3. Consider the process Gf (s ) with model (2) and a controller Cf ,1(s ) given by (41). Then the magnitude of the frequency response

open loop transfer function Lf ,1(s ) = Gf (s )⋅Cf ,1(s ) would be a decreasing function that varies from ∠Lf ,1( j⋅ω)

ω=∞

ω=0

=−

π 2

rad to ∠Lf ,1( j⋅ω)ω=∞=

−∞ if it were verified that the equation.

⎛ ⎛ ⎞ ⎞ 1 ⎟ 2⋅a ⎛ 1 ⎞⎟ a3 x2 + ⎜ a2⋅⎜⎜ 1 + + 1 − ⋅x + ⎜ ⎟ 2⎟ ⎜ ⎟ μ⋅b ⎠ μ⋅b ( μ⋅b) ⎠ τf ⎝ ⎝ ⎝ ⎠

that

varies

from

Lf ,1(j⋅ω)

ω= 0

=∞

(48)

has no real positive roots. Proof: The phase of the open loop transfer function (45) is ⎛ ω ⎞ ⎛ μ⋅b⋅ω ⎞ ⎛ ⎞ ⎟⎟ − arctan⎜ ω ⎟ − π − τf ⋅ω ∠L f ,1( j⋅ω) = arctan⎜ ⎟ + arctan⎜⎜ ⎝ a⎠ ⎝ a ⎠ 2/ τ 2 ⎝ f ⎠

(49)

Then ∠Lf ,1(j⋅ω) is decreasing if the derivative of (49) with respect to ω is lower or equal than zero:

d∠L f ,1(j⋅ω) dω

a μ⋅b

= 2

ω +

of the open loop transfer function Lf ,1(s ) = Gf (s )⋅Cf ,1(s ) would be a

Lf ,1(j⋅ω)

< 0 at ω = 0. Then

⎞ ⎛ 1 ⎛ 2⎞ 2 ⎜⎜ a + ⎟⎟ − ⎟⎟ = 0 ⋅⎜⎜ τf ⎠ τf ⎠ ⎝ μ⋅b ⎝

Table 1 shows large variations in the parameters of model (2) in function of the operation regime. Though we propose a gain scheduling control system, uncertainties in the model are always present and a robustness analysis of the designed control system is required. Moreover, the effect on the closed loop performance of actually dealing with a process with a delay instead of a process with a Padé approximation must be assessed. The following lemmas will be useful in this analysis.

function

dt 2

Lemma 4. Consider the process Gf (s ) with model (2) and a controller Cf ,1(s ) given by (41). Then the phase of the frequency response of the

6. Robustness analysis

decreasing

2

< 0 , ω ≥ 0 and Lf ,1(j⋅ω) is a strictly decreasing function. □

dt 2

0.0273⋅s 0.85 1 − 2.45⋅s ε⋅s ⋅ Gf ,1(s ) = ⋅ 1 + 1.82⋅s 0.85 1 + 2.45⋅s α

d Lf ,1

to

2

( ) a μ⋅b

It is easy to check that

= μ/( μ + 1) if it were verified that.

+

2 τf

⎛ 2 ⎞2 ω2 + ⎜ τ ⎟ ⎝ f⎠ 2 τf

⎛ 2 ⎞2 ω2 + ⎜⎜ ⎟⎟ ⎝ τf ⎠

−

τf 2

−

a − τf ≤ 0 ω2 + a2 (50)

≤ 0 , ∀ ω . Then a sufficient

condition to fulfill (50) is that a μ⋅b

⎛ 2 ⎞2 ⎛ 1 ⎞2 ⎜⎜ ⎟⎟ + ⎜ ⎟ >1 ⎝ b⋅μ ⎠ ⎝ a⋅τf ⎠

(44)

Proof: The open loop transfer function (40) using controller (41) is

(

τf

)(

1 + 2 ⋅s ⋅ μ⋅s + 2 ⋅ L f ,1(s ) = Gf (s )⋅Cf ,1(s ) = s⋅(s + a) τf ⋅( μ + 1)

a b

)

(45)

The square of the magnitude of its frequency response is

2

L f ,1 =

⎛ ⎜x + ⎝

4

⎞⎛ 2 ⎟⋅⎜ μ ⋅x + ⎝

τf2 ⎠

a2 ⎞

⎟

b2 ⎠

2

( μ + 1) ⋅x⋅( x + a2)

(46)

where ω2 = x has been substituted. This function has an extreme if

d L f ,1 dt

2

2

⎛ ⎞ 4 a2 8⋅a2 = 0 ⇒ ⎜⎜ 2 + 2 2 − a2⎟⎟⋅x2 + 2 2 2 ⋅x μ ⋅b τf ⋅μ ⋅b ⎝ τf ⎠ 4⋅a4 + 2 2 2 =0 τf ⋅μ ⋅b

ω +

2

( ) a μ⋅b

−

τf a − ≤0 , ∀ω 2 ω2 + a2

(47)

Then Lf ,1(j⋅ω) has an extreme if Eq. (47) has at least one real positive root. A sufficient condition for preventing positive roots in (47) is that all its coefficients be positive. Since the coefficients associated to the powers x1 and x0 are always positive, a sufficient

(51)

Since it can be proven that it is always verified that

τf μ⋅b 1 − < , a a 2

−τf ⋅ s

⋅e

2

(52)

inequality (51) is verified for ω = 0. Then a sufficient condition for the verification of (51) is that the function of ω does not have zero crossings. The condition of zero crossing of (51) can be written as

⎛ ⎛ ⎞ ⎞ 1 ⎟ 2⋅a ⎛ 1 ⎞⎟ 2 a3 ω4 + ⎜ a2⋅⎜⎜ 1 + + ⎜1 − ⎟⎟⋅ω + 2⎟ ⎜ μ⋅b ⎠ μ⋅b ( μ⋅b) ⎠ τf ⎝ ⎝ ⎝ ⎠ ⎛ ⎞ 1 ⎛ 2⎞ 2 ⎜⎜ a + ⎟⎟ − ⎟⎟ = 0 ⋅⎜⎜ τf ⎠ τf ⎠ ⎝ μ⋅b ⎝

(53)

which is Eq. (48). □ It can be proven that our closed loop system verifies both lemmas. Then Lf ,1(j⋅ω) has strictly decreasing magnitude and phase, and crosses the unity circumference only once, exhibiting a clear phase margin. Fig. 9 shows the Nyquist plot of Lf ,1( j⋅ω). Assume that Gf (s ) is the real process and Gf 0(s ) is a model of such process. The closed loop transfer function of the closed loop system shown in Fig. 7(b) implemented with controller (41) and

Please cite this article as: Feliu-Batlle V, et al. Wiener-Hopf optimal control of a hydraulic canal prototype with fractional order dynamics. ISA Transactions (2017), http://dx.doi.org/10.1016/j.isatra.2017.06.012i

V. Feliu-Batlle et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

Fig. 9. Nyquist plot of Lf ,1(j⋅ω) .

Fig. 10. Function Ψ ( j⋅ω) in decibels and its lower limit at high frequencies.

the model Gf 0(s ) is

M (s ) =

The slope at high frequencies of Ψ ( j⋅ω) is 20⋅α0 , which is always

F (s )⋅Cf ,1(s )⋅Gf (s )

(

)

1 + Cf ,1(s )⋅Gf 0(s ) + Cf ,1(s )⋅F1(s )⋅ Gf (s ) − Gf 0(s )

(54)

bigger than the previous value. Condition (57) is therefore always verified at high frequencies. 4) A necessary condition to fulfill (57) is that

whose characteristic equation can be expressed as

Ψmin >

⎛ ⎛ G (s ) ⎞⎞ f − 1⎟⎟⎟⎟ = 0 1 + L f0,1(s )⋅⎜⎜ 1 + F1(s )⋅⎜⎜ ⎝ Gf 0(s ) ⎠⎠ ⎝

(

−1

)

⎛ G ( j⋅ω) ⎞ f > F1( j⋅ω)⋅⎜⎜ − 1⎟⎟ , ∀ ω ⎝ Gf 0( j⋅ω) ⎠

(

⋅ F1( j⋅ω)

−1

)

−1

>

(

) −1 Kf 0⋅(1 + Tf ⋅(j⋅ω) )

Kf ⋅ 1 + Tf 0⋅(j⋅ω)α0 α

asymptotic curve that gives a lower limit of Ψ ( j⋅ω) at high freα0

(

Kf 0⋅ 1 + Tf ⋅( j⋅ωm

(58)

Moreover, based on the magnitude and phase decreasing behavior of L0f ,1( j⋅ω) demonstrated in Lemmas 3 and 4, and on exhaustive simulations, we have observed that in many the cases the verification of condition (58) implies the verification of (57).

7. Gain scheduling system 7.1. Gain scheduling laws The parameters of the fractional order model (2) depend on the operation regime. This is defined by the downstream water level. According to Table 1, the four parameters of the model α , Kf , Tf , τf vary in the following ranges: 0.8 ≤ α ≤ 1, 0.055 ≤ Kf ≤ 0.273, 1.3 ≤ Tf ≤ 2.86, 4.8 ≤ τf ≤ 5.2. The range of variation of the time delay is very small, which justifies the assumption (τf = τf 0) made to obtain condition (57). It will be assumed that the time delay is constant in all the range of operation of the canal an equal to 5 s. This value implies that errors in the estimation of the time delay are equal or below ±4% . The gain of the process is the parameter that experiences the largest variations. The relationship between the gate opening (obtained from Fig. 3) and the steady state downstream end water level is shown in Fig. 11. This data was fitted by the polynomial:

( (57)

Function Ψ ( j⋅ω) is plotted in Fig. 10, expressed in dB. An quencies is 1 + ε⋅( j⋅ω)

) −1 )) α

(56)

Assuming that the time delay is constant (τf = τf 0) and that the parameters that change are Kf , Tf and α , this condition can be expressed as

Ψ ( j⋅ω) = 1 + L f0,1( j⋅ω)

α0

(

Kf ⋅ 1 + Tf 0⋅( j⋅ωm)

(55)

where L0f ,1(s ) = Cf ,1(s )⋅Gf 0(s ). In order to assess the stability robustness of the proposed closed loop system, the following well known result on unstructured multiplicative uncertainties is recalled (e.g. [49]):

1 + L 0f ,1( j⋅ω)

11

/μ, which is also represented in this figure.

Function Ψ ( j⋅ω) has a minimum value Ψmin = 0.397 at ωm = 0.46 rad/s. Some statements can be done about this function: 1) It is independent on the actual values of the parameters of the process. It only depends on the model and the design parameters a , b , ε . 2) Since Ψ ( 0) = 1, condition (57) implies that it must be verified that 0 < Kf /Kf 0 < 2. 3) The maximum slope that can be obtained at high frequencies in the right hand side of condition (57) (expressed in dB) is given by the maximum difference between the assumed and the real fractional orders:20⋅(α0 − α ). This value is 3 dB/dec (see Table 1).

)

y^dwe (u) = 4.5745⋅10−5⋅u3 − 7.6406⋅10−3⋅u2 + 0.4646 ⋅u + 47.7155

(59)

being u the absolute gate opening, and this fitting is also shown in Fig. 11. Then the gain scheduling law is as follows. Assume an initial operation regime defined by a downstream end water level ydwe, i and a gate opening ui . Assume that the setpoint is changed to

Please cite this article as: Feliu-Batlle V, et al. Wiener-Hopf optimal control of a hydraulic canal prototype with fractional order dynamics. ISA Transactions (2017), http://dx.doi.org/10.1016/j.isatra.2017.06.012i

V. Feliu-Batlle et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

12

Fig. 13. Fractional order in function of the initial downstream end water level.

Fig. 11. Downstream water level in function of the gate opening.

* . Then the gain to be used in model (2) is ydwe

Kf =

y*

dwe −1

− ydwe, i

* − ui y^dwe ydwe

(

)

+ 0.005, (60)

where the added constant 0.005 is a correction term. By using the scheduling law (60), (59), errors in the estimation of the process gain remain below ±12% . The process time constant also experiences variations, though these are lower than in the case of the gain. The relationship between the initial downstream end water level ydwe, i and the time constant is shown in Fig. 12. This data was fitted by the polynomial: 2 Tf ( ydwe, i ) = − 0.08465⋅ydwe + 9.207⋅ydwe, i − 248.22 ,i

(61)

Fig. 12 shows that the fitting of the time constant is worse than the fitting of the gain. However, by using the scheduling law (61), errors in the estimation of the process time constant remain below ±25% . The fractional order of the process experiences relatively small variations. The relationship between the initial downstream end water level ydwe, i and the fractional order is shown in Fig. 13. This data was fitted by the polynomial: 2 α( ydwe, i ) = 0.00332⋅ydwe − 0.34⋅ydwe, i + 9.511 ,i

(62)

By using the scheduling law (62), errors in the estimation of the process time constant remain below ±7% .

Fig. 14. Bound of the condition (58) for robust stability.

7.2. Stability robustness analysis of the gain scheduling Condition (58) is used to assess the robust stability of our gain scheduling control system. Fig. 14 shows the upper bound Ψmin and the frequency ωm at which is produced, in function of the design parameter b of the prefilter. This figure shows that Ψmin is an increasing function in the range 0.15 ≤ b ≤ 0.95 (the closed loop system of the process with the real delay is unstable in the range 0 ≤ b ≤ 0.15). Frequency ωm is also an increasing function that remains approximately equal in all this range. In order to fulfill the robustness condition (58), it is convenient to have a value Ψmin as big as possible. Then it is appropriate to use big values of b in order to increase robustness. However, big values of b yield responses with a pronounced overshoot. Then a compromise between these two considerations must be attained. This has led us to use an intermediate value, b = 0.5, through this paper. The level of robustness achieved by our controllers is assessed by defining: (1) the region in which the process parameters may change and (2) a scalar index that scales such region to another region in which the closed loop system is always stable. The region in which the parameters Kf , Tf , α may change was defined in the previous subsection:

(

Kf − Kf 0 Kf 0 where

Fig. 12. Time constant in function of the initial downstream end water level.

(K

f 0,

Tf − Tf 0

≤ δK ,

Tf 0

≤ δT ,

)

α − α0 ≤ δα α0

(63)

)

Tf 0, α0 are given by the gain scheduling laws (60),

(61) and (62) respectively, and δK = 0.12, δT = 0.25 and δα = 0.07. Subsequently, a set of regions of parametric variation are defined in function of a scalar ρ as

Please cite this article as: Feliu-Batlle V, et al. Wiener-Hopf optimal control of a hydraulic canal prototype with fractional order dynamics. ISA Transactions (2017), http://dx.doi.org/10.1016/j.isatra.2017.06.012i

V. Feliu-Batlle et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

13

Fig. 16. Experimental and simulated closed loop downstream end water level responses using the fractional order controller with the gain scheduling system. Fig. 15. Stability robustness index in function of the initial operation regime ( ydwe, i ).

Ω(ρ) ≡

Kf − Kf 0 Kf 0

≤ ρ⋅δK ,

Tf − Tf 0 Tf 0

≤ ρ⋅δ T ,

α − α0 ≤ ρ⋅δα α0

(64)

Then the region of uncertainty obtained in the previous subsection corresponds to ρ = 1. In the following, for each operation regime, the region Ω(ρ) at which the closed loop system remains stable is determined. The necessary condition (58) was used to obtain the largest region of stability. It is given by the maximum ρ that can be attained fulfilling this condition. Subsequently, condition (57) was used to check the validity of the results provided by (58) and slightly correct them if needed. Note that, since condition (57) involves much more calculations that (58) – inequality (57) has to be checked through all the frequencies while inequality (58) only once -, condition (58) was used to determine the maximum ρ while (57) was used to subsequently verify the obtained result. Fig. 15 shows the maximum ρ values attained for each operation regime (defined by ydwe, i ). These values are always higher than 1, which proves that our gain scheduling system with our fractional order controller is robustly stable for all the range of parametric uncertainties obtained in the previous subsection. Moreover, this figure shows that always ρ > 1.07. Then our control system provides a security margin higher than 7% in all the stability regions associated to all the possible operation regimes. Fig. 15 plots the maximum ρ obtained using conditions (57) and (58) for each operation regime. It shows that (58) yields maximum values of ρ that are very close to the ones obtained using condition (57), and justifies the procedure of obtaining first ρ from (58) and later verifying and correcting the obtained values using (57).

8. Experimental results In this section, experimental results have been obtained using the control system developed in this work and the experimental prototype of hydraulic canal described in Section 2. Fig. 16 shows the responses of the closed loop system to several changes of the setpoint. The range of variation of the operation regime in canals (defined by the downstream end water level) is often small. In our prototype, Table 1 shows that the range of variation of the operation regime is from ydwe = 51.5 mm to ydwe = 57.5 mm (only 6 mm). Then 1 or 2 mm are reasonable amplitudes for the setpoint changes. The simulated responses of the process assuming a Padé approximation of order n = 1 of the delay and the same gain scheduling fractional order controller are also overlapped in this

figure. The ±5% band of the final value is represented in this figure and in the following ones. Fig. 16 shows that: 1. The response follows the command signal ydwe (t ) in the setpoint changes. 2. The experimental response is close to the simulated response of the fractional order model with a Padé approximation n = 1 of the delay. The differences are apparent, basically, only in the first τf s in which the experimental response remains in its initial steady state (the delay effect) while the simulated response exhibits variations as consequence of the Padé approximation. Then this approximation, which has been utilized in all the control design procedure, can be regarded as enough accurate. 3. The gain scheduling system reduces the variations in the closed loop response produced by the change of the operation regime. It can be observed that the responses to the first and second set point changes are similar. However, the response to the third set point change is slower though it exhibits about the same damping than the two previous ones. This discrepancy is caused by the difference between the real dynamics and the model given by the gain scheduling system. 4. In the first two set point changes, the experimental responses last τf s more than the simulated ones to reach the command signal, while in the third set point change, the experimental response lasts 4⋅τf s more than the simulated one. 5. A high frequency oscillation appears in the experimental responses. It represents some high frequency dynamics caused by the propagation of a wave on the surface of the canal. This dynamics is of little importance and is not taken into account in the control of the flow of water. Moreover, the amplitude of this wave is about 0.1 − 0.2 mm, which is close to the resolution of the water level sensor. Fig. 17 shows the experimental responses of our canal prototype performing a manoeuvre that starts in an operation regime close to the nominal one (it starts in ydwe, i = 54 mm while the nominal regime is ydwe, i = 54.5 mm), using the integer order and fractional order control schemes. It illustrates that designing an integer order controller from an integer order model of our process leads to a slower response than the one obtained designing the controller from the fractional order model (which is in accordance with what was shown in the simulations of Fig. 8). Fig. 18 shows the experimental closed loop downstream end water level responses in an operating regime different from the nominal one ( ydwe, i = 53 mm) in the cases of using the fractional order controller with and without the gain scheduling system. In the case of not using the gain scheduling, the controller was tuned

Please cite this article as: Feliu-Batlle V, et al. Wiener-Hopf optimal control of a hydraulic canal prototype with fractional order dynamics. ISA Transactions (2017), http://dx.doi.org/10.1016/j.isatra.2017.06.012i

14

V. Feliu-Batlle et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

Fig. 17. Experimental closed loop downstream end water level responses using the integer order controller and the fractional order controller, starting in approximately the nominal operation regime.

Fig. 18. Experimental closed loop downstream end water level responses in an initial operating regime different from the nominal one in the cases of using the fractional order controller with and without the gain scheduling system.

for the nominal operating regime. This figure shows that the response obtained using the gain scheduling system reaches the command signal faster than the response obtained using the fixed controller. This fact illustrates the robustness of the gain scheduling system to changes in the operation regime because it modifies all the parameters of the controller in function of the initial operation regime (in the case of the gain, the final operation regime is also taken into account).

9. Conclusions This article has addressed the control of a laboratory hydraulic canal prototype that has fractional order dynamics and a time delay. Controlling this prototype is relevant since its dynamics closely resembles the dynamics of real main irrigation canals: the principal difference is the time scaling, i.e., time constants and time delays are much smaller in this prototype than in a real canal. Moreover, the dynamics of hydraulic canals (both reals and prototypes) vary largely when the operation regime changes since they are strongly nonlinear systems, which poses a significant difficulty in the design of their controllers. The controller proposed in this article pursued a good time response to set point changes. The design criterium for this controller has been to minimize the integral performance index ISE between a given command signal (which is similar to a smoothed step in order to avoid overshoots) and the process response. The controller that minimizes this index is given by the well known Wiener-Hopf design method of optimal controllers. This method has the following advantages: (a) it minimizes a measure of the

global performance of the time response of the process, (b) it gives a closed form solution which is simple and does not require any iterative procedure, (c) it yields the controller with the optimal structure, i.e, it does not require the “a priori” knowledge of the controller structure, unlike many other design methods and (d) this performance index weights more large errors of the transient response and less small errors when reaching the steady state, which is well suited for our application in which noises and unmodelled high frequency dynamics (as the experimental responses show) make useless to take much care of these small errors. The Wiener-Hopf design method was proposed for standard rational transfer functions and was later extended to fractional order plants without time delay. The main novelties of this article have been therefore: (1) extending the Wiener-Hopf design method to the case of fractional order processes with time delay by using the Padé approximation of the time delay, (2) proposing a new control scheme that overcomes the problem of implementing the unproper controllers yielded by the Wiener-Hopf method, (3) designing its application to the control of an hydraulic canal where, for the first time, a fractional order model of the canal, which is more accurate than the standard integer order one, is used in the control design, (4) since the Wiener Hopf method yields a closed form solution for the controller, a gain scheduling fractional order controller based on the previous one has been proposed that improves the robustness of the control system (minimum ISE is attained – in theory - in all the operation regimes) and 5) demonstrating its feasibility and performance by carrying out experimentation in a laboratory prototype of an hydraulic canal. Then, for the first time, the control of a main irrigation canal taking into account both a fractional order model of its dynamics and its highly nonlinear behavior has been addressed. Besides this, some new theoretical results have been proposed: (1) an analysis of the stability of the proposed controller in the case of using the process with the real delay has been carried out, which yielded conditions for choosing the appropriate Padé approximation (Theorem 2) and (2) thorough analyses of the stability robustness of the Wiener-Hopf controller, alone and combined with the gain scheduling system have been provided. We mention that it is not difficult to implement our controllers since there are many tools and software programs for simulation and real time implementation of fractional order systems and controllers, e.g. [50]. Our next step will be to apply the developed methodology (both characterizing a fractional order model and designing the corresponding controller) to the automation of a real irrigation main canal pool.

Acknowledgements The authors would like to acknowledge the support provided in part by the Consejería de Educación, Cultura y Deportes de la Junta de Comunidades de Castilla-La Mancha (Spain) with the Project POII-2014-014-P, in part by the European Social Fund and in part by the Spanish scholarship FPU14/02256 of the FPU Program of the Ministerio de Educación, Cultura y Deporte.

References [1] Zhang D, Guo P. Integrated agriculture water management optimization model for water saving potential analysis. Agric Water Manag 2016;170:5–19. [2] Figueiredo J, Botto MA, Rijo M. SCADA system with predictive controller applied to irrigation canals. Control Eng Pract 2013;21:870–86. [3] Burt CM. Modern water control and management practices in irrigations: methodologies and criteria for evaluating the impact on performance. In:

Please cite this article as: Feliu-Batlle V, et al. Wiener-Hopf optimal control of a hydraulic canal prototype with fractional order dynamics. ISA Transactions (2017), http://dx.doi.org/10.1016/j.isatra.2017.06.012i

V. Feliu-Batlle et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

[4]

[5] [6] [7]

[8] [9]

[10]

[11] [12] [13]

[14] [15] [16]

[17] [18] [19]

[20] [21]

[22]

[23] [24] [25] [26]

proceedings of the Expert Consultation on Modernization of Irrigation Schemes: Past Experiences and Future Options, Bangkok, Thailand, p. 1–23; 1996. García-Garizábal I, Causapé J, Abrahao R. Application of the irrigation land environmental evaluation tool for flood irrigation management and evaluation of water use. Catena 2011;87(2):260–7. Dai ZY, Li YP. A multistage irrigation water allocation model for agricultural land-use planning under uncertainty. Agric Water Manag 2013;129:69–79. Chaudhry MH. Open-channel flow. 2nd ed. New York: Springer; 2008. Pedregal DJ, Rivas-Perez R, Feliu V, Sanchez L, Linares A. A non-linear forecasting system for the Ebro River at Zaragoza, Spain. Environ Model Softw 2009;24(4):502–9. Litrico X, Fromion V. Modeling and control of hydrosystems. London: Springer Verlag; 2009. Feliu-Batlle V, Rivas-Perez R, Castillo-Garcia FJ, Sanchez-Rodriguez L, Linares Saez A. Robust fractional order controller for irrigation main canal pools with time-varying dynamical parameters. Comput Electron Agric 2011;76(2):205– 17. Rivas-Perez R, Feliu-Batlle V, Castillo-Garcia FJ, Linares-Saez A. Mathematical model for robust control of an irrigation main canal pool. Environ Model Softw 2014;51(1):207–20. Malaterre PO, Rogers DC, Schuurmans J. Classification of canal control algorithms. J Irrig Drain Eng 1998;124(1):3–10. Rogers DC, Goussard J. Canal control algorithms currently in use. J Irrig Drain Eng 1998;124(1):11–5. Castillo‐Garcia FJ, Feliu‐Batlle V, Rivas‐Perez R. Time domain tuning of fractional order controllers combined with a Smith Predictor for automation of water distribution in irrigation main channel pools. Asian J Control 2013;15 (3):819–33. Clemmens AJ, Schuurmans J. Simple optimal downstream feedback canal controllers: theory. J Irrig Drain Eng 2004;130(1):26–34. Wahlin BT. Performance of model predictive control on ASCE test canal 1. J Irrig Drain Eng 2004;130(3):227–38. Schuurmans J, Clemmens AJ, Dijkstra S, Hof A, Brouwer R. Modeling of irrigation and drainage canal for controller design. J Irrig Drain Eng 1999;125 (6):338–44. Wahlin BT, Clemmens AJ. Automatic downstream water-level feedback control of branching canal networks: theory. J Irrig Drain Eng 2006;132(3):198–207. Weyer E. System identification of an open water channel. Control Eng Pract 2001;9(12):1289–99. Rivas-Perez R, Feliu-Batlle V, Sanchez Rodriguez L, Pedregal Tercero DJ, Linares Saez A, Aguilar Mariñosa JV, Langarita Garcia P. Identification of the first pool of the Imperial de Aragon main irrigation canal. Ing Hidraul En Mex 2008;23 (1):71–87. Aleem SA, Muhammad A, Nasir HA. System identification of distributory canals in the Indus basin. IFAC Proceedings Volumes; 47(3): 8743–8748; 2014. Malaterre PO, Baume JP. Modeling and regulation of irrigation canal: existing applications and on-going researches. In: Proceedings of IEEE International Conference on Systems, Man, and Cybernetics, San Diego, USA; 3850–3855; 1998. Feliu-Batlle V, Rivas-Perez R, Castillo-Garcia FG, Sanchez-Rodríguez L. Smith predictor based robust fractional order control: application to water distribution in a main irrigation canal pool. J Process Control 2009;19(3):506–19. Podlubny I. Fractional differential equations.San Diego: Academic Press; 1998. Monje CA, Chen YQ, Vinagre BM, Xue DY, Feliu V. Fractional-order systems and controls. Fundamentals and applications.London: Springer; 2010. Martinez-Gonzales R, Bolea Y, Grau A, Martinez-Garcia H. An LPV fractional model for canal control. Math Probl Eng 2009 [Article ID471540]. Calderon-Valdez SN, Feliu-Batlle V, Rivas-Perez R. Fractional-order mathematical model of an irrigation main canal pool. Span J Agric Res 2015;13(3):

15

e0212. [27] Podlubny I. Fractional order systems and PI-λ-D-m controllers. IEEE Trans Autom Control 1999;44(1):208–14. [28] Tavakoli-Kakhki M, Haeri M, Tavazoei MS. Simple fractional order model structures and their applications in control system design. Eur J Control 2010;16(6):680–94. [29] Luo Y, Chen YQ, Wang CY, Pi YG. Tuning fractional order proportional integral controllers for fractional order systems. J Process Control 2010;20(7):823–31. [30] Basiri MH, Tavazoei MS. On robust control of fractional order plants: invariant phase margin. J Comput Nonlinear Dyn 2015;10(5) [article number: 054504]. [31] Bedri V, Tavazoei MS. Some analytical results on tuning fractional-order [proportional-integral] controllers for fractional-order systems. IEEE Trans Control Syst Technol 2016;24(3):1059–66. [32] Buslowicz M. Stability of linear continuous-time fractional order systems with delays of the retarded type. Bull Pol Acad Sci 2008 2008;56(4):319–24. [33] Hwang C, Cheng YC. A numerical algorithm for stability testing of fractional delay systems. Automatica 2006;42(5):825–31. [34] Hamamci SE. An algorithm for stabilization of fractional-order time delay systems using fractional-order PID controllers. IEEE Trans Autom Control 2007;52(10):1964–9. [35] Özbaya H, Bonnet C, Fioravanti AR. PID controller design for fractional-order systems with time delays. Syst Control Lett 2012;61(1):18–23. [36] Oustaloup A. La commande CRONE: commande robuste d'ordre non entir. CHermes; 1991. [37] Maione G, Lino P. New tuning rules for fractional PIα controllers. Nonlinear Dyn 2007;49(1):251–7. [38] Monje CA, Vinagre BM, Feliu V, Chen YQ. Tuning and auto-tuning of fractional order controllers for industry applications. Control Eng Pract 2008;16(7):798– 812. [39] Valerio D, da Costa JS. Tuning of fractional PID controllers with Ziegler–Nichols-type rules. Signal Process 2006;86(10):2771–84. [40] Padula F, Vilanova R, Visioli A. optimization-based fractional-order PID controllers design. Int J Robust Nonlinear Control 2014;24(17):3009–26. [41] Rhouma A, Bouani F, Bouzouita B, Ksouri M. Model predictive control of fractional order systems. J Comput Nonlinear Dyn 2014;9(3):031011. http://dx. doi.org/10.1115/1.4026493. [42] Feliu-Batlle V, Rivas-Perez R, Sanchez-Rodriguez L. Fractional robust control of main irrigation canals with variable dynamic parameters. Control Eng Pract 2007;15(6):673–86. [43] Rivas-Perez R, Feliu-Batlle V, Castillo-Garcia FJ, Sanchez-Rodriguez L, Linares A. Robust fractional order controller implemented in the first pool of the Imperial de Aragon main canal. Tecnol Y Cienc Del Agua 2014;5(1):23–42. [44] Martínez-González R, Bolea Y, Grau A, Martinez-Garcia H. LPV non-integer order model and control of an irrigation canal using gain scheduling technique. In: Proceedings of the 35th Annual Conference of IEEE Industrial Electronics, Porto, Portugal; 1615–1620; 2009. [45] O’Dwyer A. Handbook of PI and PID controller tuning rules. Imperial College Press. London; 2003. [46] Aliev FA, Larin VB. Optimization of linear control systems. New York: Gordon & Breach; 1998. [47] Vinagre BM, Feliu V. Optimal fractional controllers for rational order systems: a special case of the Wiener-Hopf spectral factorization method. IEEE Trans Autom Control 2007;52(12):2385–9. [48] Lam J. Model reduction of delay systems using pade approximants. Int J Control 1993;57(2):377–91. [49] Sinha A. Linear systems. Optimal and robust control. CRC Press. Boca Raton, FL; 2007. [50] Vinagre BM, Feliu-Batlle V, Tejado I. Fractional Control: Fundamentals and User Guide. RIAI-Revista Iberoamericana de Automática e Informática Industrial; 2016. 13(3). p. 265–280.

Please cite this article as: Feliu-Batlle V, et al. Wiener-Hopf optimal control of a hydraulic canal prototype with fractional order dynamics. ISA Transactions (2017), http://dx.doi.org/10.1016/j.isatra.2017.06.012i