Abstract
This paper deals with the identification of two-time scale linear dynamic systems, which are an important class of multiscale systems. Classical identification processes may fail to yield accurate parameters for systems of this class and, for this reason, the authors propose two different techniques to estimate the system parameters. The first technique utilizes two prefilters that are iteratively tuned. The second one considers wavelet filters that are tuned based on the results of the first iterative algorithm. Identification and analysis results for a dynamical aircraft model are shown to demonstrate the algorithm’s performance.
1. Introduction
The development of system models is an essential task to any branch of science. Being a complex task, the construction of a model based on observations of the modelled system is most often developed from a single view or observation scale. Taking the example of the human body, models are developed considering scales of observation on the level of organs, tissues, cells, and molecules. The human body system works from the harmonious integration of subsystems observed on these various levels. However, the construction of a model that integrates all these scales of observation is not a trivial task. As exemplified with the human body system, complex dynamic systems can be modelled per scale of observation. The construction of several models that consider their various scales of observation and subsequent integration is what is called multiscale modelling [1].
Complex systems are characterized by a hierarchical multiscale nature with respect to not only space but also time [1]. Considering this nature of the systems, methodologies of study and modelling are in demand. Any multiscale methodology must consider the following issues: correlation between phenomena at different scales, trade-off between different dominant mechanisms, coupling between spatial and temporal structural changes, and critical phenomena occurring in complex systems.
Two-time scaled systems (TTSS) are a particular and important case of multiscale systems because several physical systems such as batteries [2], aircraft longitudinal dynamical model [3], and thermal building models [4] present two-scale behaviour. In this particular case, the system has two well defined time scales (one slow and the other fast) and one of the objectives is to identify these dynamics [5, 6]. Spatial scales are not considered in this paper.
Works related to two-time scaled systems have particular importance because traditional prediction error methods (PEMs) tends to overemphasize dynamical modes that affect the overall model response more heavily [2, 7–9].
A two-time scaled system identification procedure presented in [7] is one of the main related works. The authors of that paper show that the classical least squares (LS) method typically fails to give an accurate estimation because the data of a two-time scaled system is scattered in the frequency domain. In that paper, the authors propose a new technique to identify this kind of system and present a second-order two-time scaled system as example. The approach proposed by the authors is compared with the classic least squares method.
In [8, 9], the authors show that measurement noise can result in further model inaccuracies. The application considered therein was a model for a battery system (applied in electric vehicles) that exhibits two-time scales. In batteries, the main two phenomenological effects [2] are charge transfer and diffusion, two of the most important effects in battery dynamics. Charge transfer occurs on a time scale of 0.1–100 Hz, while diffusion occurs from around 1 Hz down to 0.001 Hz. The authors in [2] also compare their technique with the LS method and show that the proposed method produced better estimated models.
Both [2, 7] consider a prefiltering data mechanism to separate the slow and fast dynamics. In both papers, the prefilter design is carried out based on previous knowledge of the system’s characteristics. The cut-off frequencies of the filters are calculated offline and kept fixed during the estimation procedure.
The contribution of this paper is to present techniques that estimate two-time scaled models. The first one is an iterative algorithm that computes cut-off frequencies for classical prefilters and these frequencies are corrected during the iterations. The second one takes the result (cut-off frequencies) of the mentioned iterative algorithm and uses it as parameter of a second estimation stage. This second estimation stage is based on wavelets. A Lockheed F104G aircraft model example is shown to demonstrate the efficiency of the proposed procedures. This paper is organized as follows: the Materials and Methods section summarizes a theoretical background on two-time scaled system identification, presents the algorithm based on classic filters, and shows the framework that considers wavelet filters using the results of the iterative algorithm; the section Application Example and Results shows results obtained for an example, and the last section brings the conclusions.
2. Materials and Methods
In this section, we will briefly define TTSS. The definition is presented in continuous time because of the simplicity of the formulation and its easy extension to discrete time.
2.1. Definition of Two-Time Scaled System
TTSS is a system that can be characterized by two well defined scales, one fast and the other slow. Mathematical descriptions of this kind on system may be considered in time domain (singular perturbation theory) or in frequency domain. In frequency domain, the transfer function is written as a function of a small parameter called scale parameter. Reference [5] gives a complete characterization of two-frequency scale transfer functions. In that paper, its authors describe a transfer function matrix (MIMO) as a function of that may be characterized as two-time scaled given certain conditions. In that paper, the class of systems that can be described by that approach is very wide. In this paper, the mathematical formulation is based on [5], but is more restricted and will be described in the following.
Consider the following model in transfer function form [7, 10]:where models the slow system dynamics scaled by the parameter and corresponds to the fast dynamics. Considering the frequency response of the system, for high frequencies we can assume thatand for low frequencies
To keep these approximations, we will assume that the static gain of is nonzero and that the high-frequency gain of Ts is neither zero nor infinity. Adopting these assumptions, we shall further consider that Ts is biproper. From a numerical point of view, Ts must contain both slow zeros and poles, with some of these still left in Tf. This is to allow near cancellation in Ts for higher frequencies.
The described system (TTSS) has a frequency response whose Bode plot is monotonic over a “wide” frequency bandwidth and such that the approximations (2) and (3) are valid. In [5] the mathematical description of TTSS is not restricted to the monotonic case.
2.2. Prefiltering in the Identification Process
As mentioned before, prefilters are a classic solution to avoid inaccurate results when identifying TTSS. The use of prefilters may be illustrated by Figure 1.
Prefiltering has been known to increase the performance of identification schemes by transforming the shape of noise and reconditioning the excitation matrix as we can see in [5, 11, 12]. Considering this, we can define the low-pass and high-pass prefilters in the following way:where ns is the degree of Ts. Assumptions (4) and (5) are fulfilled by taking Fh = 0 in low frequency and Fl = 0 for sufficiently high frequencies, low and high frequencies being adequately defined in the context of the respective application.
In this paper, the first contribution considers classic filters. The proposed iterative algorithm supplies as result a good choice of the cut-off frequencies of the prefilter. The second contribution consists of wavelet filters that work in a frequency range defined by the mentioned iterative algorithm. A more detailed discussion about the prefilters and their structure will be given in the following sections.
2.3. Iterative Design of the Prefilters
The existing two-time scale identification methods [2, 7] consider that the classical identification requirements must be valid: the input signal must be persistently exciting; a previous structure information must be known; and an estimation algorithm must converge with minimum error covariance.
A basic idea of [2, 7] is to find prefilter that suits the input signal. The task of the prefilter is to supply an input signal that excites the system in order to highlight the desired dynamics and minimize the effect of nondesired dynamics. Usually, to highlight the dynamics, a low-pass filter is designed, and to highlight the dynamics, a high-pass filter is designed. The cut-off frequency of the low-pass filter should be high enough to encompass the passband of the slower mode, but should be low enough to contain as little of the passband of the faster mode as possible [2]. On the other hand, the passband of the high-pass filter needs to overlap with passband of the faster dynamics but has as little overlap with the passband of the slower dynamics as possible [2].
2.4. Description of the Iterative Algorithm
The cut-off frequency of the filter must maintain some relation with the cut-off frequencies of the system. Considering this affirmation, the proposed algorithm estimates (without filtering of the excitation signal) an initial model (possibly wrong) in order to obtain an initial cut-off frequency. This initial cut-off frequency will be considered for both, the low-pass and the high-pass filters. The system is separately excited with the two new signals (i.e., filtered signals) and the new data is then considered for a second estimation process (separately). The result of this second estimation process results in the first approximation of Ts and Tf that will be called Ts0 and Tf0. The cut-off frequencies of Ts0 and Tf0 will entail other low-pass and high-pass filters that will be considered for a third estimation pass. The algorithm will repeat the described process until a given stop criterion is satisfied. The estimation algorithm considered here is the classic least squares (LS) algorithm.
The recursive procedure of the proposed idea is described below.
Step 0. Set i0.
Step 1. Apply a persistent exciting signal u to the system and sample its output.
Step 2. Estimate an initial model T0 and compute its cut-off frequency (f0).
Step 3. Set and .
Step 4. Design a low-pass filter and a high-pass filter with cut-off frequencies fli and fhi, respectively.
Step 5. Generate two signals usi and ufi by filtering u with Fli and Fhi, respectively.
Step 6. Apply the signals usi and usi to the system (separately) and sample its outputs ysi and yfi.
Step 7. Estimate Tsi and Tfi using LS algorithm.
Step 8. Set ii+1.
Step 9. Compute the cut-off frequencies fli and fhi considering Tsi-1 and Tfi-1.
Step 10. If the stop criterion (see in the following) is not satisfied, go back to Step 4; else stop.
In this iterative procedure, the last estimated model is the model that will be considered for validation. The cut-off frequencies are calculated by the following expression:
The stopping criterion is based on the convergence of the cut-off frequencies. When both frequencies converge, the algorithm stops. To evaluate the convergence of the frequencies, the algorithm calculates the absolute difference between the current frequency and the frequency in the last step as in (7) and (8).
The tolerance value ϵ is used to define the stopping criteria; i.e., when dli < and dhi < ϵ, the iterative procedure stops. When dli < and dhi < , the prefilters will not change their cut-off frequency significantly and the new estimation will not be significantly different from the previous estimation. Furthermore, if the cut-off frequency does not change very significantly, the estimated parameters also will not. This assertion is not general and is valid to the class of systems described in this paper. The proposed algorithm is based on the classic LS. The convergence and consistency of this method are well established. In the filtered case, the convergence and consistency of the estimator are also well established [13].
2.5. Wavelet Prefilters
The Discrete Wavelet Transform (DWT) is a powerful tool in processing signals. The main issue of DWT is that the signal may be analyzed in both time and frequency domains. In identification of systems, wavelets are used mostly for nonlinear cases. In nonlinear cases, unknown time-varying coefficients are expressed as a linear combination of wavelet basis functions [14, 15]. In linear case, there are two main related papers [16, 17]. In those papers, the authors use the wavelet decomposition in order to prefilter the input/output signals. The applications related to those papers are the identification of reduced order models.
Wavelet filters are justified by their simple structure in frequency domain. They are orthogonal in time domain and can be used to specifically determine the portion of information from data without duplicity [17]. The next sections will show how wavelets are used in linear identification and the methodology proposed in this paper to identify two-time scaled systems.
2.6. System Identification considering the Frequency Domain
Considering an ARMAX model, the linear equations can be written in the following matrix form:where A is a regressors matrix, is a parameter vector, and Y is the output (measurement) vector. In the classic LS, the method tends to provide a balanced estimation over the whole frequency. In some cases, like in two-time scaled systems, it is desirable that certain frequency bands are privileged.
In order to delimit the desired frequency bands, we will transform (9) via an orthogonal transformation. Let P be an orthogonal operator; the projection of (9) onto the orthogonal basis P is given by
We can rewrite the i-th equation of (10) in the following form: where m and n are the number of regressors of y and u, respectively, and pi is the i-th row of P. We can conclude, observing (11), that each transformed equation is associated with a single eigenvector of P. The process data, after this transformation, contains only information in the subspace defined by the eigenvector pi. To control the relative effect of the information in the subspace spanned by pi, (10) may be premultiplied by a weighting matrix C.
The estimated parameters vector, after selecting the suitable operator P and the weighting matrix C, in order to fragment the spectra as desired, is obtained by classic LS expression:
The next subsection will show a possible choice of matrix P (using wavelets).
2.7. System Identification Using Wavelets
The DWT decomposes the Hilbert space (L2) into two subspaces (approximations V and details W). The resulting approximations subspace V can be decomposed further because it is closed and has a countable basis. This process may be repeated in a maximum number of levels. Denoting Vj and Wj as the space of approximations and details at the j-th level, respectively, we can derive the following expressions at the and level:or at the level:
The DWT decomposition considers a scale vector called (approximations) and a wavelet vector (details). Multiresolution analysis [18] shows us that even shifts of the same vector are orthonormal as well as shifts of different vectors. Let us define the operator shift for scale and wavelet vectors of the following form:
In practical analysis with DWT, the length of data in a level is always split in half from the previous level. Considering a data vector of length N, and j levels, we should have N=. In the approach considered by this paper (see [17]), the data of length N1=N are decomposed into approximations and details of both lengths considering the following expression:where int(x) corresponds to the integer part of x, L is the length of the filter, and S is the shift.
Considering this, we will show how to compute the matrix P defined in the previous section using wavelet filters. Analyzing the first level, the matrix P1 must consider both filters (approximations) and (details) shifted in time.
If we consider the analysis in just one level, we can see clearly that, premultiplying P1 by Y, considering that the signal Y is a real signal, we would have the data projected onto the space of approximation (lower frequencies) and details (higher frequencies). In this case, using a sampling time Ts (sampling frequency Fs), the frequency range is split in halves. In the details level, the frequency range would be fmin=fs/4 and fmax=fs/2.
In higher levels, the matrices are computed of the following form:where the dimension of Ij isand
Considering the Pj, j=1,..,p, matrices, the P matrix is computed aswhere the dimension of P, considering p levels, is given by
The weighting matrix C has some properties that we will discuss. First, the matrix C normalizes the filtered data. As we know, the gain of wavelet filters is not unitary. In this case, we can compute the matrix C in order to compensate the change of the system gain. Second, we can compute the matrix C in order to select the desired frequency ranges. As we can see in (24), the rows of P are blocks of size N2,N3,..,Nj,Nj+1. The C matrix may have the following structure:where C1 is a matrix that compensates the filter gain and C2 is a matrix defined by the user. Some proposals of its structure may be found in [16, 17]. Our approach is focused on the matrix C1 and will be shown in the next section.
2.8. Computation of the Weighting Matrix
The contribution of this paper in terms of a system identification using wavelet methodology is to provide a framework of two-time scaled systems identification that can be extended to multi-time scale systems. This extension is natural because wavelet filters split the frequency range into a given number of levels that can be computed according to the number of scales of the system. In the two-time scaled systems approach, the iterative algorithm is considered to compute the fast and the slow cut-off frequencies.
The proposed procedure is based on the structure of the matrix C1 that is given, considering p levels, by
Each identity block weights a determined level in the data decomposition (range of frequency). As we can see, it is possible to modify C1, by replacing some identity matrices by null matrices. By doing this, we disregard that frequency range (corresponding to the details of DWT). This may be very useful in the identification of two-time scaled systems. In [16, 17] the authors tune the weights according to the previous knowledge of the system. In this paper, we will consider the iterative algorithm presented in the previous Section 2.4 in order to define the number of levels considered in wavelet filtering process and which level will be considered or not by setting the identity matrices of C1 as null. To understand how to set properly C1, we must understand how the wavelet filters split the frequency range. Considering a sampled signal with sample time ( is the sample frequency), Table 1 shows the frequency range.
A more detailed explanation of the range frequencies of Table 1 may be obtained from [17]. Considering that the iterative algorithm converged, and denoting and as the cut-off frequencies of the low and high pass filters, respectively, obtained from the iterative algorithm, we will first compute the maximum number of levels from . The maximum number of levels is given bywhere round(x) is the nearest integer to .
The low-pass filter considers p levels and the matrix C1 will be set considering the last p details and approximation levels (two last rows of Table 1). The high-pass filter considers levels:
In order to tune the high-pass filter, the matrix C1 will be set considering the levels . The wavelet filter gain is usually where j is the level considered. The compensation of the wavelet filter gain is done by setting the matrix C2 with in the corresponding levels.
The algorithm that estimates the TTSS considers the wavelet prefilters. In the procedure, two wavelet matrices and two weighting matrices are computed (one for each time scale). For each time scale, a model is estimated. The pseudocode of the algorithm is shown below.
Step 1. Apply a persistent exciting signal u to the system and sample its output with sample frequency .
Step 2. Perform the iterative algorithm (see Section 2.4) and store the cut-off frequencies and .
Step 3. Compute the maximum number of the wavelet decomposition p via (27) that corresponds to the low-pass prefilter.
Step 4. Compute the maximum level of the high-pass prefilter via (28).
Step 5. Mount the regressor matrix A and the output vector Y.
Step 6. Compute the wavelet matrix P considering p levels (computed in Step 3).
Step 7. Compute the matrix considering only the levels of the wavelet decomposition (that is, set the blocks of (26), , , as null matrices).
Step 8. Compute the parameters that correspond to the fast scale model as .
Step 9. Compute matrix considering the p-th level (both details and approximation).
Step 10. Compute the parameters that correspond to the slow scale model as .
The estimated models in each time scale ( and ) considering the parameters vectors and , respectively, are the models that will be considered for validation. The overall model is obtained by (1). One important consideration that must be done is that wavelet filters are prone to noise. In order to avoid this problem, all input/output data are denoised when wavelet filters are considered.
3. Application Example and Discussion
The case study considered in this work is the identification of the longitudinal dynamics of Lockheed F104G aircraft [19]. This kind of longitudinal dynamics indeed exhibits two-time scale behaviour [3].
The excitation signal is a conventional flight test maneuver signal. Because a persistent signal is needed for identification, the flight test maneuver that has been considered is the doublet (3-2-1-1), which is a persistent excitation signal composed of square pulses. The design of doublet maneuvers consists in choosing how long the square waves will be and then choosing the signal amplitude. More details about these maneuvers are found in [20]. An example is shown in Figure 2.
3.1. Aircraft Longitudinal Dynamics
The equations of motion of the aircraft are a set of nonlinear equations (longitudinal and lateral). In our case, we will consider decoupled longitudinal and lateral dynamics, which are described by linearized equations.
An illustration of the quantities involved in the longitudinal motion is shown in Figure 3.
The linearized matrix equation of the longitudinal motion is shown in the following:where VA is the horizontal-velocity deviation in m/s; γ is the flight-path angle in rad; is the angle of attack in rad; qk is the pitch rate in rad/s; and is the elevator deflection in rad. The matrices A and B are given by
As we can see, the eigenvalues of the system matrix A are -2.4148 ±10.1386j and -0.0362 ± 0.0989j that correspond to fast and slow dynamics, respectively. In this example, we will consider VA as output (slower variable). The next subsections will show the identification results considering the methods proposed in this paper.
3.2. Results considering Classic Prefilter Tuned by the Iterative Algorithm
In this subsection, we will compare the estimations using a standard LS method and the method proposed in this paper (recursive algorithm) using prefilters. The choice of the filter structure usually is done by the user. There is a trade-off between the order of the filter and the sharpness of cut-off. In our case, a second-order filter produced good results.
In all simulations used in this paper, the system output is contaminated with a Gaussian zero mean noise. The considered sample time is =0.1s. The temporal comparison between the methods is shown in Figure 4. In this simulation, a doublet 3-2-1-1 signal produces the persistent excitation.
As we can see in Figure 4, the standard LS failed in the parameter accuracy in the estimation process. The proposed method, in time domain, estimated a model with good precision. The frequency domain validation is also presented in order to compare the identified fast and slow dynamics. The Bode plot is shown in Figure 5. The comparison in domain frequency also shows that the proposed method is more accurate mainly in high frequencies. The evolution of the cut-off frequencies for the low- and high-pass filters is shown in Figures 6 and 7.
The computational effort of the algorithm is not significant and its convergence, in this simulation, was fast (six steps). The same output signal obtained from the doublet excitation signal is used in all the iterations. When the estimation method is offline LS, the time consumption depends on the regressors matrix.
3.3. Results considering Wavelet Prefilters
The previous results from the proposed algorithm show that the cut-off frequencies are and . In this case, the maximum level of analysis, considering (27), is and the low-pass filter will be computed considering level j=10 (approximations and detail) by (28). The high-pass filter in this example will consider , so the high-pass filter is obtained from the levels .
The results, in time and frequency domains, are shown in Figures 8 and 9, respectively. An overall and quantitative comparison, considering the mean squared error (MSE), is shown in Table 2. The graphical and quantitative comparison shows that the recursive algorithm (considering classical filters) and the wavelet method produced similar results.
The time consuming in this case is more significative than the classical filter case. The product of the matrices considering P and C results in a matrix with larger dimensions. In this case, the computational effort is greater.
The models identified with the classical LS, with the iterative algorithm, and with the algorithm using wavelets are shown in (32), (33), and (34), respectively, all in discrete time.
In order to compare the fast scale and slow scale behaviour, we will present the estimated models in their respective scales.
3.4. Numerical Sensibility Analysis
In order to analyze the sensibility of the proposed method, we will present a numerical simulation considering (i) how sensitive the convergence of the iterative algorithm is in relation to its initial condition and (ii) how sensitive the filtered method (considering both classical and wavelet filters) is in relation to the cut-off frequency.
The first simulation (Figures 10 and 11) shows a deviation of ±10% in relation to the initial cut-off frequency of the iterative algorithm versus the converged cut-off frequency.
As we can see in the simulation, the algorithm demonstrated low sensibility in relation to the initial frequency. It indicates that the initial condition proposed in the algorithm was a good choice in this example.
The second numerical analysis considers a variation in the converged cut-off frequency in relation to the performance index considered in Table 2 (MSE). For both filters, a deviation of ±10% in relation to the converged algorithm is considered. A set with 50 combinations of filters into the mentioned interval were evaluated and no significant variation in MSE was observed. The simulation was performed for both the classic filter case and the wavelet filter case. The result is shown in Table 3.
As we can see, a small standard deviation was observed in the simulations. In this case, we can conclude that, in this example, a small sensibility in relation to the MSE was observed in the identifier.
4. Conclusions
This paper presents an iterative method that tunes classical prefilters in order to identify two-time scaled systems. The tuning procedure of the prefilters may be quite time consuming [17] because the user does not have enough information about the process, a priori, to accomplish this task. The results of the iterative algorithm showed that the estimated model represented with accuracy the system in whole frequency range.
The wavelet filters have a more intuitive tuning. Basically, the user must choose the maximum number of levels and choose the levels that correspond to low and high frequency. The second contribution of this paper is to take the cut-off frequencies computed by the recursive algorithm to compute the suitable maximum number of levels and the levels that correspond to low and high frequencies.
The simulation case study, which is commonly considered in aerospace industry, is a good example of TTSS and the proposed methods resulted in accurate estimated models. The classical LS method failed in estimating an accurate model as expected. Both iterative algorithm and wavelet method produced a suitable model with good accuracy in both low and high frequencies.
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research has been supported through grants 88887.092490/2015-00 (CAPES, Brazil), 309331/2015-3 (CNPq, Brazil), and EP/K03877X/1 (EPSRC, UK).