Printer Friendly

Identification of Two-Time Scaled Systems Using Prefilters.

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 e 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 e 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]:

T(s) = [T.sub.s] (s / [epsilon]) [T.sub.f] (s) (1)

where [T.sub.s] (s/[epsilon]) models the slow system dynamics scaled by the parameter e and [T.sub.f] (s) corresponds to the fast dynamics. Considering the frequency response of the system, for high frequencies we can assume that

T (j[omega]) = [T.sub.s] (j[omega]/[epsilon]) [T.sub.f] (j[omega]) [approximately equal to] [T.sub.s] (j[infinity]) [T.sub.f] (j[omega]) (2)

and for low frequencies [epsilon][omega]

T ([epsilon]j[omega]) = [T.sub.s] (j[omega]) [T.sub.f] ([epsilon]j[omega]) [approximately equal to] [T.sub.s] (j[omega]) [T.sub.f] (0) (3)

To keep these approximations, we will assume that the static gain of [T.sub.f] is nonzero and that the high-frequency gain of [T.sub.s] is neither zero nor infinity. Adopting these assumptions, we shall further consider that [T.sub.s] is biproper. From a numerical point of view, [T.sub.s] must contain both slow zeros and poles, with some of these still left in [T.sub.f]. This is to allow near cancellation in [T.sub.s] 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:

[mathematical expression not reproducible] (4)

[mathematical expression not reproducible] (5)

where [n.sub.s] is the degree of [T.sub.s]. Assumptions (4) and (5) are fulfilled by taking [F.sub.h] = 0 in low frequency and [F.sub.l] = 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 Ts dynamics, a low-pass filter is designed, and to highlight the T1 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 cutoff 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 [T.sub.s] and [T.sub.j] that will be called [T.sub.s0] and [T.sub.f0]. The cut-off frequencies of [T.sub.s0] and [T.sub.f0] 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 i [left arrow] 0.

Step 1. Apply a persistent exciting signal u to the system and sample its output.

Step 2. Estimate an initial model [T.sub.0] and compute its cut-off frequency ([f.sub.0]).

Step 3. Set [] [left arrow] [f.sub.0] and [f.sub.hi] [left arrow] [f.sub.0].

Step 4. Design a low-pass filter and a high-pass filter with cutoff frequencies [] and [f.sub.hi], respectively.

Step 5. Generate two signals [] and [] by filtering u with [] and [F.sub.hi], respectively.

Step 6. Apply the signals [] and [] to the system (separately) and sample its outputs [] and [].

Step 7. Estimate [] and [] using LS algorithm.

Step 8. Set I [left arrow] i+1.

Step 9. Compute the cut-off frequencies [] and [f.sub.hi] considering [] and [].

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:

[mathematical expression not reproducible] (6)

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).

d[l.sub.i] = [absolute value of f[l.sub.i] - f[l.sub.i-1]] (7)

d[h.sub.i] = [absolute value of f[h.sub.i] - f[h.sub.i-1]] (8)

The tolerance value e is used to define the stopping criteria; i.e., when d[l.sub.i] < [member of] and d[h.sub.i] < [member of], the iterative procedure stops. When d[l.sub.i] < [member of] and d[h.sub.i] < [member of], 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:

Y-A[theta] = 0 (9)

where A is a regressors matrix, [theta] 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

P(Y-A[theta]) = 0 (10)

We can rewrite the i-th equation of (10) in the following form:

(<[p.sup.T.sub.i], [y.sub.n]> + [a.sub.1] <[p.sup.T.sub.i], [y.sub.n-1] + ... + <[p.sup.T.sub.i], [y.sub.0]>}


