A parametric frequency response method for non-linear time-varying systems

A new parametric frequency response algorithm is introduced to investigate linear and non-linear dynamic systems with time-varying parameters. In the new algorithm the time-varying parameters are regarded as additional inputs of the systems and the non-linear generalised frequency response functions for multi-input-single-output systems are then employed to obtain Zadeh's system functions from a differential equation representation. The parametric frequency response method reveals how the time-varying parameters affect the behaviour of the systems through a time-varying term. The new method can be applied to both linear and non-linear time-varying systems.


Introduction
Time-varying systems occur frequently in engineering systems, including control systems and communications. Sometimes the time variation of parameters can be ignored, such as in ageing and the deterioration of electronic components. These systems can be considered to be essentially time invariant. When the parameters change quickly and the speed of change is comparable to the time constants of the system the time variation of the systems cannot be ignored and new methods for time-varying systems are needed. In the past, different representations and techniques have been introduced for the investigation of linear time-varying systems (Zadeh 1950a(Zadeh , 1950bPuri and Weygandt 1963;Wierwille 1963;Kozek and Hlawatsch 1991;Shenoy and Parks 1994;Hlawatsch 1998) and to a much lesser degree non-linear time-varying systems (Bansal 1969;Sandberg 1982;Sandberg 1983;Yuan and Opal 2001). The stability of time-varying systems has also been studied by several authors (Bongiorno 1964;Sandberg 1964;Zames 1966;Brockett and Lee 1967;Zhu and Johnson 1988). These representations and methods for time-varying dynamic systems mainly fall into three classes: differential equation or Volterra series based time domain methods, Zadeh's system functions-based frequency domain methods and the Wigner distribution and ambiguity function-based time-frequency domain methods.
The system functions introduced by Zadeh (1950b) are a natural extension of the traditional frequency response function to linear time-varying systems. Zadeh's system function is a two-dimensional function of both frequency and time which defines a time-varying frequency response function. The system function represents a frequency * Corresponding author. Email: s.billings@shef.ac.uk relationship between the input and output at an instant in time. Zadeh showed that a system function satisfies a differential equation and the system function can be expanded as a series of fixed (time-invariant) frequency response functions by approximating the solution using perturbation methods. Taking the Fourier transform with respect to the time parameter of the system functions then yields bi-frequency system functions. The bi-frequency system functions have subsequently been developed by Bayan and Erfani (2005), and have been extended to non-linear time-varying cases by some authors (Bansal 1969;Yuan and Opal 2001). Correlation functions and power spectra in variable networks have also been studied (Zadeh 1950a).
Although Zadeh's system functions inherit many merits of frequency response functions and have a clear physical meaning, the application of this representation is limited. The primary drawback to this method is that it is limited to slowly varying systems or small variations so that perturbation methods can be used (Gibson 1963). These difficulties become much worse when this idea is applied to non-linear time-varying systems. In this paper a new parametric frequency response (PFR) method is introduced to study time-varying systems based on results from GFRF (generalised frequency response function) methods. A new algorithm is introduced where the time-varying parameters are regarded as additional inputs to the systems and results for MISO (multi-input-single-output) time-invariant systems are used to obtain the system function representations for time-varying systems. A related idea can be found in Lawrence's paper (1979) which was used to solve timevarying linear differential equations. In this paper the relationships between the system function representation, the C 2013 The Author(s). Published by Taylor & Francis. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Volterra series and the associated GFRFs will be studied. Generalised frequency response functions represent extensions of the classical linear frequency response function to non-linear systems. It is shown that the probing method for MISO non-linear systems (Worden, Manson, and Tomlinson 1997;Swain and Billings 2001) can be used to derive the system functions. One advantage of the new parametric frequency response method is that it can provide deep insight into time-varying systems. The new parametric method separates the effects of inputs and time-varying parameters to reveal how time-varying parameters imposed at different positions in a differential equation representation affect the behaviour of the system through the nonlinear terms. In addition, the new parametric frequency response method applies to non-linear time-varying systems because all the techniques used are based on non-linear representations.
The new parametric frequency response method could be applied for example in the design of aircraft control systems based on the gain-scheduling method (Leith and Leithead 2000). Using this method a family of linear timeinvariant controllers are designed at different operating points and the parameters in the controllers are often fitted as a function of the scheduling variables, such as altitude and Mach number. Information on how a control parameter affects the performance of the systems can be very useful and is often crucial for system design. The new parametric frequency response method provides a clear path to an explanation of this kind of information.
This paper is organised as follows: Section 2 briefly reviews Zadeh's system functions for linear time-varying dynamic systems and the extension to non-linear timevarying systems. The new parametric frequency response method is introduced in Section 3 to exploit the useful information in the system function representations.
Two illustrative examples are given in Section 4 to illustrate the new method. Conclusions are finally drawn in Section 5.

System functions and bi-frequency system
functions Consider a continuous time LTV (linear time-varying) system denoted as where L (t, •) represents the time-varying system and u (t) and y (t) are the input and output of the system, respectively. The relationship between the input and output of the system can be described by a convolution of the input and the weighting function For a causal system with a zero initial condition, Equation (2) becomes because h(t, τ ) = 0, when t < τ. In the remainder of this paper, all the systems under study are assumed to be causal. The weighting function h (t, τ ) is also known as the unit impulse response function which defines the response of a linear time-varying system at time t to an impulse input imposed at τ . This provides a two-dimensional representation of the system involving the observed time t and the application time τ . For an LTI (linear time invariant) system L (•) the two time-dimensions are related and the behaviour of the system depends only on the difference t − τ . The Fourier transform of the time-invariant impulse response function with respect to the value t − τ , denoted as H (jω), is the well-known linear frequency response function of the system. Signal e jωt is an eigenfunction of the LTI system and H (jω) is the associated eigenvalue, that is When the system is linear time-varying although the superposition applies the system exhibits behaviour quite unlike LTI systems. A new representation, a system function, was introduced by Zadeh (Zadeh 1950b) for LTV systems The system function representation transfers the application time τ of the impulse response function into the frequency domain providing a time-frequency description of the system. Zadeh's system function is a natural extension of the frequency response function for LTI systems.
The response of the LTV system with respect to input e jωt is Therefore, the system function can be alternatively defined as: response of an LTV system to e jωt e jωt .
Substituting the input in (2) with an inverse Fourier transform, the output of an LTV system can then be expressed using the new defined system function as: Now, the output of the LTV system becomes an inverse Fourier transform of the product of the system function and the Fourier transform of the input W (t, j ω) U (jω).
The frequency domain output of the system can be given as follows by applying the Fourier transform to both sides of (8) Equation (10) gives the definition of the bi-frequency system function. The bi-frequency system function is essentially the Fourier transform of W (t, j ω) e jωt .
System functions have been extended to non-linear time-varying systems by several authors (Bansal 1969;Yuan and Opal 2001). Consider the time-varying non-linear system which can be expanded in a convergent Volterra series where w k (t, τ 1 , . . . , τ k ) represent the time-varying Volterra kernels. The study of Volterra series expansions of timevarying non-linear systems can be found in the papers by Sandberg (1982Sandberg ( , 1983).
The generalised system functions (or high-order system functions) can then be defined as: Similarly the corresponding bi-frequency system functions can be defined as: Notice that although there are more than two frequencies involved, the name 'bi-frequency system function' is used for convenience of expression.

The parametric frequency response method
In this section a new general form of the high-order system functions for non-linear time-varying systems will be derived using the new parametric frequency response method. System functions for LTV systems are particular cases of general results when k = 1.
Consider a time-varying system expressed as the following model: represents the rth order non-linearities composed of the input and output and the associated derivatives; α r,p (t, c 1 , c 2 , . . . , c r ) represents timevarying coefficients. The term reduces to a pure input non-linearity when p = 0 and a pure output non-linearity when p = r. System (14) reduces to a linear time-varying system when R = 1 and the associated multi-input time-varying system becomes a bi-linear form.
The non-linear model is composed of a linear combination of non-linear terms and is linear-in-the parameters. In practical applications, the time-varying model is often derived from the first principles or obtained using a system identification method such as the wavelet-NARMAX algorithm (Billings and Wei 2005). In this model, the time-varying parameters α r,p (t, c 1 , c 2 , . . . , c r )'s can be considered as additional inputs so that the system becomes a time-invariant system with more than one input.
For example, a specific instance of the non-linear ODE (ordinary deferential equation) model in (14) y +ẏ + ρ cos (μt) yẏ + y = bu (15) can be produced from the general form by setting the coefficients as Assume the time-varying system admits a Volterra series. Volterra series representation has been widely studied and many results on the properties, including convergence analysis are available (Ku and Chin Chun 1967;Brockett 1976;Gilbert 1977;Lesiak and Krener 1978;Sandberg 1982;Sandberg 1983). GFRFs are defined as the Fourier transform of the Volterra kernels and GFRFs of MISO time-invariant systems have been studied by Worden et al. (1997) and Swain and Billings (2001). GFRFs can easily be obtained from non-linear differential equation models using a harmonic probing method. Recent studies have shown that the Volterra series representation can also be obtained using an Adomian method when the timevarying system is analytic with respect to the output and linear with respect to the input (Guo, Guo, Billings, Coca, and Lang 2012).
The new parametric frequency response method exploits the connection between the direct representation for time-varying systems and representations for associated multi-input time-invariant systems. This allows results for MISO time-invariant systems to be used to investigate the behaviour of time-varying systems.

Derivation of system functions from the Volterra series representation
Consider the time-varying system (14) and assume the system has a zero equilibrium point, that is, the response of the system is zero when the input u(t) is set as zero. Notice that this means that zero is also an equilibrium point of the MISO system when the time-varying parameters are treated as additional inputs. It is also assumed that both input and the time-varying parameters have finite energy so that the Fourier transform can be applied.
Assume the time-varying system admits a Volterra series representation and the time-varying Volterra series is expressed as where the kernels w n (t, τ 1 , . . . , τ n ) are time-dependent. Assume that there is a total of d time-varying parameters in system (14) denoted as α 1 (t), . . ., α d (t) and one input u (t). The system is also assumed to have a zero equilibrium point, which means the response of the system is zero when the input is set as zero. When the time-varying parameters are considered as additional inputs the system becomes a time-invariant system with many inputs. The time-invariant Volterra series representation can then be given as a sum of direct-kernels of each input and the cross-kernels between different inputs.
where k 1 ∼ k m can be any integer between 1 and d.
Rearrange the Volterra series in (18) according to the order of the input.
The first term on the right-hand side of (19) is which is composed of direct and cross kernels of the total of d-varying parameters. This term will be zero because zero is an equilibrium point of the system. The second term, which is composed of direct kernels with respect to the input u (t), represents the contribution purely from the input. The last terms represent the combined contributions from both the input and the time-varying parameters.
Therefore, Equation (19) can further be written as follows by combining the second and the third terms.
Comparing Equations (21) and (17), the relationship between representations (17) and (18) can be expressed as Generalised system functions can then be defined as and the output of the system can be expressed using the system functions as For the linear time-varying system the only system function will be W 1 (t; jω 1 ), usually written as W (t, j ω), which agrees with Zadeh's definition.
Substituting (22) into the definition of the system function in (23) yields W k (t; jω 1 , . . . , jω k ) Equation (25) provides an explicit expansion for the new generalised systems functions. The computation of the expansion is much simpler compared with perturbation methods when the Volterra kernels of the multi-input time-invariant system are known. Furthermore, the new expansion provides more detail about the position of the timedependent parameters within the model.
Equation (25) shows that the response of the system consists of two parts: the time-invariant response and the time-varying response. The first term on the right-hand side of (25) presents the time-invariant part of the system when m = 0 and represents the output of the system purely under the effect of the input when the time-varying parameters are zero. Notice that the first term in (25) is H k (jω 1 , . . . , jω k ) when m = 0. This means that the generalised system functions reduce to frequency response functions when all the time-varying parameters are equal to zero. This agrees with the generalised frequency response function representation of a time-invariant system. When the parameters change periodically the expansion becomes much simpler and easy to understand In this case it is easy to observe that the system function can be expanded as a summation of a series of timeinvariant frequency response functions. This is coincident to the conclusions in Zadeh's paper.
However, the calculation of the Volterra kernels directly from a differential equation model could be very difficult, whereas the computation of the associated generalised frequency response functions is much simpler using results based on the harmonic probing method (Worden et al. 1997;Swain and Billings 2001;Worden and Tomlinson 2001;Peyton Jones 2007). A new expansion of the generalised system functions for non-linear time-varying systems with generalised frequency response functions will therefore be derived next.
Express the time-varying parameters in a Fourier transform form: and substitute (27) into (25). The generalised system function becomes Equation (28) provides a different form of expansion of the generalised system functions where now the generalised frequency response functions are used instead of Volterra kernels. The new expansion of system functions can now be interpreted. The new expansion shows how the timedependent information from the time-varying parameters is packed into the kth order system function. The generalised system functions can also be regarded as a kind of timedependent generalised frequency response function. Additionally, Equation (28) is a generalisation of Equation (26) for a more general form of time-varying parameters.
According to the definition (13), the kth order bifrequency system functions can be defined as and the output of system can be given in the frequency domain as Although the definition of the bi-frequency system function was given by Zadeh he did not give much explanation of the new introduced frequency υ. The (generalised) bi-frequency system functions are recognised as important for the frequency domain study of time-varying systems especially for the study of the output properties of systems. The next section will reveal the importance underlying the (generalised) bi-frequency system functions using the new parametric frequency response method.

Derivation of bi-frequency system functions from
the generalised frequency response functions In this subsection the kth bi-frequency system function will be derived from the generalised frequency response function representations.
Following the discussion as the end of Section 2, the bi-frequency system function k (jυ, jω 1 , . . . , jω k ) and W k (t; jω 1 , . . . , jω k ) e j k i=1 ω i t are a pair of Fourier transforms, that is Substituting (28) into (24), the output of the system can be calculated by Define a new frequency denoted as υ The frequency represents a combination of the frequencies from both input and time-varying parameters.
Replace the variable μ m = υ − k i=1 ω i − m−1 i=1 μ i and the left-hand side of (31) can be written as Compare Equation (34) with (31). The kth order bifrequency system function can be written as k (jυ, jω 1 , . . . , jω k ) According to the results in Lang and Billings (1996), the integral in (35) can be written as a surface integral on the hyper-plane D : μ = μ 1 + μ 2 + · · · + μ m−1 + k i=1 ω i , that is, k (jυ, jω 1 , . . . , jω n ) The new parametric frequency response method provides additional insight into bi-frequency system functions. The frequency υ represents the output frequency that is a combination of all the frequencies from both the input and the time-varying parameters. Equations (35) and (36) provide an explicit expansion of the kth order bi-frequency system function in the generalised frequency response functions, which give an interpretation about how the frequency components from the input and the time-varying parameters make a contribution to the final output frequency through the time-varying terms.
As a result, the output of the system can be written as This is coincident with the results in (30). When the input and time-varying parameters have a discrete spectrum the integrals in (35) and (36) become a summation of these discrete frequencies which satisfy (33) and the results become simpler. These results are not given here but can be seen in the first example in Section 4. Finally, the system functions and the bi-frequency system function can be approximated with a truncated series when the generalised frequency response function representation of the MISO system converges. Results on the convergence of the GFRF representations can be found in the references (Tomlinson, Manson, and Lee 1996;Peng and Lang 2007).

Illustrative examples
In this section two illustrative examples are studied to demonstrate the new parametric response method for a time-varying system with a discrete spectrum, that is, the time-varying parameters take a simple form as the summation of a finite number of frequencies.

A linear time-varying system
A first-order linear time-varying system which was used in Zadeh's paper will be considered.
Rewriting the system in a more general form by replacing the time-varying parameters with a function of time yields where e 1 (t), e 2 (t) represent the input and output of the system, respectively and α (t) represents the time-dependent parameters.
Consider the time-varying parameter as an additional input so that the system becomes a time-invariant MISO system. Using the harmonic probing method for MISO systems (Worden et al. 1997;Worden and Tomlinson 2001), the first few direct and cross generalised frequency response functions can be calculated as follows: H e 1 e 1 2 (jω 1 , jω 2 ) = 0.
Observe that the higher order direct generalised frequency response functions about the input H e 1 ···e 1 k (jω 1 , . . . , jω k ), k ≥ 2 are zero. This is because the system is a linear system and the non-linearity arises from the multiplication of the time-varying parameters with the input and output. All the direct generalised frequency response functions with respect to the time-varying parameter α (t) are zero because e 2 (t) = 0 when e 1 (t) = 0 no matter how the parameter α (t) changes.
The time-varying parameter can be decomposed as a pair of frequency components.
According to the results in (26) the system function can be expanded as where μ k i = μ 0 or −μ 0 .
Substituting the first few GFRFs into (47) yields Observe that the result obtained using the parametric frequency response method is exactly the same as the result that Zadeh obtained (Zadeh 1950b) where two different perturbation methods were used. However, the parametric frequency response method is more straightforward even when a higher order approximation is needed. The new parametric frequency response method also provides a serial expansion of bi-frequency system functions which reveal an important insight into the output properties of the system.
The bi-frequency system function can be calculated using the expansion given in the last section (jυ, jω) The results shows that, the output can only take a discrete spectrum with the output frequency υ = ω + pμ 0 where p is any integer when the input frequency ω is fixed. The magnitude and phase of the output of a certain frequency is determined by the expansion of the bi-frequency system function.

A non-linear example: the Duffing equation
The Duffing equation is often used as a typical prototype of non-linear dynamic systems and has been widely used as an example for the study of the Volterra approximation of non-linear systems. The Duffing equation admits a Volterra series representation around the zero equilibrium point when the input is within the radius of convergence.
Since only a third-order non-linearity is included in the Duffing equation, all the Volterra kernels of an even-order are zeros. A Duffing oscillator with a time-varying non-linear restoring force will be considered where α (t) represents the time-varying coefficient of the non-linear restoring force. The time-varying Volterra series expansion of the Duffing oscillator is given as h k (t, τ 1 , . . . , τ k ) × u (τ 1 ) · · · u (τ k ) dτ 1 . . . dτ k . Equation (69) shows that the steady-state response of the system to a harmonic input can be approximated using a truncated series of generalised system functions. Based on this idea the steady-state output of the Duffing oscillator is synthesised using the obtained generalised system functions. The comparison of the simulated output and the synthesised output is shown in Figure 1 where the parameters are set as m = 1, c = 1.4, k 1 = 1, k 3 = 4, A = 2.4, ω 0 = 2.5 rad/s, ρ = 2, μ 0 = 1.7 rad/s. Figure 1 shows the synthesised output approaching the simulated output when N increases.
From the top to bottom Figure 1 shows how the response of the system distorted gradually and approached the final output. The topmost subplot shows the linear part of the system response which is a stationary sinusoidal signal with a phase delay. Under the effects of the time-varying parameters the response of the system becomes non-stationary. Notice that in this example the timevarying parameter α (t) = ρ cos (μ 0 t) is neither small nor slowly varying. However the parametric frequency response method works very well.

Conclusions
A new methodology has been introduced to investigate both linear and non-linear time-varying dynamic systems. The new parametric frequency response method which treats time-varying parameters as additional inputs and uses welldeveloped techniques for multi-input-single-out non-linear systems is easy to operate. The new method provides an explicit expansion of the system functions. The new form expansion encapsulates the time-varying information scattered in the traditional generalised frequency response function representation and extends the generalised frequency response function representation to a wide class of timevarying dynamic systems. The new method gives a remarkable insight into how time-varying parameters affect the output of a system through time-varying terms.
The new algorithm has been compared with the perturbation methods used by Zadeh. Results show that the new method yields the same results as the perturbation methods but is more efficient avoiding the solution of differential equations. In addition, the time-varying parameters need not to be periodic in the new method. Another obvious advantage is that the new method can be applied to non-linear time-varying systems.
The new parametric frequency response method also reveals information underlying in the bi-frequency system functions. The expansion of the bi-frequency system functions in this paper is a powerful tool for the analysis of the output frequency of time-varying systems. The new parametric frequency response method applies to a wide class of systems and can be an important and generic tool for the analysis and design of time-varying systems.