# Identification of Two-Time Scaled Systems Using Prefilters

Authors: Anderson L. O. Cavalcanti; Karl H. Kienitz; Visakan Kadirkamanathan
Journal: Journal of Control Science and Engineering
 (2018)

---

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


---

## Body

## 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]:(1)Ts=TssεTfswhere Tss/ε models the slow system dynamics scaled by the parameter ε and Tfs corresponds to the fast dynamics. Considering the frequency response of the system, for high frequencies we can assume that(2)Tjω=TsjωεTfjω≈Tsj∞Tfjωand for low frequencies εω(3)Tεjω=TsjωTfεjω≈TsjωTf0To keep these approximations, we will assume that the static gain ofTf is nonzero and that the high-frequency gain ofTs is neither zero nor infinity. Adopting these assumptions, we shall further consider thatTs is biproper. From a numerical point of view,Ts must contain both slow zeros and poles, with some of these still left inTf. This is to allow near cancellation inTs 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 Figure1.Figure 1
Identification using prefilters.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:(4)limε→0⁡maxω∈R⁡Tsjωε-Tsj∞Fhjω=0(5)limε→0⁡maxω∈R⁡Tfεjω-Tf0ωkFljω=0,0≤k≤nswherens is the degree ofTs. Assumptions (4) and (5) are fulfilled by takingFh = 0 in low frequency andFl = 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 Tf 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 ofTs andTf that will be calledTs0 andTf0. The cut-off frequencies ofTs0 andTf0 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⟵0.Step 1.
Apply a persistent exciting signalu to the system and sample its output.Step 2.
Estimate an initial modelT0 and compute its cut-off frequency (f0).Step 3.
Setfli⟵f0 and fhi⟵f0.Step 4.
Design a low-pass filter and a high-pass filter with cut-off frequenciesfli andfhi, respectively.Step 5.
Generate two signalsusi andufi by filteringu withFli andFhi, respectively.Step 6.
Apply the signalsusi andusi to the system (separately) and sample its outputsysi andyfi.Step 7.
EstimateTsi andTfi using LS algorithm.Step 8.
Set i⟵i+1.Step 9.
Compute the cut-off frequenciesfli andfhi consideringTsi-1 andTfi-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:(6)ωcutoff=-max⁡ω⁡realpolesGjω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).(7)dli=fli-fli-1(8)dhi=fhi-fhi-1The tolerance valueϵ is used to define the stopping criteria; i.e., whendli< ϵ anddhi <ϵ, the iterative procedure stops. Whendli< ϵ anddhi < ϵ, 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:(9)Y-Aθ=0whereA is a regressors matrix, θ is a parameter vector, andY 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. LetP be an orthogonal operator; the projection of (9) onto the orthogonal basisP is given by(10)PY-Aθ=0We can rewrite the i-th equation of (10) in the following form: (11)piT,yn+a1piT,yn-1+…+piT,y0-b1piT,un-1-…-bmpiT,un-m=0wherem andn are the number of regressors ofyandu, respectively, andpi is the i-th row ofP. We can conclude, observing (11), that each transformed equation is associated with a single eigenvector ofP. The process data, after this transformation, contains only information in the subspace defined by the eigenvectorpi. To control the relative effect of the information in the subspace spanned bypi, (10) may be premultiplied by a weighting matrixC.(12)CPY=CPAθThe estimated parameters vector, after selecting the suitable operatorP and the weighting matrixC, in order to fragment the spectra as desired, is obtained by classic LS expression:(13)θ^=ATPTCPA-1ATPTCPYThe next subsection will show a possible choice of matrixP (using wavelets).


### 2.7. System Identification Using Wavelets

