Dynamic Response Simulation for a Non-linear System

Laboratory simulation testing has for many years contributed signiﬁcantly to the durability and quality of motor vehicles. Most sophisticated test rigs use an iterative algorithm that generates the input drive ﬁles that reproduce service environments under laboratory conditions. Essentially the algorithm solves a non-linear, multiple channel dynamic system. In this paper, the non-linear problem is recast as a system of algebraic equations. This mathematical framework allows the application of alternative but well understood solution techniques. Using mathematical simulations, conclusions are drawn concerning the choice of iteration gain in the current algorithm and the better performance of alternative numerical solution procedures.


Introduction
In laboratory simulation testing, a structure is mounted in a test rig and is excited in such a way that the service environment, as represented by a set of responses from transducers, is reproduced.It is believed that, when these responses are replicated, the complex stress field within the structure that occurs in service is also reproduced.The test rig and the test structure form a non-linear dynamic system and the problem to be solved is to determine the input to this unknown multiple channel non-linear system.The technology that achieves this was developed in the 1970's -see e.g.[1] -following the introduction of the hydraulic servo-valve, the construction of algorithms for quickly processing random data in terms of Fourier Transforms, and of course the development of more powerful computers.The technology is well summarised by Dodds and Plummer [2].Generally, the procedure is that the system is treated as linear and measured using spectral analysis.An inverse system is then defined before an iterative algorithm determines the required drive files.Work to improve the performance of the iteration algorithm has been carried on over the years by, among others, Raath [3] who has developed a time domain version of the algorithm as an alternative to the usual frequency domain implementation, and also by de Cuyper et al [4] who examine improvements in the identification of the non-linear system.
The work presented here reports on the realisation that the problem may be recast mathematically as a system of non-linear algebraic equations.The conventional iteration algorithm is in fact an example of more general computational techniques for solving such systems.In the paper, these more general methods are introduced, and an application of them is then demonstrated in simulations using a single degree of freedom non-linear mathematical model for the system, the Duffing equation.The new viewpoint involves both time and frequency domain considerations.Note that, for this paper, the single degree of freedom system employed differs from the multiple channel physical laboratory simulation test system.Cost of equipment and control of parameters were considerations, but also using a single channel meant that the work could concentrate on the non-linearity rather than interaction between channels.The latter will be studied at a later date.
Before introducing the new approach, the current algorithm is applied to the chosen simulation model, demonstrating the method and its characteristics in the face of various degrees of severity of non-linearity.The situation is then studied mathematically and it is shown how discretisation leads to a system of non-linear equations.After presenting some general methods for solving systems of non-linear equations, the current algorithm is then shown to be a particular case.Finally, the feasibility of the more general approach is explored by comparing the success of the results of alternative solution methods.

Current Algorithm
The current algorithm for achieving drive signals exists in several commercial software programs.For a description, the reader is referred to Dodds et al [2].
The procedure may be summarised as follows: • Measurements of the response of the system are made during normal operation or specified operating conditions.These measurements are edited to provide a target response.In this paper, the target response is generated by exciting the system with band-limited white noise.
• The frequency response of the test rig and specimen is measured using spectral analysis.
• The validity of the frequency response measurement and the test rig design is then established using multiple and partial coherence functions e.g.Bendat et al [5].
• An inverse frequency response function is computed and, from this, an initial drive file is derived using the target response.
• The drive file excites the system and produces a response, which is compared with the target response.The difference is then used to create a new drive file and the process continues as an iteration until an acceptable level of error is achieved.
The excitation data used for measuring the system consist of bandlimited white noise, represented by the components of a vector x := (x 0 , x 1 , . . ., x N −1 ).
The system response is sampled, yielding another vector y := (y 0 , y 1 , . . ., y N −1 ), where, for signals of period T, y i := y(t i ) with t i := iT /N for i = 0, 1, . . ., N − 1.In the system measurement, spectral analysis uses the Discrete Fourier Transform of these signals, for which the k th components are denoted by X k and Y k , respectively, for k = 0, 1, . . ., N − 1, represented by the transform pairs: The frequency response is based on the Cross Spectral Density estimate of the input and output signals as given by Bendat et al [5], p138, where T is the period of duration of the signals and ω k is the k th discrete frequency, and • • • denotes an expectation value.The auto power spectral estimates S xx (ω k ), S yy (ω k ) are defined in a similar manner and the frequency response function is then given by In the simulations to be presented here, a target response signal y D is determined by exciting the system using a sequence x D of random numbers generated as bandlimited white noise.The iteration process is described more mathematically in Figure 1.The fraction of the drive signal increment p (n)   which is fed back is stipulated by the iteration gain λ n , a positive scalar quantity not greater than unity, which is chosen manually.In practice, the full drive signal is not normally used in determining the first drive file since the approximations in the estimate of the model may lead to the system being damaged.Similarly, the gain during the iteration is generally less than one to ensure convergence of the iteration and is again chosen manually.

