Time-Invariant Dynamic Systems Identification Based on ... - CiteSeerX
Time-Invariant Dynamic Systems Identification Based on the Qualitative Features of the Response
Juan J. Flores a , Nelio Pastor a,∗ a Divisi´ on
de Estudios de Postgrado, Facultad de Ingenier´ıa El´ectrica, Universidad Michoacana de San Nicol´as de Hidalgo, Morelia, M´exico
Abstract The problem of Systems Identification starts with a time-series of observed data and tries to determine the simplest model capable of exhibiting the observed behavior. This optimization problem searches the model from a space of possible models. In this paper, we present the theory and algorithms to perform Qualitative and Quantitative Systems Identification for Linear Time Invariant Dynamic Systems. The methods described here are based on successive elimination of the components of the system’s response. Sinusoidals of high frequencies are eliminated first, then their carrying waves. We continue with the process until we obtain a non-oscillatory carrier. At that point, we determine the order of the carrier. This procedure allows us to determine how many sinusoidal components, and how many exponential components are found in the impulse response of the system under study. The number of components determines the order of the system. The paper is composed of two important parts, the statement of some mathematical properties of the responses of Linear Time Invariant Dynamic Systems, and the proposal of a set of filters that allows us to implement the recognition algorithm. We present the aplication of the proposed methodology to analysis and modeling of electrical circuits and electrical power systems.
Key words: Dynamic systems, systems identification, differential equations, qualitative reasoning
This research has been supported by CIC-UMSNH grant No. 10.3 ∗ Correspondig author. Address: Maestro de Avila 90, Fracc. Fray Antonio de San Miguel, Morelia, Michoac´an, 58270, M´exico Email addresses: [email protected] (Juan J. Flores), [email protected] (Nelio Pastor).
Preprint submitted to Engeenering Applications of Artificial Intelligence 10 March 2004
1 Introduction
System identification deals with the problem of building mathematical models of dynamical systems based on observed data from the system and tries to determine the simplest model capable of exhibiting the observed behavior. This optimization problem searches the model from a space of possible models. In traditional methods, the process of structural identification has received less attention than the parametrization of the model. In most of cases the structural estimation is not generally made. Typically, a certain model structure is chosen by the user which contains unknown parameters. The choice of an appropiate model structure is most crucial for a successful identification application. This choice must be based both on an understanding of the identification procedure and on insights and knowledge about the system to be identified. An alternating way to infer a suitable structure, guided by system knowledge and the collected data set, is presented in this paper. In our approach we limit the search space to Linear Time-Invariant models.
The problem we are assessing in this paper is that of Qualitative Systems Identification for Linear Time Invariant (LTI) Dynamic Systems. This kind of systems can be represented by Linear ordinary Differential Equations with Constant Coefficients. Although this problem may seem over-constrained, there are many important problems in engineering and physics that can be expressed in mathematical terms by this kind of differential equations. We could even assess time-varying or even non-linear systems, if we consider them as piece-wise decomposed by linear approximations of the original systems. So, this is a limited yet interesting domain to start with.
The Systems Identification process can be decomposed in two steps: the first step called structural (or qualitative) identification, involves determining the qualitative features of the mechanism, i.e. the qualitative form of the system inside the black box [Bradley, Easley, & Stolle2001] [Bellazzi, Guglielmann, & Ironi1999] [Bradley & Stolle1996] [Kay, Rinner, & Kuipers2000]; once we know the nature of the mechanism, in the second step we proceed to determine the numerical value of the parameters of the model determined in the first step.
This second part can be done by any optimization process (e.g. minimum square error [Ljung1987] [P.1989], genetic algorithms [Goldberg1998] [Haupt & Haupt1998] [Hunt1993] [Pastor2000] [Zhang Zibo1987], etc.). Some algorithms use an optimization process to determine both faces at the same time (e.g. genetic algorithms [Downing et al.1996] [Kristinsson & Dumont1992]). 2
2 Linear Time-Invariant Dynamic Systems
Linear ODEs with constant coefficients are the most studied kind of differential equations; they have complete analytical solutions. So, if we restrict the kind of systems we are to analyze to those that can be expressed by an nth-order ordinary differential equation with constant coefficients, we know the kind of responses we are to get. We can express the behavior of a system in terms of the exponential and sinusoidal components in the response. We define En1 (t) =
1≤i≤n1
ai eri t
(1)
as a summation with at most n exponential terms, and ESn2 (t) =
1≤i≤n2
ai eri t sin ωi t
(2)
as a summation of exponentially decreasing sinusoidal functions. Note that we are not interested in giving analytical solutions to the differential equation, but a qualitative description of its behaviors. Theorem 1 The response of an nth order system can be expressed as in Equation 3. X(t) = En1 (t) + ESn2 (t)
(3)
where n1 + 2n2 = n. This result is evident from the definitions of Equations( 1) and (2). We see that if the second term of Equation 3 does not exist, the response will be acyclic. Otherwise, it is a sinusoidal wave, where En1 (t) represents its axis or attractor, and ESn2 (t) its sinusoidal components. Note that if we include a forcing function, the system’s response would be decomposed into Natural (the solution to the homogeneous equation) and Forced responses. If we restrict the forcing functions to be of the form e αt sin βt (i.e. constant, exponential, or sinusoidal), the forced response always has the same qualitative form as the forcing function [Boyce & DiPrima1969]. The force response would preserve the qualitative form of the response, and only adds one more exponential or sinusoidal term to the response. Let us analyze the qualitative form of the responses, as expressed by Equation (3). This qualitative form can be derived from the qualitative form of its exponential and sinusoidal components. 3
2.1 Exponential Components
The qualitative behavior of the exponential part of the response is characterized by Theorem 2.
Theorem 2 X(t) = En (t) = 1≤i≤n ai eri t exhibits at most n extrema (maxima or minima), including the one when t → ∞. We will assume, without loss of generality, that r 1 > r2 > . . . > rn . Theorem 2 is equivalent to saying that the derivative X (t) has at most n different zeroes. That is, X (t) = −r1 a1 e−r1 t − . . . − rn an e−rn t = 0 = r1 a1 e−r1 t + . . . + rn an e−rn t = 0
(4)
Performing the variable change z = e−t , Equation (4) becomes −X (t) = r1 a1 z r1 + . . . + rn an z rn = 0
(5)
Equation (5) is a polynomial in z of degree r 1 . Based on Descartes’ theorem [Kurosch1977,Rees, Sparks, & Re the number of roots of a polynomial are determined by the sign changes in the polynomial coefficients. Since we have n terms in polynomial (4), we can have at most n − 1 changes of sign, and therefore n − 1 roots. There is one more root placed where the function becomes zero, and that is when t → ∞. Given we have at most n roots for Equation (4), we can therefore have at most n extrema for En (t).
2.2 Sinusoidal Components
If the frequencies of two sinusoidal components of Equation (3) are equal, their shapes are reduced to one [Flores & Farley1995]. If their frequencies are different, they can be seen as the faster sinusoidal mounted on the slower one. If two frequencies of the sinusoidal components are very close together, the resulting wave beats. The results presented in this section fully characterizes all possible responses of a LTI Dynamic System. In the next section we will describe how to use these results to produce a framework for performing Systems Identification at the qualitative and quantitative levels. 4
3 QSI Algorithm
The identification algorithm presented in this section is based on the fact that the response of a LTI system can be decomposed in a sumation of exponential terms. If some of those exponential terms are complex, in which case they are conjugate complex pairs, each pair forms a sinusoidal. The algorithm that we propose in this paper takes as input a time-series representing the response of an LTI dynamic system and is capable of separating the terms of Equation (3) to determine the structure or qualitative form of the system exhibiting the observed behavior. Separating the terms of the system’s response is performed by a filtering process. Figure (1) shows a simplified version of the QSI algorithm. QSI determines the order of the system by adding the order of all eliminated components. This algorithm analyses the response and eliminates each component at a time, there are two main functions in QSI algorithm: T P AF ilter and ExpF ilter. The function T P AF ilter (see next section) eliminates the sinusoidal components and returns the order corresponding to those components, the remainder signal and the parameters. The function Expf ilter eliminates the exponential componentes and returns the order of the model, the remainder signal and the parameters. At same time that we eliminate each component, we isolate it to determine its parameters (quatitative or parametric identification). QSI(X) degree=0 P=0; Parameter Matrix (X,k,P)=TPAFilter(X,degree,P) (X,k,P)=ExpFilter(X,degree,P) return Model(degree,X,P) Fig. 1. QSI Algorithm
The function T P AF ilter eliminates one sinusoidal component, at a time, starting with the highest frequency component, the remainder X ∗ (t) contains the summation of all the previous components except the eliminated one. After the elimination of j sinusoidal components, the remainder is: Xj∗ (t) = En1 (t) + ESn2 −j (t)
(6)
The elimination of components continues until the rest of the signal is non-oscillatory. The remainder signal, after extracting the oscillatory components, is a sum of exponential terms and is filtering by ExpF ilter function, this filter also performs a 5
successive elimination of components, starting with the component with the slowest decay rate. Figure 2 illustrates the application of QSI to a sample signal. 6
4
magnitude
2
0
−2
−4
−6
0
2
4
6 time
8
10
12
Fig. 2. Component Elimination Process.
4 Filters
Digital filters have been used mainly for two general purposes: signal separation, for signals that have been previously mixed, and signal restoration, for signals that have been distorted. In our application, we use filters to separate the different components that form the impulse response of a LTI system.
4.1 Two-Point Average Filter
A modification of a two-point difference filter [Smith1999] is proposed as an alternative to perform the first step in the filtering process. This first step of the filtering process detects the points where the derivative of the observed time series changes. The derivative can be computed using Equation 7 x (t) =
x(t + 1) − x(t − 1) 2∆t
(7)
where ∆t is the sampling period, x (t) is the derivative of x at point t, and x(t) is the input signal at time point t. Using the approximation given by Equation 7 we determine the time points where the sign of the derivative changes, i.e., we are detecting the extrema of x. Once computed the vectors of maxima and minima, we model the envelope by a cubic spline interpolation using the points at the maxima and minima, respectively, 6
to smoothen the upperhigher and lower envelopes. Once determined the envelopes, we compute the carrier by averaging the envelopes at each point of the original signal. This is a simple implementation of a low-pass filter. Using this filter, the midpoint vector m is computed. Vector m represents the remainder X ∗ (t) defined in Equation 6. A cubic spline interpolation approximation is used to smoothen the form of the carrier and has proven to provide better results than the bare filtered data. These ideas were yield to implement the filter that eliminated the sinusoidal componentsof the response. Figure 3 shows the TPA filtering algorithm. (TPAFilter).
TPAFilter(X, k,P) Xtrmax =determinem axima(X) Xtrmin =determinem inima(X) if oscillatory(X) k =k+2 //Compute Envelope Envelopemax =CubicInterp(Xtrmax ) Envelopemin =CubicInterp(Xtrmin ) //Compute m, middle point vector X ∗ =CubicInterp(m) S = X − X∗ P=P+SinP aram(S) (X, k,P)=TPAFilter(X ∗ , k,P) else X∗ = X return (X, k,P) Fig. 3. TPAFilter Algorithm
Figure 4 illustrates one step in the elimination process as performes by the algorithm in figure 3. Extrema are marked with plus signs and the computed middle points with circles. The signal formed by the middle points constitutes the carrier of the original wave m (i.e. the remainder). 7
0.5
Maxima envelope
Middle points
0
−0.5
−1
Remainder X*(t) −1.5
Minima envelope −2
1.8
2
2.2
2.4
2.6
2.8
3
3.2
Fig. 4. One step in the sinusoidal component elimination process.
4.2 Exponential Filter
Once we have filtered the oscillatory components of the signal, the remainder will be a summation of exponentials. In this step we apply a filter to extract each exponential component. In this case, to make the extraction of the exponentialcomponents, it becomes necessary to compute the parameters of the component that will be extracted. In the following section a complete description of the parametrizacin process is developed. Figure 5 presents the algorithm of this filter, and the Figure 6 shows its application. The function ExpFilter eliminates the exponential components at a time starting with the slowest decay rate component. Using the approximation given by equation 7 we determine the extrema of the TPAFilter’s remainder, i.e., the number of exponential components (see theorem 2). The exponential components are eliminated by means of the following process: first we detected the last inflection point (Pi ) and the points P1 and P2 are selected (see Figure 6). The points P1 and P2 are used to calculate the exponential component, that will be eliminate. The remainder X ∗ (t) contains the summation of all the exponentials components except of the eliminated one. In this case, after the elimination of i exponential component, the remainder is Xi∗ (t) = En1 −i (t)
(8)
This process repeats until all the components have been eliminated. In case that the extrema of remainder be zero, ExpFilter verifies if the residual has exponential form or not. For each eliminated exponential component the algorithm adds an unit to the order 8
ExpFilter(X,k,P) Xtr=count extrema of X if Xtr > 0 do k=k + 1 Pi =InflectionPoint(X) P=ExpParam(X,Pi ) X=X − Exponential(P ) Xtr=Xtr − 1 while Xtr > 0 else k=k + 1 P=ExpParam(X) X=X − Exponential(P ) return (X, k,P) Fig. 5. Exponential Filter Algorithm
of the model and at the same time that filter the components their parameters are calculated (see next section).
5 Quantitative Identification
To extend the capabilities of QSI we can determine the parameters of the response at the same time we are filtering each component. The next two subsections show how to perform this task for the exponential and sinusoidal cases.
5.1 Exponentials components
If the remainder signal is non-periodical or non-oscillatory it is a summation of exponentials. The basis to separate the remainder exponential is theorem 3: Theorem 3 If x1 , x2 , ..., xn are components of a sum of exponentials E(t) with decay rates r1 , r2 ...rn respectively, and r1 > r2 > ... > rn , then x1 (t∗ )
Juan J. Flores a , Nelio Pastor a,∗ a Divisi´ on
de Estudios de Postgrado, Facultad de Ingenier´ıa El´ectrica, Universidad Michoacana de San Nicol´as de Hidalgo, Morelia, M´exico
Abstract The problem of Systems Identification starts with a time-series of observed data and tries to determine the simplest model capable of exhibiting the observed behavior. This optimization problem searches the model from a space of possible models. In this paper, we present the theory and algorithms to perform Qualitative and Quantitative Systems Identification for Linear Time Invariant Dynamic Systems. The methods described here are based on successive elimination of the components of the system’s response. Sinusoidals of high frequencies are eliminated first, then their carrying waves. We continue with the process until we obtain a non-oscillatory carrier. At that point, we determine the order of the carrier. This procedure allows us to determine how many sinusoidal components, and how many exponential components are found in the impulse response of the system under study. The number of components determines the order of the system. The paper is composed of two important parts, the statement of some mathematical properties of the responses of Linear Time Invariant Dynamic Systems, and the proposal of a set of filters that allows us to implement the recognition algorithm. We present the aplication of the proposed methodology to analysis and modeling of electrical circuits and electrical power systems.
Key words: Dynamic systems, systems identification, differential equations, qualitative reasoning
This research has been supported by CIC-UMSNH grant No. 10.3 ∗ Correspondig author. Address: Maestro de Avila 90, Fracc. Fray Antonio de San Miguel, Morelia, Michoac´an, 58270, M´exico Email addresses: [email protected] (Juan J. Flores), [email protected] (Nelio Pastor).
Preprint submitted to Engeenering Applications of Artificial Intelligence 10 March 2004
1 Introduction
System identification deals with the problem of building mathematical models of dynamical systems based on observed data from the system and tries to determine the simplest model capable of exhibiting the observed behavior. This optimization problem searches the model from a space of possible models. In traditional methods, the process of structural identification has received less attention than the parametrization of the model. In most of cases the structural estimation is not generally made. Typically, a certain model structure is chosen by the user which contains unknown parameters. The choice of an appropiate model structure is most crucial for a successful identification application. This choice must be based both on an understanding of the identification procedure and on insights and knowledge about the system to be identified. An alternating way to infer a suitable structure, guided by system knowledge and the collected data set, is presented in this paper. In our approach we limit the search space to Linear Time-Invariant models.
The problem we are assessing in this paper is that of Qualitative Systems Identification for Linear Time Invariant (LTI) Dynamic Systems. This kind of systems can be represented by Linear ordinary Differential Equations with Constant Coefficients. Although this problem may seem over-constrained, there are many important problems in engineering and physics that can be expressed in mathematical terms by this kind of differential equations. We could even assess time-varying or even non-linear systems, if we consider them as piece-wise decomposed by linear approximations of the original systems. So, this is a limited yet interesting domain to start with.
The Systems Identification process can be decomposed in two steps: the first step called structural (or qualitative) identification, involves determining the qualitative features of the mechanism, i.e. the qualitative form of the system inside the black box [Bradley, Easley, & Stolle2001] [Bellazzi, Guglielmann, & Ironi1999] [Bradley & Stolle1996] [Kay, Rinner, & Kuipers2000]; once we know the nature of the mechanism, in the second step we proceed to determine the numerical value of the parameters of the model determined in the first step.
This second part can be done by any optimization process (e.g. minimum square error [Ljung1987] [P.1989], genetic algorithms [Goldberg1998] [Haupt & Haupt1998] [Hunt1993] [Pastor2000] [Zhang Zibo1987], etc.). Some algorithms use an optimization process to determine both faces at the same time (e.g. genetic algorithms [Downing et al.1996] [Kristinsson & Dumont1992]). 2
2 Linear Time-Invariant Dynamic Systems
Linear ODEs with constant coefficients are the most studied kind of differential equations; they have complete analytical solutions. So, if we restrict the kind of systems we are to analyze to those that can be expressed by an nth-order ordinary differential equation with constant coefficients, we know the kind of responses we are to get. We can express the behavior of a system in terms of the exponential and sinusoidal components in the response. We define En1 (t) =
1≤i≤n1
ai eri t
(1)
as a summation with at most n exponential terms, and ESn2 (t) =
1≤i≤n2
ai eri t sin ωi t
(2)
as a summation of exponentially decreasing sinusoidal functions. Note that we are not interested in giving analytical solutions to the differential equation, but a qualitative description of its behaviors. Theorem 1 The response of an nth order system can be expressed as in Equation 3. X(t) = En1 (t) + ESn2 (t)
(3)
where n1 + 2n2 = n. This result is evident from the definitions of Equations( 1) and (2). We see that if the second term of Equation 3 does not exist, the response will be acyclic. Otherwise, it is a sinusoidal wave, where En1 (t) represents its axis or attractor, and ESn2 (t) its sinusoidal components. Note that if we include a forcing function, the system’s response would be decomposed into Natural (the solution to the homogeneous equation) and Forced responses. If we restrict the forcing functions to be of the form e αt sin βt (i.e. constant, exponential, or sinusoidal), the forced response always has the same qualitative form as the forcing function [Boyce & DiPrima1969]. The force response would preserve the qualitative form of the response, and only adds one more exponential or sinusoidal term to the response. Let us analyze the qualitative form of the responses, as expressed by Equation (3). This qualitative form can be derived from the qualitative form of its exponential and sinusoidal components. 3
2.1 Exponential Components
The qualitative behavior of the exponential part of the response is characterized by Theorem 2.
Theorem 2 X(t) = En (t) = 1≤i≤n ai eri t exhibits at most n extrema (maxima or minima), including the one when t → ∞. We will assume, without loss of generality, that r 1 > r2 > . . . > rn . Theorem 2 is equivalent to saying that the derivative X (t) has at most n different zeroes. That is, X (t) = −r1 a1 e−r1 t − . . . − rn an e−rn t = 0 = r1 a1 e−r1 t + . . . + rn an e−rn t = 0
(4)
Performing the variable change z = e−t , Equation (4) becomes −X (t) = r1 a1 z r1 + . . . + rn an z rn = 0
(5)
Equation (5) is a polynomial in z of degree r 1 . Based on Descartes’ theorem [Kurosch1977,Rees, Sparks, & Re the number of roots of a polynomial are determined by the sign changes in the polynomial coefficients. Since we have n terms in polynomial (4), we can have at most n − 1 changes of sign, and therefore n − 1 roots. There is one more root placed where the function becomes zero, and that is when t → ∞. Given we have at most n roots for Equation (4), we can therefore have at most n extrema for En (t).
2.2 Sinusoidal Components
If the frequencies of two sinusoidal components of Equation (3) are equal, their shapes are reduced to one [Flores & Farley1995]. If their frequencies are different, they can be seen as the faster sinusoidal mounted on the slower one. If two frequencies of the sinusoidal components are very close together, the resulting wave beats. The results presented in this section fully characterizes all possible responses of a LTI Dynamic System. In the next section we will describe how to use these results to produce a framework for performing Systems Identification at the qualitative and quantitative levels. 4
3 QSI Algorithm
The identification algorithm presented in this section is based on the fact that the response of a LTI system can be decomposed in a sumation of exponential terms. If some of those exponential terms are complex, in which case they are conjugate complex pairs, each pair forms a sinusoidal. The algorithm that we propose in this paper takes as input a time-series representing the response of an LTI dynamic system and is capable of separating the terms of Equation (3) to determine the structure or qualitative form of the system exhibiting the observed behavior. Separating the terms of the system’s response is performed by a filtering process. Figure (1) shows a simplified version of the QSI algorithm. QSI determines the order of the system by adding the order of all eliminated components. This algorithm analyses the response and eliminates each component at a time, there are two main functions in QSI algorithm: T P AF ilter and ExpF ilter. The function T P AF ilter (see next section) eliminates the sinusoidal components and returns the order corresponding to those components, the remainder signal and the parameters. The function Expf ilter eliminates the exponential componentes and returns the order of the model, the remainder signal and the parameters. At same time that we eliminate each component, we isolate it to determine its parameters (quatitative or parametric identification). QSI(X) degree=0 P=0; Parameter Matrix (X,k,P)=TPAFilter(X,degree,P) (X,k,P)=ExpFilter(X,degree,P) return Model(degree,X,P) Fig. 1. QSI Algorithm
The function T P AF ilter eliminates one sinusoidal component, at a time, starting with the highest frequency component, the remainder X ∗ (t) contains the summation of all the previous components except the eliminated one. After the elimination of j sinusoidal components, the remainder is: Xj∗ (t) = En1 (t) + ESn2 −j (t)
(6)
The elimination of components continues until the rest of the signal is non-oscillatory. The remainder signal, after extracting the oscillatory components, is a sum of exponential terms and is filtering by ExpF ilter function, this filter also performs a 5
successive elimination of components, starting with the component with the slowest decay rate. Figure 2 illustrates the application of QSI to a sample signal. 6
4
magnitude
2
0
−2
−4
−6
0
2
4
6 time
8
10
12
Fig. 2. Component Elimination Process.
4 Filters
Digital filters have been used mainly for two general purposes: signal separation, for signals that have been previously mixed, and signal restoration, for signals that have been distorted. In our application, we use filters to separate the different components that form the impulse response of a LTI system.
4.1 Two-Point Average Filter
A modification of a two-point difference filter [Smith1999] is proposed as an alternative to perform the first step in the filtering process. This first step of the filtering process detects the points where the derivative of the observed time series changes. The derivative can be computed using Equation 7 x (t) =
x(t + 1) − x(t − 1) 2∆t
(7)
where ∆t is the sampling period, x (t) is the derivative of x at point t, and x(t) is the input signal at time point t. Using the approximation given by Equation 7 we determine the time points where the sign of the derivative changes, i.e., we are detecting the extrema of x. Once computed the vectors of maxima and minima, we model the envelope by a cubic spline interpolation using the points at the maxima and minima, respectively, 6
to smoothen the upperhigher and lower envelopes. Once determined the envelopes, we compute the carrier by averaging the envelopes at each point of the original signal. This is a simple implementation of a low-pass filter. Using this filter, the midpoint vector m is computed. Vector m represents the remainder X ∗ (t) defined in Equation 6. A cubic spline interpolation approximation is used to smoothen the form of the carrier and has proven to provide better results than the bare filtered data. These ideas were yield to implement the filter that eliminated the sinusoidal componentsof the response. Figure 3 shows the TPA filtering algorithm. (TPAFilter).
TPAFilter(X, k,P) Xtrmax =determinem axima(X) Xtrmin =determinem inima(X) if oscillatory(X) k =k+2 //Compute Envelope Envelopemax =CubicInterp(Xtrmax ) Envelopemin =CubicInterp(Xtrmin ) //Compute m, middle point vector X ∗ =CubicInterp(m) S = X − X∗ P=P+SinP aram(S) (X, k,P)=TPAFilter(X ∗ , k,P) else X∗ = X return (X, k,P) Fig. 3. TPAFilter Algorithm
Figure 4 illustrates one step in the elimination process as performes by the algorithm in figure 3. Extrema are marked with plus signs and the computed middle points with circles. The signal formed by the middle points constitutes the carrier of the original wave m (i.e. the remainder). 7
0.5
Maxima envelope
Middle points
0
−0.5
−1
Remainder X*(t) −1.5
Minima envelope −2
1.8
2
2.2
2.4
2.6
2.8
3
3.2
Fig. 4. One step in the sinusoidal component elimination process.
4.2 Exponential Filter
Once we have filtered the oscillatory components of the signal, the remainder will be a summation of exponentials. In this step we apply a filter to extract each exponential component. In this case, to make the extraction of the exponentialcomponents, it becomes necessary to compute the parameters of the component that will be extracted. In the following section a complete description of the parametrizacin process is developed. Figure 5 presents the algorithm of this filter, and the Figure 6 shows its application. The function ExpFilter eliminates the exponential components at a time starting with the slowest decay rate component. Using the approximation given by equation 7 we determine the extrema of the TPAFilter’s remainder, i.e., the number of exponential components (see theorem 2). The exponential components are eliminated by means of the following process: first we detected the last inflection point (Pi ) and the points P1 and P2 are selected (see Figure 6). The points P1 and P2 are used to calculate the exponential component, that will be eliminate. The remainder X ∗ (t) contains the summation of all the exponentials components except of the eliminated one. In this case, after the elimination of i exponential component, the remainder is Xi∗ (t) = En1 −i (t)
(8)
This process repeats until all the components have been eliminated. In case that the extrema of remainder be zero, ExpFilter verifies if the residual has exponential form or not. For each eliminated exponential component the algorithm adds an unit to the order 8
ExpFilter(X,k,P) Xtr=count extrema of X if Xtr > 0 do k=k + 1 Pi =InflectionPoint(X) P=ExpParam(X,Pi ) X=X − Exponential(P ) Xtr=Xtr − 1 while Xtr > 0 else k=k + 1 P=ExpParam(X) X=X − Exponential(P ) return (X, k,P) Fig. 5. Exponential Filter Algorithm
of the model and at the same time that filter the components their parameters are calculated (see next section).
5 Quantitative Identification
To extend the capabilities of QSI we can determine the parameters of the response at the same time we are filtering each component. The next two subsections show how to perform this task for the exponential and sinusoidal cases.
5.1 Exponentials components
If the remainder signal is non-periodical or non-oscillatory it is a summation of exponentials. The basis to separate the remainder exponential is theorem 3: Theorem 3 If x1 , x2 , ..., xn are components of a sum of exponentials E(t) with decay rates r1 , r2 ...rn respectively, and r1 > r2 > ... > rn , then x1 (t∗ )