The DWT decomposes the Hilbert space (L2) into two subspaces (approximationsV and detailsW). The resulting approximations subspaceV can be decomposed further because it is closed and has a countable basis. This process may be repeated in a maximum number of levels. DenotingVj andWj 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:(14)L2=V1⊕W1(15)V1=V2⊕W2or at the pth level:(16)L2=Vp⊕Wp⊕Wp-1⊕…⊕W1The 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:(17)q-izk=zk-iIn 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 lengthN, andj levels, we should haveN=2j. In the approach considered by this paper (see [17]), the data of lengthN1=N are decomposed into approximations and details of both lengths considering the following expression:(18)Nj=1+intNj-1-LSwhere int(x) corresponds to the integer part ofx,L is the length of the filter, andS is the shift.Considering this, we will show how to compute the matrixP defined in the previous section using wavelet filters. Analyzing the first level, the matrixP1 must consider both filters φ (approximations) and ψ (details) shifted in time.(19)P1=ψq-2ψ⋮q-2N2ψ---φq-2φ⋮q-2N2φ2N2×N1If we consider the analysis in just one level, we can see clearly that, premultiplyingP1 byY, considering that the signalYis 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 timeTs (sampling frequencyFs), the frequency range is split in halves. In the details level, the frequency range would befmin=fs/4andfmax=fs/2.In higher levels, the matrices are computed of the following form:(20)Pj=Ij00Pj-∑i=2jNi+2Ni+1×∑i=2j-1Ni+2Niwhere the dimension ofIj is(21)dimIj=∑i=2jNi×∑i=2jNiand(22)Pj-=ψq-2ψ⋮q-2Nj+1-1ψ---φq-2φ⋮q-2Nj+1-1φ2Nj×Nj-1Considering thePj,j=1,..,p, matrices, theP matrix is computed as(23)P=∏i=p1Piwhere the dimension ofP, consideringp levels, is given by (24)dim⁡P=∑i=2jNi+2Nj+1×N1The weighting matrixC has some properties that we will discuss. First, the matrixCnormalizes the filtered data. As we know, the gain of wavelet filters is not unitary. In this case, we can compute the matrixC in order to compensate the change of the system gain. Second, we can compute the matrixC in order to select the desired frequency ranges. As we can see in (24), the rows ofP are blocks of sizeN2,N3,..,Nj,Nj+1. TheC matrix may have the following structure:(25)C=C2C1whereC1 is a matrix that compensates the filter gain andC2 is a matrix defined by the user. Some proposals of its structure may be found in [16, 17]. Our approach is focused on the matrixC1 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 matrixC1 that is given, consideringp levels, by(26)C1=IN200⋱…0⋯00000INp+1⋮⋯INp+1Each identity block weights a determined level in the data decomposition (range of frequency). As we can see, it is possible to modifyC1, 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 ofC1 as null. To understand how to set properlyC1, we must understand how the wavelet filters split the frequency range. Considering a sampled signal with sample time ts (fs is the sample frequency), Table 1 shows the frequency range.Table 1
Frequency range in levels of wavelet decomposition.
Level Analysis fmin fmax 1 Details f s/4 f s/2  2 Details f s/8 f s/4  3 Details f s/16 f s/8  ⋮ ⋮ ⋮ ⋮  p t h Details f s / 2 p + 1 f s / 2 p  p t h Approximation 0 f s / 2 p + 1A more detailed explanation of the range frequencies of Table1 may be obtained from [17]. Considering that the iterative algorithm converged, and denoting fl and fh 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 fl. The maximum number of levels is given by(27)p=roundlog2⁡fsfl-1whereround(x) is the nearest integer to x.The low-pass filter considersp levels and the matrixC1 will be set considering the lastp details and approximation levels (two last rows of Table 1). The high-pass filter considers jh levels:(28)jh=roundlog2⁡fsfh-1In order to tune the high-pass filter, the matrixC1 will be set considering the levels j=1,⋯,jh. The wavelet filter gain is usually 2j wherej is the level considered. The compensation of the wavelet filter gain is done by setting the matrixC2 with 1/2j in the corresponding j=1,⋯,jh 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 frequencyfs.Step 2.
Perform the iterative algorithm (see Section2.4) and store the cut-off frequencies fl and fh.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 prefilterjh via (28).Step 5.
Mount the regressor matrixA and the output vector Y.Step 6.
Compute the wavelet matrixP consideringp levels (computed in Step 3).Step 7.
Compute the matrixCh considering only the levels j=1,⋯,jh of the wavelet decomposition (that is, set the blocks of (26), Ij, j=jh+1,⋯,p, as null matrices).Step 8.
Compute the parameters that correspond to the fast scale model asθ^f=(ATPTChPA)-1ATPTChPY.Step 9.
Compute matrixCl considering the p-th level (both details and approximation).Step 10.
Compute the parameters that correspond to the slow scale model asθ^s=(ATPTChPA)-1ATPTClPY.The estimated models in each time scale (Ts and Tf) considering the parameters vectors θ^s and θ^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.