Example of Current Iteration
The behaviour of the current algorithm is illustrated using a model of the Duffing equation constructed in MATLAB/Simulink.The system being simulated represents a mechanical single degree of freedom, damped spring-mass system Measure H ? Initial drive x (1) ↔ X (1)   n := 0 subject to the initial conditions y(0) = ẏ(0) = 0.The mass is taken to be 100kg, the damping ratio ζ = 0.1 and the natural frequency is normalised to unity.The right-hand side is chosen so that the input and output signals have similar magnitudes.
When identifying the physical system, the normal practice is to use a large number of averages to improve the expectation value of equation( 3) and achieve a smooth frequency response function.Here, a small number of averages are taken, but the function is smoothed using a least squares fit.Numerical experiments suggest that the least squares fitting is as good as employing a large number of averages.The iteration algorithm is applied with λ n = 0.5, n = 0, 1, . .., to produce a sequence of response vectors y (n) , n = 0, 1, . . .which converge to y D .The results are summarised in Figure 9.The algorithm stops if the fractional Euclidean norm of the error vector i.e.
falls below 5%.The error in the response achieved is shown in Figure 4 as a function of time.In practice, when the iterations fail to converge, the operator is free to adjust the iteration gain.For example, at the higher non-linearity of k = 0.25, the gain would be reduced and the iteration restarted, at the expense of slowing the convergence.

System of Algebraic Equations
In this section the problem is restated in terms of a system of algebraic equations.This opens up the possibility of applying well known numerical techniques for solving such systems.In addition, it is shown how the conventional approach appears as a particular case.One such computational strategy is applied to the simulation introduced in the previous section.
The point of view proposed in this paper is to note that the sampled response vector y is a function of the input signal x as symbolised in Figure 5: or, in component form: . . .To illustrate this, the model problem of the previous section is considered.
Equation( 4) is discretised to produce a system of equations, thus providing explicit information about the vector valued function f and the corresponding Jacobian.
First of all, consider the linear system obtained by setting k = 0 in equation( 4).
The response is related to the input by a convolution in the time domain: where the discretised impulse response h := (h 0 , h 1 , . . ., h N −1 ), is the inverse Discrete Fourier Transform of the frequency response function, This convolution may be written as a matrix product: in which the The i th component of the vector equation( 11) yields the approximate value of the response y(t) at t = t i .Equation( 11) may be rewritten: which may be regarded as a discretisation of equation( 4) with k = 0.
This process is extended to approximate the whole of the left hand side of equation(4) at t = t i for non-zero k : for some appropriate circulant matrix C, such that the n th component of the DFT of the vector [Cy] is given by The discretisation of equation( 4), after division by k, may now be expressed as a vector equation: where [g(y)] i := k y 3 i , thus yielding an explicit form for the function inverse of f in equation (7).
The mathematical problem may be stated as follows: given a vector y D = (y D 0 , y D 1 , . . ., y D N −1 ) and a particular function f , determine a vector x, such that This is a system of non-linear algebraic equations for which the solution is readily seen to be f −1 (y D ).In practice, the explicit form of f is not known, but for a given vector x, the value of y = f (x) may be obtained by "running the system".

Iterative Solutions
This type of problem is common, and there are well-known computational techniques for solving equation( 16).For a survey of practical algorithms which may be used to solve systems of non-linear algebraic equations the reader is referred to a review by Martinez [6].All the methods considered are iterative.Starting from some initial approximation x (0) , a sequence of iterates, x (0) x (1) , x (2) . .., is generated which converge, ideally, to the desired solution.
In order to understand these techniques, a brief account of Newton's method for systems of non-linear equations is given.This algorithm follows from the Taylor expansion in several variables of f (x) about the current approximation where J f (x) denotes the Jacobian matrix of order N × N for the vector-valued function f (x) in equation( 7) the partial derivatives being evaluated at x.In the context of matrix algebra, vectors are considered as column matrices.
The Jacobian for the model problem may be constructed from equation(15) where The dependence of the Jacobian on y and, therefore, on x is clear.We note that where x and y are related by equation(15).
For the linear system (9), it may be seen, from its definition, that the Jacobian is given by i.e. a constant matrix.
For the general non-linear system, if x D is a solution to equation( 16), then, setting y = y D in (17), where Ignoring the error term in equation( 23) leads to the following iteration scheme: provided the Jacobian is non-singular at x (n) .This is Newton's method which is locally quadratically convergent -see e.g.[6].It may be rewritten as where However, since f is not known explicitly, the Jacobian matrix cannot be constructed.Hence, "quasi-Newton" methods are considered which generalise equation(24) to in which the matrix B n plays the role of J f (x (n) ).This may be recast as equation( 25), where The idea is that, starting from some initial estimate of the Jacobian, B 0 , this is then updated using a simple formula.A very common approach is based on a version of the secant method and was suggested by Broyden [7], in which the updated inverse matrix [B n+1 ] −1 may be expressed in terms of [B n ] −1 , thus enabling a computationally efficient implementation of ( 27), provided we can where δx (n) := x (n+1) − x (n) , and δey (n) := ey (n+1) − ey (n) .