-[b.sub.1] {<[p.sup.T.sub.i], [u.sub.n-1]> - ... - [b.sub.m] {<[p.sup.T.sub.i], [u.sub.n-m]>) = 0

where m and n are the number of regressors of y and u, respectively, and [p.sub.i] 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 [p.sub.i]. To control the relative effect of the information in the subspace spanned by [p.sub.i], (10) may be premultiplied by a weighting matrix C.

CPY = CPA[theta] (12)

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:

[??] = [([A.sup.T][P.sup.T]CPA).sup.-1] [A.sup.T][P.sup.T] CPY (13)

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 ([L.sup.2]) 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 [V.sub.j] and [W.sub.j] as the space of approximations and details at the j-th level, respectively, we can derive the following expressions at the 1st and 2nd level:

[L.sup.2] = [V.sub.1] [direct sum] [W.sub.1] (14)

[V.sub.1] = [V.sub.2] [direct sum] [W.sub.2] (15)

or at the pth level:

[L.sup.2] = [V.sup.p] [direct sum] [W.sup.p] [direct sum] [W.sup.p-1] [direct sum] ... [direct sum] [W.sub.1] (16)

The DWT decomposition considers a scale vector called [phi] (approximations) and a wavelet vector [psi] (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:

[q.sup.-1]z (k) = z(k - i) (17)

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=[2.sup.j]. 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:

[N.sub.j] = 1 + int ([N.sub.j-1] - L / S) (18)

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 [P.sub.1] must consider both filters [phi] (approximations) and [psi] (details) shifted in time.

[mathematical expression not reproducible] (19)

If we consider the analysis in just one level, we can see clearly that, premultiplying [P.sub.1] 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 [T.sub.s] (sampling frequency [F.sub.s]), the frequency range is split in halves. In the details level, the frequency range would be [f.sub.min] = [f.sub.s]/4 and [f.sub.max]=[f.sub.s]/2.

In higher levels, the matrices are computed of the following form:

[mathematical expression not reproducible] (20)

where the dimension of [I.sub.j] is

dim ([I.sub.j]) = [j.summation over (i=2)][N.sub.i] x [j.summation over (i=2)][N.sub.i], (21)


[mathematical expression not reproducible] (22)

Considering the [P.sub.j], j=1, ..., p, matrices, the P matrix is computed as

[mathematical expression not reproducible] (23)

where the dimension of P, considering p levels, is given by

dim (P) = ([j.summation over (i=2)][N.sub.i]) + 2[N.sub.j+1] x [N.sub.1] (24)

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 [N.sub.2],[N.sub.3], ..., [N.sub.j], [N.sub.j+1]. The C matrix may have the following structure:

C = [C.sub.2][C.sub.1] (25)

where [C.sub.1] is a matrix that compensates the filter gain and [C.sub.2] 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 [C.sub.1] 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 cutoff frequencies.

The proposed procedure is based on the structure of the matrix [C.sub.1] that is given, considering p levels, by

[mathematical expression not reproducible]. <26>

Each identity block weights a determined level in the data decomposition (range of frequency). As we can see, it is possible to modify [C.sub.1], 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 [C.sub.1] as null. To understand how to set properly [C.sub.1], we must understand how the wavelet filters split the frequency range. Considering a sampled signal with sample time [t.sub.s] ([f.sub.s] 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 [f.sub.l] and [f.sub.h] 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 [f.sub.l]. The maximum number of levels is given by

p = round ([log.sub.2] ([f.sub.s] / [f.sub.l])-1) (27)

where round(x) is the nearest integer to x.

The low-pass filter considers p levels and the matrix [C.sub.1] will be set considering the last p details and approximation levels (two last rows of Table 1). The high-pass filter considers [j.sub.h] levels:

[j.sub.h] = round ([log.sub.2] ([f.sub.s]/[f.sub.h])-1) (28)

In order to tune the high-pass filter, the matrix Cx will be set considering the levels j = 1, ..., [j.sub.h]. The wavelet filter gain is usually [([square root of (2)]).sup.j] where j is the level considered. The compensation of the wavelet filter gain is done by setting the matrix [C.sub.2] with 1 / [([square root of (2)]).sup.j] in the corresponding j = 1, ..., [j.sub.h] 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 [f.sub.s].

Step 2. Perform the iterative algorithm (see Section 2.4) and store the cut-off frequencies [f.sub.l] and [f.sub.h].

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 [j.sub.h] 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 Ch considering only the levels j = 1, ..., [j.sub.h] of the wavelet decomposition (that is, set the blocks of (26), [I.sub.j], j = [j.sub.h] + 1, ..., p, as null matrices).

Step 8. Compute the parameters that correspond to the fast scale model as [[??].sub.f] = [([A.sup.T][P.sup.T][C.sub.h]PA).sup.-1] [A.sup.T][P.sup.T][C.sub.h]PY.

Step 9. Compute matrix [C.sub.l] considering the p-th level (both details and approximation).

Step 10. Compute the parameters that correspond to the slow scale model as [[??].sub.S] = [([A.sup.T][P.sup.T][C.sub.h]PA).sup.-1] [A.sup.T][P.sup.T][C.sub.l]PY.

The estimated models in each time scale ([T.sub.s] and [T.sub.f]) considering the parameters vectors [[??].sub.s] and [[??].sub.f], 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:

[mathematical expression not reproducible] (29)

where [V.sub.A] is the horizontal-velocity deviation in m/s; [gamma] is the flight-path angle in rad; [alpha] is the angle of attack in rad; [q.sub.k] is the pitch rate in rad/s; and [eta] is the elevator deflection in rad. The matrices A and B are given by

[mathematical expression not reproducible] (30)

[mathematical expression not reproducible] (31)

As we can see, the eigenvalues of the system matrix A are -2.4148 [+ or -]10.1386/ and -0.0362 [+ or -] 0.0989/ that correspond to fast and slow dynamics, respectively. In this example, we will consider [V.sub.A] 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 [t.sub.s] = 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 [f.sub.l] = 0.0362 rad/s and [f.sub.h] = 9.9995 rad/s. In this case, the maximum level of analysis, considering (27), is p = 10 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 [j.sub.h] = 2, so the high-pass filter is obtained from the levels j = 1,2.

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.


= 0.3144[z.sup.3] + 1.003[z.sup.2] + 0.7177z + 0.1543 / [z.sup.4] - 0.8338[z.sup.3] - 0.3956[z.sup.2] - 0.05931z + 0.2926 (32)


= -0.05899[z.sup.2] + 0.5891z - 0.4589 / 4.121[z.sup.4] - 10.43[z.sup.3] + 7.176[z.sup.2] + 0.4489z - 1.319 (33)


= 1.224[z.sup.2] + 2.399z-3.014 / 5.976[z.sup.4] - 16.88[z.sup.3] + 19.38[z.sup.2] - 11.98z + 3.509 (34)

In order to compare the fast scale and slow scale behaviour, we will present the estimated models in their respective scales.

[G.sub.IAfast] = -0.0712z + 0.6504 / [z.sup.2] - 0.5375z - 0.3224 (35)

[G.sub.IAslow] = 0.8285z - 0.7056 / 4.121[z.sup.2] -8.212z + 4.091 (36)

[G.sub.wavfast] = z + 2.83 / [z.sup.2] - 0.8321z +0.591 (37)

[G.sub.wavslow] = 1-224z - 1-065 / 5.976[z.sup.2] - 11.91Z+ 5.937 (38)

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 [+ or -]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 [+ or -]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 maybe 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.


This research has been supported through grants 88887.092490/2015-00 (CAPES, Brazil), 309331/2015-3

(CNPq, Brazil), and EP/K03877X/1 (EPSRC, UK).


[1] P. Bhattacharya and M. Viceconti, "Multiscale modeling methods in biomechanics," Wiley Interdisciplinary Reviews: Systems Biology and Medicine, vol. 9, no. 3, p. e1375, 2017.

[2] L. Liu, L. Y. Wang, Z. Chen, C. Wang, F. Lin, and H. Wang, "Integrated system identification and state-of-charge estimation of battery systems," IEEE Transactions on Energy Conversion, vol. 28, no. 1, pp. 12-23, 2013.

[3] K. D. Pham, "New risk-averse control paradigm for stochastic two-time-scale systems and performance robustness," Journal of Optimization Theory and Applications, vol. 146, no. 2, pp. 511-537, 2010.

[4] P. Malisani, F. Chaplais, N. Petit, and D. Feldmann, "Thermal building model identification using time-scaled identification methods," in Proceedings of the 2010 49th IEEE Conference on Decision and Control, CDC 2010, pp. 308-315, 2010.

[5] D. W. Luse and H. K. Khalil, "Frequency domain results for systems with slow and fast dynamics," Institute of Electrical and Electronics Engineers Transactions on Automatic Control, vol. 30, no. 12, pp. 1171-1179,1985.

[6] V. R. Saksena, J. O'Reilly, and P. V. Kokotovic, "Singular perturbations and time-scale methods in control theory: survey 1976-1983," Automatica, vol. 20, no. 3, pp. 273-293, 1984.

[7] J. Li, W. Ge, J. Zhang, and M. Kwauk, "Multi-scale compromise and multi-level correlation in complex systems," Chemical Engineering Research and Design, vol. 83, no. A6, pp. 574-582, 2005.

[8] F. Chaplais and K. Alaoui, "Two time scaled parameter identification by coordination of local identifiers," Automatica, vol. 32, no. 9, pp. 1303-1309,1996.

[9] M. Sitterly, L. Y. Wang, G. G. Yin, and C. Wang, "Enhanced identification of battery models for real-time battery management," IEEE Transactions on Sustainable Energy, vol. 2, no. 3, pp. 300308, 2009.

[10] Y. Hu and Y.-Y. Wang, "Two time-scaled battery model identification with application to battery state estimation," IEEE Transactions on Control Systems Technology, vol. 23, no. 2, pp. 1180-1188, 2015.

[11] R. H. Middleton and G. C. Goodwin, Digital Control and Estimation, a Unified Approach, Prentice-Hall, Englewood Cliffs, NJ, USA, 1990.

[12] I. D. Landau, Identification et Commande des Systems, Hermes, Paris, France, 1988.

[13] L. Ljung, System Identification: Theory for the User, Prentice-Hall, Englewood Cliffs, NJ, USA, 1987

[14] M. K. Tsatsanis and G. B. Giannakis, "Time-Varying System Identification and Model Validation Using Wavelets," IEEE Transactions on Signal Processing, vol. 41, no. 12, pp. 3512-3523, 2002.

[15] H.-L. Wei and S. A. Billings, "Identification of time-varying systems using multiresolution wavelet models," International Journal of Systems Science, vol. 33, no. 15, pp. 1217-1228, 2002.

[16] J. F. Carrier and G. Stephanopoulos, "Wavelet-Based Modulation in Control-Relevant Process Identification," AIChE Journal, vol. 44, no. 2, pp. 341-360, 1998.

[17] Z. Vam and H. A. Preisig, "System identification in frequency domain using wavelets: Conceptual remarks," Systems & Control Letters, vol. 61, pp. 1041-1051, 2012.

[18] C. S. Burrus, R. A. Gopinath, and H. Guo, Introduction to Wavelets and Wavelet Transforms, Prentice Hall, New Jersey, NJ, USA, 1998.

[19] R. Brockhaus, Flugregelung, Springer, Berlin, Germany, 1994.

[20] N . S . B. Neto, B. C . O. Maciel, P. H . I . A . Oliveira, E . M . Hemerly, and L. C. S. Goes, "Comparison Between Conventional Flight Test Maneuvers for Aircraft Aerodynamic Derivatives Estimation," in Proceedings of the XII Brazilian Congress of Mechanical Engineering, 2005.

Anderson L. O. Cavalcanti, (1) Karl H. Kienitz, (2) and Visakan Kadirkamanathan (3)

(1) Department of Computation and Automation, Federal University of Rio Grande do Norte, Natal, RN, Brazil

(2) Division of Electronic Engineering, Technological Institute of Aeronautics, Sao Jose dos Campos, SP, Brazil

(3) Department of Automatic Control & Systems Engineering, University of Sheffield, Sheffield, UK

Correspondence should be addressed to Anderson L. O. Cavalcanti;

Received 15 April 2018; Revised 24 July 2018; Accepted 17 September 2018; Published 3 October 2018

Academic Editor: Mario Russo

Caption: FIGURE 1: Identification using prefilters.

Caption: FIGURE 2: Doublet (3-2-1-1) persistent excitation signal.

Caption: FIGURE 3: Aircraft illustration (longitudinal motion).

Caption: FIGURE 4: Comparison between identified models: horizontal velocity deviation (m.[s.sup.-1]).

Caption: FIGURE 5: Model identified on full range frequency: Bode plot.

Caption: FIGURE 6: Evolution of cut-off frequencies: low-pass filter.

Caption: FIGURE 7: Evolution of cut-off frequencies: high-pass filter.

Caption: FIGURE 8: Comparison between identified models: horizontal-velocity deviation (m.[s.sup.-1]).

Caption: FIGURE 9: Model identified on full range frequency: Bode plot.

Caption: FIGURE 10: Initial cut-off frequency of the iterative algorithm versus the converged cut-off frequency: low-pass filter.

Caption: FIGURE 11: Initial cut-off frequency of the iterative algorithm versus the converged cut-off frequency: high-pass filter.
TABLE 1: Frequency range in levels of wavelet decomposition.

Level     Analysis              fmin                    fmax

1          Details           [f.sub.s]/4             [f.sub.s]/2
2          Details           [f.sub.s]/8             [f.sub.s]/4
3          Details          [f.sub.s]/16             [f.sub.s]/8
...          ...                 ...                     ...
pth        Details      [f.sub.s]/[2.sup.p+1]    [f.sub.s]/[2.sup.P]
pth     Approximation             0             [f.sub.s]/[2.sup.P+1]

TABLE 2: Comparison of the mean squared error (MSE).

Method                  MSE

Standard LS           276.0001
Recursive Algorithm   42.8698
Wavelet               46.5788

TABLE 3: Sensibility analysis considering a deviation in the
cut-off frequency.

Method                Mean of MSE   Standard Deviation of MSE

Recursive Algorithm     43.1935              0.9911
Wavelet                 47.3041               1.023
COPYRIGHT 2018 Hindawi Limited
No portion of this article can be reproduced without the express written permission from the copyright holder.
Copyright 2018 Gale, Cengage Learning. All rights reserved.

Article Details
Printer friendly Cite/link Email Feedback
Title Annotation:Research Article
Author:Cavalcanti, Anderson L.O.; Kienitz, Karl H.; Kadirkamanathan, Visakan
Publication:Journal of Control Science and Engineering
Article Type:Report
Date:Jan 1, 2018
Previous Article:Multiple Feature Vectors Based Fault Classification for WSN Integrated Bearing of Rolling Mill.
Next Article:Research on Integrated Guidance and Control of Distributed Cooperation of Multi-Interceptor with State Coupling.

Terms of use | Privacy policy | Copyright © 2021 Farlex, Inc. | Feedback | For webmasters |