## 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]:(1)Ts=TssεTfswhere Tss/ε models the slow system dynamics scaled by the parameter ε and Tfs corresponds to the fast dynamics. Considering the frequency response of the system, for high frequencies we can assume that(2)Tjω=TsjωεTfjω≈Tsj∞Tfjωand for low frequencies εω(3)Tεjω=TsjωTfεjω≈TsjωTf0To keep these approximations, we will assume that the static gain ofTf is nonzero and that the high-frequency gain ofTs is neither zero nor infinity. Adopting these assumptions, we shall further consider thatTs is biproper. From a numerical point of view,Ts must contain both slow zeros and poles, with some of these still left inTf. This is to allow near cancellation inTs 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 Figure1.Figure 1
Identification using prefilters.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:(4)limε→0⁡maxω∈R⁡Tsjωε-Tsj∞Fhjω=0(5)limε→0⁡maxω∈R⁡Tfεjω-Tf0ωkFljω=0,0≤k≤nswherens is the degree ofTs. Assumptions (4) and (5) are fulfilled by takingFh = 0 in low frequency andFl = 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 Tf 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 ofTs andTf that will be calledTs0 andTf0. The cut-off frequencies ofTs0 andTf0 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⟵0.Step 1.
Apply a persistent exciting signalu to the system and sample its output.Step 2.
Estimate an initial modelT0 and compute its cut-off frequency (f0).Step 3.
Setfli⟵f0 and fhi⟵f0.Step 4.
Design a low-pass filter and a high-pass filter with cut-off frequenciesfli andfhi, respectively.Step 5.
Generate two signalsusi andufi by filteringu withFli andFhi, respectively.Step 6.
Apply the signalsusi andusi to the system (separately) and sample its outputsysi andyfi.Step 7.
EstimateTsi andTfi using LS algorithm.Step 8.
Set i⟵i+1.Step 9.
Compute the cut-off frequenciesfli andfhi consideringTsi-1 andTfi-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:(6)ωcutoff=-max⁡ω⁡realpolesGjω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).(7)dli=fli-fli-1(8)dhi=fhi-fhi-1The tolerance valueϵ is used to define the stopping criteria; i.e., whendli< ϵ anddhi <ϵ, the iterative procedure stops. Whendli< ϵ anddhi < ϵ, 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:(9)Y-Aθ=0whereA is a regressors matrix, θ is a parameter vector, andY 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. LetP be an orthogonal operator; the projection of (9) onto the orthogonal basisP is given by(10)PY-Aθ=0We can rewrite the i-th equation of (10) in the following form: (11)piT,yn+a1piT,yn-1+…+piT,y0-b1piT,un-1-…-bmpiT,un-m=0wherem andn are the number of regressors ofyandu, respectively, andpi is the i-th row ofP. We can conclude, observing (11), that each transformed equation is associated with a single eigenvector ofP. The process data, after this transformation, contains only information in the subspace defined by the eigenvectorpi. To control the relative effect of the information in the subspace spanned bypi, (10) may be premultiplied by a weighting matrixC.(12)CPY=CPAθThe estimated parameters vector, after selecting the suitable operatorP and the weighting matrixC, in order to fragment the spectra as desired, is obtained by classic LS expression:(13)θ^=ATPTCPA-1ATPTCPYThe next subsection will show a possible choice of matrixP (using wavelets).


## 2.7. System Identification Using Wavelets