Alternative Iteration Schemes
There are many variations of the basic iteration (27) -for other, similar, approaches see Martinez [6].One possibility, which is relevant to our interests, is to keep B n constant at some value B 0 .The conventional approach discussed earlier -which treats the system as if it were linear -fits into this scheme in which the matrix C h is based on the vector h, the impulse response as in equation( 12).This impulse response corresponds to the measured frequency response function.
Another possibility consists of using C h as an initial approximation to the Jacobian in the non-linear system.
and then using equation( 29) to produce the updates for (28).It may be noted here, that, for example, in the computation of x in equation( 11), the Fast Fourier Transform may be employed i.e. there is NO matrix multiplication performed.Indeed, all matrix multiplications are avoided by implementing the algorithm described in [8].This algorithm is a memory efficient approach which only requires scalar products of vectors to be computed.
However, these techniques are only locally convergent.That is, the initial approximations to the solution and the Jacobian must be good enough for convergence to follow.Even Newton's method may fail to converge for cases where there is a unique solution.

Improving Global Convergence
It was noted earlier, that, in the current algorithm, the iteration gain is reduced if the iterations diverge.In fact, a search can be conducted to determine a suitable iteration gain.In an attempt to achieve global convergence the basic iteration (25) is modified to allow a variable step in the search direction : where λ n , for n = 0, 1, . . .are real numbers lying between 0 and 1.The general iteration process is shown in Figure 6.The values of λ n may be constant, or, depending on the situation, the operator may vary them e.g.some circumstances may warrant a moderate step reduction (λ ≈ 0.5), while others may require larger reductions (λ 0.5).The value of λ yielding the minimum error may be estimated using a backward line search.To do this a merit function is defined as follows: or, to avoid a square root, a common choice is: where x is the estimated drive at the last iteration,and p is the current search direction.The vector ey is the error in the response to the input x + λp.
The idea behind a backward line search is to model the merit function using a polynomial -typically a quadratic or a cubic.Quadratic interpolation is employed in this paper.
In the simulations presented later, the merit function of equation( 33) is used to start off the process, i.e. to compute λ 0 , and hence x (0) .A search is made for a value of λ that minimises the merit function.As an illustration of the behaviour of the error, φ is plotted as a function of λ for k = 0.30 in Figure 7.
The alternative merit function, equation(34) is employed during the iteration.
Again, for illustration, Figure 8 is a plot of ψ as a function of λ for the case of k = 0.3 in the fourth iteration of Broyden's method.For a full explanation of these and other search algorithms the reader is referred to Dennis et al [9], The alternative methods of iteration are now examined.The validity of the mathematical methods is established and their performance is compared.There are alternative choices for parameters and so, in the simulations presented here, each one uses the same Duffing model, the same desired solution and the same starting point.The same seven levels of non-linearity are chosen for each alternative iteration method, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40 and 0.45.
In each case, the iteration is stopped when the response is within 5% of the target response, or after a specified number of iterations.
For reasons of clarity of presentation, only four cases are plotted.However, all cases are represented in the tables.
The first method to be tested is the basic conventional algorithm as shown in Figure 9.The graphs figure shows that the conventional algorithm fails to converge at levels of non-linearity 0.25 and greater.In industrial practice, the engineer would reduce the iteration gain at the expense of running the system more often and also would examine the spectral densities in which troublesome frequencies may be detected.
Figure 10 and Table 1 demonstrate the validity of manually reducing the An early conclusion of considering the algorithm in the context of solving a system of algebraic equations, was that an appropriate iteration gain might be computed from the progress of the iteration, using a search for an appropriate iteration gain.This has been implemented and the results are presented in Figure 11 and Table 2.The divergent behaviour of Figures 9 and 10 are eliminated, although manually choosing a small iteration gain initially performs better.However, even at a reduced gain, the convergence of the iteration stagnates for k = 0.25 and 0.45, with a response error of about 11% after 35 iterations (about 100 system runs).The same behaviour is also observed for k = 0.30 as indicated in Table 2.
The above results used the conventional algorithm, and the conventional algorithm with search.These use a constant approximation to the Jacobian.The effect of updating the approximation, using Broyden's method, is now considered.Figure 12 illustrates the results of this approach without a search.
The method works well for low levels of non-linearity where it produces faster convergence.There is also convergence for levels of non linearity that failed to converge using conventional iteration.At higher levels of non-linearity, the method still fails to converge.The progress of Broyden's method improves when the iteration gain is reduced, Figure 13, but has no great advantage over the conventional method as measured by the number of runs required to achieve a tolerance of 5%.
However, it was noted that, as the number of runs were increased the error dropped faster for the updated technique -as indicated in Table 1.
The last of the sets of simulations presents, in Figure 14, Broyden's method with a search.The method is successful in achieving convergence at all the levels of non-linearity that were considered but at the expense of running the system more often.tion to the Jacobian, by showing a faster convergence, and also by achieving convergence when the conventional algorithm fails.
The performance of search routines depends on chosen parameters and this requires further study.For example, the work presented does not use restarts, nor does it consider the effect of the many different forms of line search.In addition, there are many other types of update -including updating the frequency response function itself -for others see e.g.Martinez [6].Other approaches take advantage of the particular structure of the Jacobian.The aim of this work is to indicate that the particular point of view presented can be advantageous, and that recourse can be made to an arsenal of tried and tested techniques.

