Abstract
Canal systems are complex nonlinear, distributed parameter systems with changing parameters according to the operating point. In this paper, a linear parameter-varying (LPV) state-space canal control model is obtained by identification in a local way using a multimodel approach. This LPV identification procedure is based on subspace methods for different operating points of an irrigation canal covering the full operation range. Different subspace algorithms have been used and compared. The model that best represents the canal behavior in a precise manner has been chosen, and it has been validated by error functions and analysis correlation of residuals in a laboratory multireach pilot canal providing satisfactory results.
1. Introduction
Water is one of the most used resources by industrial and agricultural sectors, and obviously by population. One fundamental use of water is the irrigation activity, and one the main challenges in this area is to prevent water losses and to permit an efficient use of this scarce and vital resource. These aspects have led to the usage of automatic control systems and the implementation of different advanced control algorithms for the regulation of open-flow irrigation canals. Hence, those control techniques will allow fulfilling the desired performance and the ecological flow in irrigation as well as saving water at the same time. To design an effective controller, a good control model is needed. Therefore, advanced process modeling techniques are required to make an accurate control model.
Modeling and control of nonlinear complex systems is a challenging task. Nonlinear effects can no longer be neglected to meet the specifications imposed on today’s complex control systems. Unfortunately, the higher the complexity, the lower our ability to deal with it and to understand it. Open-flow canals are complex systems, that is, they are large distributed parameter systems that have the following main characteristics: nonlinear behavior and dependence of the parameters with the operating point and coupling among pools [1]. This type of systems can be fully described by Saint-Venant’s equations [2]. This representation is the most used model to describe the physical dynamics of a real irrigation canal. However, this complex model is based on a nonlinear hyperbolic partial differential equations system that has analytical solution only in very special cases, requiring the use of numerical methods to solve it properly [3]. This complex representation of the system is suitable for simulation models, but it is not suitable to design controllers that fulfill the control design needs. Then, linearization and simplification of Saint-Venant’s equations are currently studied by the irrigation researchers’ community of control [4] to develop simpler control models.
Distributed parameter systems with a very large number of states, that is, systems with coupling, have been approximated by decoupled low-order linear time invariant (LTI) models in order to use classical linear control design tools, as a usual practice in control engineering. In fact, control researchers’ community has usually used linear control techniques (such as fuzzy control [5], robust control [6], etc.), and even nonlinear control approaches (such as, sliding control [7–9], etc.) for this kind of systems. LTI control models widely used are Hayami model [10], Muskingum model [3], IDZ model [4, 11], or black-box models identified using parameter estimation by classical identification methods [1, 12, 13]. However, these systems are not completely amenable using conventional linear modeling approaches due to the lack of precise, formal knowledge about the system; strongly nonlinear behavior; high degree of uncertainty; time varying characteristics; dynamic parameters changing over the operating point and coupling between pools. Then, simplified control model structures are needed preserving their information. Taking into account these previous properties, a linear parameter varying (LPV) model is required, which consists in a model that regards both the parameter and delay variations with respect to the operating points. In this way, the system information is preserved while it would be lost with a linear control model. These LPV control models permit the design and computation of LPV controllers that rigorously guarantee the system stability and performance [14] for smooth variations of system parameters as well as abrupt ones [15]; this is the case of irrigation canals. The preferred representation scheme for complex plants (multivariable systems involving large system orders) is a state-space model. Then, subspace-based system identification methods are a branch that has been recently developed in system identification attracting much attention thanks to their computational simplicity and effectiveness in identifying dynamic state-space linear multivariable systems. These algorithms are numerically robust and do not involve nonlinear optimization techniques, being fast and accurate.
Due to applications of large dimensions commonly found in industrial processes, subspace identification methods are very promising in this field. In the basis of the aspects explained above, in our case (a multireach canal system) a state-space model representation is suitable instead of a transfer function system description. For this reason, an LPV state-space control model has been developed through LPV identification techniques.
Besides, since system canals are nonlinear, a common engineering approach to deal with this complexity is the divide-and-conquer strategy: decompose the complex problem into several subproblems easier to solve. According to this previous strategy, a method to model complex nonlinear systems has arisen. It is based on partitioning the whole operating range of the nonlinear system into multiple, smaller operating regimes and modeling the system for each and every of these regimes. The task of finding a complete global model for the system is thus replaced by determining linear local models and subsequently combining these local models into a global description for the system obtaining the LPV model (local LPV identification). This multiple model approach is often referred to as operating regime decomposition [16] or multimodel identification. Interactions between the relevant system phenomena are less complex locally than globally. From the divide-and-conquer point of view, it is desirable to choose the local models such that they are less complex than the global nonlinear model. It is expected that the simpler the local models are, the more models and thus the more operating regimes are needed to describe the global system accurately enough. A trade-off has to be made, because too simple local models lead to an explosion of the number of operating regimes needed. Another manner to identify the nonlinear systems is by the use of specific identification methods in a global way in one-shot through optimization [17] (global LPV identification). Although it is possible to use nonlinear local models (see, e.g., [18]) a common choice is to use local linear models. The main reasons for this choice are a solid theory for linear systems has been developed over the years, linear models are easy to understand, and they are widely used by engineers.
Specifically, the main contribution of the paper is the design of a low-order LPV state-space multivariable control model describing the water flow dynamics in a multireach irrigation canal. The model is estimated over the full operation point range using local LPV identification. Several subspace identification methods are applied and their performances are compared in order to select the subspace algorithm that yields the best control model. This model will be suitable to design LPV controllers that will warranty the stability and the desired performance around all operating points of the system with rigorous formality [14, 15].
The structure of this paper is as follows: in the next section, the main issues related to LPV identification and subspace identification methods used in this study are presented. Section 3 briefly presents the two reach irrigation canal used in this research. Section 4 discusses the important steps (generation and pretreatment of the data set, order estimation, performance quality criteria, and global model obtained by interpolation) in developing a suitable LPV subspace model for the system and compares the performance of the used subspace identification algorithms (N4SID, MOESP, and CVA) carrying out the model validation. Finally, Section 5 provides conclusions.
2. LPV Subspace Identification Methods
Linear time-invariant (LTI) models are not suitable to control systems such as open-flow canals with coupled pools and distributed parameters with nonlinear behaviour that depend on the operating points. However, by an LPV model (with varying parameters depending on the operating points), to preserve the aforementioned information of the system (nonlinearity, coupling, etc.) is also possible, obtaining a more accurate and faithful behaviour with the reality. As it has been emphasized in Section 1, there are two procedures to carry out LPV identification: (i) multimodel identification (local LPV identification) and (ii) one-shot LPV identification (global LPV identification). The former approach consists in a two-step procedure where (1) LTI models at several different equilibrium (operating condition) are identified by classical methods [13, 19]; (2) a global multi-model is obtained by interpolation among the local LTI models, and different interpolation techniques can be used such as membership fuzzy functions [20], polynomial interpolation [21], among others. The latter approach consists in carrying out a one-shot identification in a global way as proposed in [17]. The local approach has the important practical advantage that many engineers are well experienced in LTI identification experiments and that the local LTI models can be estimated using a wide variety of well-established and widely spread LTI identification algorithms. To properly interpolate these local models, all local LPV identification techniques require that the local models are represented in a consistent state-space form.
The discrete-time subspace identification methods refer to a kind of algorithms which allow identifying a robust and reliable state-space model of MIMO linear systems estimating state sequences directly from the input-output measurements. Based on orthogonal or oblique projected subspaces generated by the rows or columns of Hankel matrices of the input-output data, the process is followed by a singular value decomposition so as to determine the order of the model and the observability matrices. Finally, the state-space model is obtained through the solution of a least squares problem. Subspace-based methods for state-space modeling have their origin in state-space realization, as developed by [22]. The term “subspace identification method” was introduced in the early 90s. The subspace identification can use many different versions of subspace methods such as Canonical Variate Analysis (CVA), Multivariable Output-Error State-Space model identification (MOESP), State-Space System Identification (N4SID), Canonical Correlation Analysis (CCA), and Deterministic and Stochastic Subspace System Identification and Realization (DSR) [10]. These algorithms attract much attention because they present many advantages: their computational simplicity and effectiveness to determinate dynamic linear multivariable systems. Nevertheless, a drawback which can be noticed is that these algorithms require a large amount of data to build accurate models. So the experiments to collect data can be large and time consuming. For this reason in control problems usually an off-line identification is used.
LPV models obtained using subspace identification methods are mathematically described by the following form: where the vectors and are the observations at the discrete time of inputs and outputs of the process, respectively. The vector represents the state vector of the process at discrete time instant and contains the numerical values of states, and is the parameter vector.
LPV system can be viewed as a nonlinear system that is linearized along a time-varying trajectory determined by the time-varying parameter vector . Hence, the time-varying parameter vector of an LPV system corresponds to the operating point of the nonlinear system. In the LPV framework, it is assumed that this parameter is measurable for control. In many industrial applications, such as process control, the operating point can indeed be determined from measurements, making the LPV approach viable. Control design for LPV systems is an active research area. Within the LPV framework, systematic techniques for designing gain-scheduled controllers can be developed. Such techniques allow tighter performance bounds and can deal with fast variations of the operating point. Furthermore, control design for LPV systems has a close connection with modern robust control techniques based on linear fractional transformations [23]. The important role of LPV systems in control system design motivates the development of identification techniques that can determine such systems from measured data, the LPV identification [23, 24].
Local linear modeling is one of the many possibilities to approximate a nonlinear dynamic system. It is based on partitioning the full operating range of the nonlinear system into multiple, smaller operating regimes and modeling the system for each of these regimes by a linear model. By making a weighted combination of these linear models, it is expected to describe the complete nonlinear behaviour accurate enough. The local lineal model structure is where is the number of local models and and the other variables have the same meaning that in (2.1). The weighting vectors are unknown functions of the scheduling vector . This scheduling vector corresponds to the operating point of the system. This system is closely related to the LPV system in (2.1). The weighting functions can be interpreted as model validity functions: they indicate which model or combination of models is active for a certain operating regime of the system. A weighted combination of local linear models can be used to approximate a smooth nonlinear system up to an arbitrary accuracy by increasing the number of local models. As it is stated, the local linear state-space system is closely related to the LPV system: consider that time-varying parameters are known, while for the local linear model structure the weighting functions have to be determined from input and output data [24].
In this work, the plant is identified by several local LPV subspace identification methods (cited above) to estimate the space-state representation that describes the system dynamics suitably. The LPV identification method used for the experimental modeling of our pilot two-pool canal (presented in the next section) is a two-step procedure where (1) linear state-space models are identified at several different operating points by subspace identification methods over the full range of operation; (2) a global state-space multi-model is obtained at the end interpolating the local state-space models using polynomials [21]. In this paper, the following identification methods have been used: N4SID, the standard method (i.e., N4SID1 from now on) and the robust method (i.e., N4SID2 from now on), CVA algorithm, and MOESP procedure [10]. These methods are used to estimate the model in each operating point. The local identification method forces the local models to fit the system separately and locally. The steps of the identification procedure are explained in Section 4.
3. Description of the Process
An experimental canal prototype is used in this research. This canal consists in two tanks TT and TD with cross section A. The full structure of the canal prototype is presented in Figure 1 (up) and Figure 2. The two tanks are serially connected with pipes, and they are slightly tilted to allow the flow of the water. There is a reservoir at the bottom of the plant to supply water to the tanks. A first pump (with a flow of 3,800 liters/h), named in Figure 1 (up), permits to collect water from the reservoir to fill the upper tank TT. An ultrasound sensor attached to the metallic structure at the end of the first tank, named , measures the water level. A second pump, named , (with a flow of 1,300 liters/h) allows to drain the upper tank and to fill the lower tank TD.
A second ultrasound sensor positioned also at the end of the tank, named , enables to measure the water level. Finally, a third pump, called , gives the possibility to drain away the last tank to the reservoir. The sensors and the pumps are directly connected to an electronic board, which allows to power them and to exchange data between a Matlab program (executed on a PC shown in Figure 2) and the canal. The program uses the Real Time Windows Target (Matlab toolbox) to communicate with the electronic board and the canal. Finally, it is important to emphasize that the prototype canal constitutes a closed system. The water leaves from the bottom reservoir to the tank TT and arrives to the tank TD via the pump . Then, the water returns to the reservoir via the pump . That constitutes a coupled system where the first tank has a big influence on the second one.
There are many methacrylate plates along the two tanks TT and TD shown in Figure 1 (down). These plates are 2 centimeters apart creating a zigzag path in each tank. That provides a delay in the water to reach the other extremity of the tank where the sensor is located. This delay has to be taken into account for the identification. The delay changes depending on the water level in the canal. The more the water in the canal, the smaller the delay. This fact justifies the use of LPV identification to deal with this problem. Each sensor measures the canal level at the end of its path with a precision of 1 mm. The maximum allowed level is 15 cm to avoid the overflow.
4. Identification of an LPV Subspace Model for Two-Reach Pilot Canal
In this section, the different subspace estimation methods in an LPV local way (see Section 2) are applied to the multivariable pilot canal plant. The identification process is carried out following the steps: (1) design of the experiment, collection of input-output data in each operating point (taking into account a suitable selection of the excitation input of the system), and pretreatment of these data; due to the system delays, those delays are estimated and removed from the identification data; (2) selection of the model order via different criteria, singular value decomposition (SVD), and Akaike information criterion (AIC); (3) estimation of the local space-state models for each operating point by the aforementioned subspace identification methods and their interpolation to obtain a global model by Nearest Neighbour interpolation [10, 25]; (4) validation of the model in all the operation range by error functions (MRSE and MVAF) and correlation analysis of residuals. The identification results of each estimation algorithm (N4SID1, N4SID2, CVA, and MOESP) are compared and studied. When canal models are identified, two problems have to be considered: the problem of large and variable delays and the nonlinearity of the dynamics and its variability with the operating points [26, 27]. These problems are separately treated. The parameters models are estimated using the state-space algorithms in each operating point without the delay effect, previously calibrated by correlation method [13].
4.1. Generation and Pretreatment of the Data Set
It is not an easy task to select either the input or the output variables of a process. In this experiment, two inputs () and two outputs () are considered. The input corresponds to the voltage of the first canal pump; corresponds to the voltage of the second pump; the output is the downstream level of pool 1, and the downstream level of the pool 2. Pseudorandom Binary Sequences (PBRS) are widely used in the identification process [13]. These signals are persistent input signals that contain a large number of frequencies representative of the dynamics of the plant. In order to choose the number of operation points (OPs) of the canal plant in a rigorous manner, the optimized OP multipoint technique is used [28]. In this paper, four equidistant operating points have been used. The local model identification will be performed in every operating point because the system is not linear.
In our experiment, a pulse generator creates a train of PBRS as pump input voltage signals which adequately excite the system at different operating points. For the first pump, the PBRS signal, , changes the pump opening at intervals of 800 s and for a period of 10 s. The identification procedure was carried out off-line using 3200 samples of the data set. In Figure 3, the train of PBRS signals of each pump is shown.
(a)
(b)
(c)
4.2. Calibration of the Delay and Order Model Selection
In this subsection, the delay in both pools will be estimated. It is known that open-flow canals present large delays that change with the operating point (in this case, the pump operation, i.e., the upstream level of each canal pool) [26]. In this work, the delay estimation has been derived using correlation method [1, 13]. The delay estimation error is equal or less than 1% in all cases. Next, the identification will be performed having removed the delays.
In subspace-based algorithms, the determination of the model order () can be complex. This order can be determined calculating the number of singular values (SVD) different from zero of the orthogonal (or oblique) projections of row spaces of data block-Hankel matrices. However, it is difficult to calculate it when the system data are corrupted by noise. It is also not straightforward to calculate this number, so that the decision is taken by detecting a gap in the spectrum of the singular values. As it can be seen in Figure 4 for N4SID2, the gap is difficult to determinate and hence the application of this strategy becomes really subjective. Therefore, the decision of the model order will be taken with the following criterion.
A reliable technique is Akaike’s final prediction error (FPE) criterion and his related Akaike Information Criterion (AIC). AIC procedure allows determining the order of a system, and it is defined as where is the loss function (quadratic fit) for the structure, is the length of the data, and is the total number of estimated parameters. Using the AIC criterion, the best order model is given by the minimum AIC() value. In Figure 4, the results of this criterion using N4SID2 algorithm on each operating point can be seen. A second-order model has been selected by the AIC’s criterion for all the above subspace identification methods in the full operation range of the canal. Note that in Figure 5 a substantial difference between first and second order models can be seen, but choosing higher order models is irrelevant for the results using N4SID2. Therefore, a second-order model is enough to achieve a suitable control model.
The above fact demonstrates why engineers widely accept a first-order model with delay (IDZ model) or a second-order model with delay (Hayami model) for an irrigation canal approximation. As it is stated in [11, 29], the chosen model structure depends on the celerity coefficients, diffusion, and length of the canal. In our case, the canal is large enough to consider a second-order model, corroborating the model order selection by the chosen statistical method.
4.3. Results and Performance Quality Criteria
After having chosen the best order , the goal is to determine the best algorithm to obtain the final model. In this study, various subspace algorithms have been tested in the full operation range along the identification process. In order to select the best algorithm for identification, these methodologies have been compared among them to see which one fits best with data. The degree of adaptation of each one has been quantified by means of a cross-validation method, using the following typical performance indicators: MRSE (mean relative square error) and MVAF (mean variance-accounted for).
MRSE:
MVAF: With being the th real output, the th simulated output produced by the model, and is the number of repetitions of the experiment. The MRSE index given by (4.2) is used to measure the mean relative square error between the real process outputs and the outputs predicted by the model. As stated by (4.2), an MRSE index of 0 indicates a perfect model. MVAF in (4.3) is a measure for evaluating the dynamic properties of the produced models. If the ratio of variance /variance () is small, the MVAF is close to 1. This index constitutes a quantitative measure of the model quality.
An experiment around the full the operating point range has been carried out and the best method has been selected. In Table 1, the goodness of the each algorithm (N4SID1, N4SID2, MOESP, and CVA) is shown. We can observe that all the methods provide a suitable prediction to obtain a control model.
In average, N4SID2 is the most precise method in the full system operating range. However, in the central operating points, MOESP and CVA have a 4% more of precision, and in the highest operating point N4SID1 is 1% more precise. In Figures 6 and 7, the comparison between the measured and predicted downstream level of both pools is shown for a specific set-up around the full system operation range demonstrating the goodness and precision of the selected subspace identification method.
(a)
(b)
(c)
(a)
(b)
(c)
Those results were already expected because MOESP algorithm has the inherent drawback that it estimates the state sequence using a certain past window, possibly leading to biased results. Similar approximations are made in the subspace LTI algorithm N4SID; however, by making the past window larger, and larger this bias will tend to zero.
Apart from this error function used to validate the identified model, a correlation analysis of residuals is required. One of the most basic tests [2] is to compute the correlation between the regressors, the past inputs in this case, and the residuals:
It is usual to plot these estimates as a function of and compare with their standard deviations to check if they are significantly different from zero. If not, there is no significant influence of input in , so it is not possible to say the estimated model has not picked up all the influence of on (the input on the output). It is supposed the assumption that is a white noise with variance and zero mean. The result is typically presented as a plot of the autocorrelation of the residuals and a plot of the cross-correlation between the inputs and the residuals. In Figure 8, auto- and cross-correlation functions with uncertainty regions for both pools are presented. It can be observed that the auto- and cross-correlation are within the regressors standard deviations providing a model that reproduces pretty well the main characteristics of dynamics of the pilot canal complex process.
(a)
(b)
Finally, the LPV canal global model is obtained by the use of N4SID2 method. The local models obtained with this algorithm are combined by interpolation to create the global model. Therefore, the parameters of the estimated model, that is, the matrices , and are interpolated by Nearest Neighbor Interpolation algorithm obtaining the (2.2). This algorithm is a numerical method widely extended by the scientific community [10, 25]. This method sets the value of an interpolated point to the value of the nearest data point. The results of this interpolation can be seen in Figure 9. It can be observed that each parameter of the system depends on the operating point, the gate openings ().
| (a) |
| (b) |
| (c) |
| (d) |
5. Conclusions
This paper introduces an approach for approximate modeling of distributed parameter processes using LPV identification. The use of LPV models allows the system to be approximated by multiple local low-order models combined by interpolation. Here, specifically an LPV MIMO state-space model for an irrigation canal is identified. The LPV identification is carried out with several subspace-based algorithms (N4SID, robust N4SID, MOESP, and CVA) in a local way. These methods have been compared and the most accurate in the full operation range has been selected. The model has been validated in a laboratory pilot multi-reach canal obtaining very good results.
Acknowledgment
This work has been funded by the Spanish Ministry of Science Project DPI2010-17112.