The DWT decomposes the Hilbert space (L2) into two subspaces (approximationsV and detailsW). The resulting approximations subspaceV can be decomposed further because it is closed and has a countable basis. This process may be repeated in a maximum number of levels. DenotingVj andWj 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:(14)L2=V1⊕W1(15)V1=V2⊕W2or at the pth level:(16)L2=Vp⊕Wp⊕Wp-1⊕…⊕W1The 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:(17)q-izk=zk-iIn 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 lengthN, andj levels, we should haveN=2j. In the approach considered by this paper (see [17]), the data of lengthN1=N are decomposed into approximations and details of both lengths considering the following expression:(18)Nj=1+intNj-1-LSwhere int(x) corresponds to the integer part ofx,L is the length of the filter, andS is the shift.Considering this, we will show how to compute the matrixP defined in the previous section using wavelet filters. Analyzing the first level, the matrixP1 must consider both filters φ (approximations) and ψ (details) shifted in time.(19)P1=ψq-2ψ⋮q-2N2ψ---φq-2φ⋮q-2N2φ2N2×N1If we consider the analysis in just one level, we can see clearly that, premultiplyingP1 byY, considering that the signalYis 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 timeTs (sampling frequencyFs), the frequency range is split in halves. In the details level, the frequency range would befmin=fs/4andfmax=fs/2.In higher levels, the matrices are computed of the following form:(20)Pj=Ij00Pj-∑i=2jNi+2Ni+1×∑i=2j-1Ni+2Niwhere the dimension ofIj is(21)dimIj=∑i=2jNi×∑i=2jNiand(22)Pj-=ψq-2ψ⋮q-2Nj+1-1ψ---φq-2φ⋮q-2Nj+1-1φ2Nj×Nj-1Considering thePj,j=1,..,p, matrices, theP matrix is computed as(23)P=∏i=p1Piwhere the dimension ofP, consideringp levels, is given by (24)dim⁡P=∑i=2jNi+2Nj+1×N1The weighting matrixC has some properties that we will discuss. First, the matrixCnormalizes the filtered data. As we know, the gain of wavelet filters is not unitary. In this case, we can compute the matrixC in order to compensate the change of the system gain. Second, we can compute the matrixC in order to select the desired frequency ranges. As we can see in (24), the rows ofP are blocks of sizeN2,N3,..,Nj,Nj+1. TheC matrix may have the following structure:(25)C=C2C1whereC1 is a matrix that compensates the filter gain andC2 is a matrix defined by the user. Some proposals of its structure may be found in [16, 17]. Our approach is focused on the matrixC1 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 matrixC1 that is given, consideringp levels, by(26)C1=IN200⋱…0⋯00000INp+1⋮⋯INp+1Each identity block weights a determined level in the data decomposition (range of frequency). As we can see, it is possible to modifyC1, 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 ofC1 as null. To understand how to set properlyC1, we must understand how the wavelet filters split the frequency range. Considering a sampled signal with sample time ts (fs is the sample frequency), Table 1 shows the frequency range.Table 1
Frequency range in levels of wavelet decomposition.
Level Analysis fmin fmax 1 Details f s/4 f s/2  2 Details f s/8 f s/4  3 Details f s/16 f s/8  ⋮ ⋮ ⋮ ⋮  p t h Details f s / 2 p + 1 f s / 2 p  p t h Approximation 0 f s / 2 p + 1A more detailed explanation of the range frequencies of Table1 may be obtained from [17]. Considering that the iterative algorithm converged, and denoting fl and fh 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 fl. The maximum number of levels is given by(27)p=roundlog2⁡fsfl-1whereround(x) is the nearest integer to x.The low-pass filter considersp levels and the matrixC1 will be set considering the lastp details and approximation levels (two last rows of Table 1). The high-pass filter considers jh levels:(28)jh=roundlog2⁡fsfh-1In order to tune the high-pass filter, the matrixC1 will be set considering the levels j=1,⋯,jh. The wavelet filter gain is usually 2j wherej is the level considered. The compensation of the wavelet filter gain is done by setting the matrixC2 with 1/2j in the corresponding j=1,⋯,jh 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 frequencyfs.Step 2.
Perform the iterative algorithm (see Section2.4) and store the cut-off frequencies fl and fh.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 prefilterjh via (28).Step 5.
Mount the regressor matrixA and the output vector Y.Step 6.
Compute the wavelet matrixP consideringp levels (computed in Step 3).Step 7.
Compute the matrixCh considering only the levels j=1,⋯,jh of the wavelet decomposition (that is, set the blocks of (26), Ij, j=jh+1,⋯,p, as null matrices).Step 8.
Compute the parameters that correspond to the fast scale model asθ^f=(ATPTChPA)-1ATPTChPY.Step 9.
Compute matrixCl considering the p-th level (both details and approximation).Step 10.
Compute the parameters that correspond to the slow scale model asθ^s=(ATPTChPA)-1ATPTClPY.The estimated models in each time scale (Ts and Tf) considering the parameters vectors θ^s and θ^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.Figure 2
Doublet (3-2-1-1) persistent excitation signal.