Conclusions
A new mathematical framework for the derivation of drive files for laboratory simulation test systems is demonstrated.The conventional algorithm is shown to be part of a broad mathematical area for which established mathematical techniques are available.This approach can achieve convergence in systems that do not readily converge with the conventional algorithm.It has also been shown that there is potential for improving the convergence of the latter using a backward line search.
Thus, this paper reports on a beginning, not a completion, of an investigation.
The authors regard the work as the opening up of an area for further research.
Consideration will be given to other solution techniques and the sensitivity of these to measurement errors.Systems with multiple channels, physical systems

Fig. 1 .m
Fig. 1.Current iteration algorithm -setting x (0) := 0. comprised of a mass m, a viscous damper with coefficient c, and a non-linear spring.The stiffness of the spring increases with amplitude as described by a linear stiffness coefficient k, and a non-linear factor kk .Such systems are usually described in terms of natural frequency 1 2π

Figure ? ?
Figure ?? illustrates the magnitude of the measured frequency response function -averaged over ten records -and a smoothed version obtained from a least squares fit to produce a rational function which has as numerator a linear polynomial and as denominator a quadratic polynomial in frequency.In addition, the frequency response function corresponding to the linear part of equation(4) is also shown for comparison.

Fig. 4 .
Fig. 4. Response error against time using the current algorithm, for k = 0.2.

Fig. 9 .
Fig. 9. Progress of conventional iteration for various levels of non linearity and an iteration gain of 0.5

Fig. 10 .
Fig. 10.Progress of conventional iteration for various levels of non linearity and an iteration gain of 0.2 iteration gain to 0.2.Convergence is achieved, for all bar the most non-linear case, at the expense of a slower rate of convergence.The choice of iteration gain is generally left to the experience of the operators of industrial systems.

Fig. 11 .
Fig. 11.Progress of iteration for various levels of non linearity -conventional algorithm with search

Fig. 12 .
Fig. 12. Progress of iteration for various levels of non linearity -Broyden's method with an iteration gain of 0.5

Fig. 13 .
Fig. 13.Progress of iteration for various levels of non linearity -Broyden's method with an iteration gain of 0.2

Fig. 14 .
Fig. 14.Progress of iteration for various levels of non linearity -Broyden's method with search.For k = 0.15 the iterations are shown by '+', for k = 0.25 by a dash-dot line , for k = 0.35 by ' * ' , and for k = 0.45 as a dashed line.

Table 2
compares the results for the conventional algorithm(constant Jacobian approximation) and Broyden's method(updated Jacobian approximation) using backward line searches for all the cases of non-linearity.It indicates that the Broyden update has an advantage over the use of a constant approxima-

Table 1
Table of percentage error in the response for various levels of non-linearity, for 35 runs, with λ = 0.2.