### 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 Figure3.Figure 3
Aircraft illustration (longitudinal motion).The linearized matrix equation of the longitudinal motion is shown in the following:(29)q˙Kα˙V˙Aγ˙=AqKαVAγ+B.ηwhereVA is the horizontal-velocity deviation inm/s; γ is the flight-path angle inrad; α is the angle of attack in rad;qk is the pitch rate inrad/s; and η is the elevator deflection inrad. The matricesA andB are given by(30)A=-1.980-103.00001.000-2.85.00000.00700.0000-0.00010.00000.000-17.20700.0002.8500.0720-9.80700.00010.0000(31)B=-60.2-0.339300.3393As 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 considerVA 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 ists=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.Figure 4
Comparison between identified models: horizontal velocity deviation (m.s−1).As we can see in Figure4, 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.Figure 5
Model identified on full range frequency: Bode plot.Figure 6
Evolution of cut-off frequencies: low-pass filter.Figure 7
Evolution of cut-off frequencies: high-pass filter.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 arefl=0.0362rad/s and fh=9.9995rad/s. In this case, the maximum level of analysis, considering (27), is p=10 and the low-pass filter will be computed considering levelj=10 (approximations and detail) by (28). The high-pass filter in this example will consider jh=2, so the high-pass filter is obtained from the levels j=1,2.The results, in time and frequency domains, are shown in Figures8 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.Table 2
Comparison of the mean squared error (MSE).
Method MSE Standard LS 276.0001  Recursive Algorithm 42.8698  Wavelet 46.5788Figure 8
Comparison between identified models: horizontal-velocity deviation (m.s−1).Figure 9
Model identified on full range frequency: Bode plot.The time consuming in this case is more significative than the classical filter case. The product of the matrices consideringP andC 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. (32)GLS=0.3144z3+1.003z2+0.7177z+0.1543z4-0.8338z3-0.3956z2-0.05931z+0.2926(33)GIA=-0.05899z2+0.5891z-0.45894.121z4-10.43z3+7.176z2+0.4489z-1.319(34)Gwav=1.224z2+2.399z-3.0145.976z4-16.88z3+19.38z2-11.98z+3.509In order to compare the fast scale and slow scale behaviour, we will present the estimated models in their respective scales.(35)GIAfast=-0.0712z+0.6504z2-0.5375z-0.3224(36)GIAslow=0.8285z-0.70564.121z2-8.212z+4.091(37)Gwavfast=z+2.83z2-0.8321z+0.591(38)Gwavslow=1.224z-1.0655.976z2-11.91z+5.937


### 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 (Figures10 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.Figure 10
Initial cut-off frequency of the iterative algorithm versus the converged cut-off frequency: low-pass filter.Figure 11
Initial cut-off frequency of the iterative algorithm versus the converged cut-off frequency: high-pass filter.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 Table2 (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.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.023As 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.


## 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 Figure3.Figure 3
Aircraft illustration (longitudinal motion).The linearized matrix equation of the longitudinal motion is shown in the following:(29)q˙Kα˙V˙Aγ˙=AqKαVAγ+B.ηwhereVA is the horizontal-velocity deviation inm/s; γ is the flight-path angle inrad; α is the angle of attack in rad;qk is the pitch rate inrad/s; and η is the elevator deflection inrad. The matricesA andB are given by(30)A=-1.980-103.00001.000-2.85.00000.00700.0000-0.00010.00000.000-17.20700.0002.8500.0720-9.80700.00010.0000(31)B=-60.2-0.339300.3393As 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 considerVA 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 ists=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.Figure 4
Comparison between identified models: horizontal velocity deviation (m.s−1).As we can see in Figure4, 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.Figure 5
Model identified on full range frequency: Bode plot.Figure 6
Evolution of cut-off frequencies: low-pass filter.Figure 7
Evolution of cut-off frequencies: high-pass filter.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 arefl=0.0362rad/s and fh=9.9995rad/s. In this case, the maximum level of analysis, considering (27), is p=10 and the low-pass filter will be computed considering levelj=10 (approximations and detail) by (28). The high-pass filter in this example will consider jh=2, so the high-pass filter is obtained from the levels j=1,2.The results, in time and frequency domains, are shown in Figures8 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.Table 2
Comparison of the mean squared error (MSE).
Method MSE Standard LS 276.0001  Recursive Algorithm 42.8698  Wavelet 46.5788Figure 8
Comparison between identified models: horizontal-velocity deviation (m.s−1).Figure 9
Model identified on full range frequency: Bode plot.The time consuming in this case is more significative than the classical filter case. The product of the matrices consideringP andC 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. (32)GLS=0.3144z3+1.003z2+0.7177z+0.1543z4-0.8338z3-0.3956z2-0.05931z+0.2926(33)GIA=-0.05899z2+0.5891z-0.45894.121z4-10.43z3+7.176z2+0.4489z-1.319(34)Gwav=1.224z2+2.399z-3.0145.976z4-16.88z3+19.38z2-11.98z+3.509In order to compare the fast scale and slow scale behaviour, we will present the estimated models in their respective scales.(35)GIAfast=-0.0712z+0.6504z2-0.5375z-0.3224(36)GIAslow=0.8285z-0.70564.121z2-8.212z+4.091(37)Gwavfast=z+2.83z2-0.8321z+0.591(38)Gwavslow=1.224z-1.0655.976z2-11.91z+5.937


## 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 (Figures10 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.Figure 10
Initial cut-off frequency of the iterative algorithm versus the converged cut-off frequency: low-pass filter.Figure 11
Initial cut-off frequency of the iterative algorithm versus the converged cut-off frequency: high-pass filter.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 Table2 (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.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.023As 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.



---

*Source: 3138149-2018-10-03.xml*