Dynamical system analyser

IS 2 149 172 T3

DESCRIPTION

Dynamic systems analyzer.

This invention relates to a dynamic systems analyzer, and more particularly to a device of this type applicable for the analysis of dynaomic systems that can be non-linear or chaotic in the mathematical sense.

The investigation of nonlinear dynaomic systems presents particular difficulties because the traditional tools of electronic spectral analysis are not appropriate. As an example, consider a chaotic nonlinear dynamical system arranged to produce an electronic signal characteristic of its behavior. It is known from research into dynaomic systems that such a signal can have a substantial continuous component in its spectrum. Filtering, either in time or frequency to improve signal-to-noise relationships can distort or alter the dynamics of the system as perceived from a filtered signal. This is exemplified in the writings of Badii and others in Phys. Rev. Lett. 60 (1988) page 979 and Milschske et al. In Phys. Rev. A37 (1988) page 4518. Furthermore, the power of spectrum analysis is insufficient to characterize the dynamics of a system when the data show a deterministic chaos that results in a continuous spectrum (as opposed to a set of harmonics with a simple relationship). This was covered by FCMoon in a paper entitled "Chaotic Vibrations", Wiley-Interscience. Spectrally selective signal processing of this type intended to simplify analysis can make the apparent behavior of the system more complex.

The behavior of a dynamic system is normally represented by a curve in a multidimensional phase space. Successive points on the curve indicate the evolution of the system as a function of time. The line is referred to as a "trajectory", and the region of phase space to which it is confined after an arbitrary time is called an "attractor". If the dynamical system is chaotic, the region is called a "foreign attractor". The attractor is the region of the phase space in which the behavior of the system was located and towards which it can be said that it is "attracted".

A device is needed that is applicable for the analysis of non-linear systems, since traditional techniques simply complicate the apparent behavior rather than reducing it to identifiable components. Traditional techniques, such as Fourier spectral analysis, are applicable only to systems that can be modeled as the linear superposition of a limited number of harmonic modes.

There is a particular need for a device that can detect a change from normal system behavior to chaotic behavior. This is of particular importance in the field of aircraft engines and other sources of mechaonic motive power. Such energy sources are characterized by very regular oscillatory and / or rotational behavior in normal operation, but show irregular characteristics before catastrophically failing. A device capable of detecting the onset of uneven wear-related behavior would provide a means to anticipate and avoid catastrophic failure by prior closure and repair. It could also provide a reduction in service costs, since maintenance could be deferred until justified by the behavior of the system.

It is an object of the present invention to provide a dynaomic systems analyzer.

The present invention provides a dynaamic systems analyzer and a dynaamic systems analysis method as set forth in claim 1.

The invention provides the advantage that it produces a mathematic model of a dynaomic system, and is applicable to the analysis of non-linear systems. As will be described later, a mathematical model created by heuristic processing means can be used to generate predictions of the future behavior of the dynamical system based on the preceding behavior. Alternatively, mathematical models derived from different dynamical systems can be compared with each other in order to analyze similarity or variation.

The heuristic processing means can also be arranged to operate in a test mode to transform another sequence of phase space coordinates and provide predictable output values ​​based on the mathematical model.

The invention provides the advantage of being trainable according to reference data to generate a mathematical model. This allows the reference data to be selected and also to be replaced if retraining is required.

The heuristic processing means may be arranged to employ in the training mode reference data which is one or both of the time series and the sequence of phase space coordinate sets derived from it. The reference data is entered into the heuristic processing means with a time interval ahead of the input of said sequence for the QR decomposition and the adjustment by square mononyms; and the lead time interval is arranged not to be greater than the predictable time scale inherent in the dynamic system as indicated by its maximum Lyapunov component. In this case, the mathematic model that is generated in the training mode is the one required to extrapolate the trajectory of the phase space and / or the time series of the dynamical system in advance of time to the extension of the lead time interval. In test mode, the output values ​​correspond to predictions of the future behavior of a test system that produced the other sequence of phase space coordinates referenced above. These values ​​can be the coordinates of the phase space and / or the time series of said test system predicted based on the mathematical model obtained in the training mode.

IS 2 149 172 T3

Alternatively, or additionally, the heuristic processing means may be arranged to compare the predictions of the mathematical model of the phase space coordinates and / or the time series of the test system with actual values ​​of these parameters measured later in time. In this case, error values ​​are presented between each prediction and a respective measured parameter. A large error value indicates a significant deviation from a predicted extrapolation of a phase space trajectory. In a mechanical system, this could be associated with impending catastrophic failure or other untimely behavior changes.

The generation means can be a digital computer programmed to carry out a singular value decomposition of the time series of the dynamical system. It can be a one-time-for-all computer, in which case your results can be stored in a memory for use when needed.

The dynamic system can be a system (such as a mechanical or electronic device) in a functionally acceptable initial state. A test system to be compared to the dynamic system can be the same system after deterioration in service has occurred over a period of time. This would allow maintenance requirements on large machines (for example) to be evaluated based on behavior change.

The transformation means can be a finite impulse response (FIR) filtering device arranged to introduce the time series in the form of Takens vectors (to be defined later) and calculate the projections of each Takens vector on each of the a set of singular vectors obtained by decomposing into singular values ​​to provide a respective set of coordinates in a phase space for each such vector of Takens.

The heuristic processing means are preferably substantially as described in international patent application No. PCT / GB90 / 00142 published under the Patent Cooperation Treaty as No.^° WO 90/09643.

The dynamical systems analyzer of the invention can be employed with the heuristic processing means operating in training mode to provide mathematical models of respective dynamical systems. The models can then be compared to provide an indication of the degree of similarity between the systems. The singular vectors obtained from a series of times of the dynamic system can be used in the production of coordinates of the phase space for each of the time series of the systems. Alternatively, each time series of the system can be used in singular value decomposition and transformation by heuristic processing means.

In order for the invention to be understood more fully, an embodiment thereof would now be described, with reference to the attached drawings, in which:

Fig. 1 is a schematic block diagram of a dynamical systems analyzer of the invention in the form of a prediction device;

Fig. 1a shows a part of Fig. 1 with my detail;

Fig. 2 is a graph of a series of signal times for the input to the device of the

Fig. 1;

Fig. 3 is a graph of an autocorrelation function of the time series of Fig.

2;

Fig. 4 provides the results of a singular value decomposition of the time series of Fig. 2 and indicates a system noise level above digital quantization noise;

Fig. 5 is equivalent to Fig. 4 except that the relevant noise level of the system is below the digital quantization noise;

Fig. 6 schematically illustrates an alternative embodiment of the invention arranged to determine the predictability of the system; Y

Fig. 7 is a graph showing our results of the kind that could be obtained from the realization of Fig. 6.

With reference to Figures 1 and 1a, a schematic block diagram of a dynamical systems analyzer device of the invention is shown generally indicated 10.

The device 10 incorporates a transducer sensor 12 mounted on a non-linear system 14 and which derives an analog electrical signal A (t) characteristic of the behavior of said system. The output signals from sensor 12 are fed to an analog-to-digital converter (ADC) 16. This has an S / H input arranged to sample and hold the signal A (t) at successive times t1 ... ti ... tM with a constant spacing τ. Here you are given by:

ti = ti + (í-1) t (1)

The analog signal A (t) and the sampling times ti etc. are shown in a graph 18. Discrete samples of the signal form a time series A (ti), where i = 1 to M. These are digitized by the ADC 16 to produce a digital time series V (ti), where V (ti) corresponds to A (ti) in each case for i = 1 to M.

The ADC operates under the control of a sampling clock 20 that supplies timing pulses at separate times τ. The magnitude τ is adjusted by a computer 22 which supplies a digital word to the clock 20 indicating the required sampling interval.

The digital time series with the general term V (ti) is output from the ADC 16 to the computer 22 in parallel bit form via an input bus 24 to the computer. Computer 22 stores this series

ES 2 149 172 T3 of times on a hard disk, and calculates the mean and variance of the terms of the time series V (t_i) to V (t_M). Produces a normalized version of the time series, and then processes it to provide singular value decomposition. The decomposition will be described in more detail later. Produces a set of singular vectors f1 to fd, each of which is a mathematical function expressed as a series of n successive values. Some of the vectors are schematically illustrated in boxes such as 28, and their corresponding singular value spectrum is schematically indicated at 29.

The computer 22 loads the singular vectors f1 afd into a random access memory (RAM) 30 via a first output bus 32. The RAM 30 has respective portions 301 a30d for the vectors f1 afd. Each portion, such as 301, stores n values ​​per vector, for example f11 to f1n for f1. RAM 30 is incorporated into a multidimensional finite impulse response (FIR) filter 34, which also incorporates a chain 36 of n multi-bit delay latches 361 to 36n. Parts of RAM 30 and latch chain 36 that are not shown They are explicitly indicated by dotted lines. The FIR filter 34 is of a generally known and commercially available type. One such incorporating a RAM 30 and a chain of bolts 36 was offered for sale in May 1985 by Calmes Systems Inc., a Canadian company. The FIR filter 34 may alternatively be based on the devices set forth in US Patent Nos. 4,885,715 and 4,833,635 referred to as correlators and convolutioners.

Lock chain 36 receives input of stored digital time series values, such as V (ti), from computer 22 to a branch 38a of a second output bus 38. Bus 38 also supplies a clock signal. to the FIR filter 34. This makes successive values ​​V (ti) etc. in the series of digital times they are timed in a final lock 36n of the chain 36. Each value is a digital word of multiple bits in parallel and its bits are inserted synchronously. After entry, each time series value is progressively timed through the chain to a first latch 361 at the rate of one latch per clock cycle. After n clock cycles, the digital time series portion V (t1) to V (tn) is stored in locks 361 to 36n, respectively. This is known in electronics as a record pre-fill operation. After the pre-fill, the computer 22 supplies a clock signal drive to the FIR filter 34 to initiate its filtering operation. The operation of filter 34 is to multiply each time series value by a respective horizontally opposite singular vector in Fig. 1 to form products such as V (ti) f11 through V (ti + n-1) f1n. The products formed in each clock cycle are added to provide a respective total for each singular vector.

Defining the first filter clock cycle as the one during which the digital time series signal or element V (t1) is loaded into the FIR filter 34, then the filter outputs on the i-th filter clock cycle coordinates g1 (i) agd (i) given by:

n g1 (i) = V (ti + j-1) f1j (2.1) j = 1 n

gk (i) = V (ti + j-1) fkj (2.k) j = 1 n

gd (i) = V (ti + j-1) fdj (2.d) j = 1

Equations (2.1) to (2.d) indicate that the ith digital time series value V (ti) is used by the FIR filter 34 to form a set of coordinates g1 (i) to gd (i). These coordinates are produced through the subset of the time series of n successive values ​​V (ti) to V (ti + n-1), and represent the position in time ti of a vector that represents the non-linear system in a phase space having dimension d, defined by the respective base vectors f1 to fd The general coordinate gk (i) is the projection of the time series subset V (ti) to V (ti + n-1) on the base vector general (singular vector) fk. The subset of n values ​​V (ti) to V (ti + n-1) is referred to as Takens vector, by F. Takens, "Detecting Stange Attractors in Turbulence", Lecture Notes in Mathematics, Ed.

DA Rand and LS Young, Springer, Berlin 1981 page 266. The ith Takens vector has an ith element T¿ (i) defined by:

T, (i) = V (ti +, _ x) (3)

Varying the element number é from 1 to n generates V (ti) to V (ti + n-1). The general coordinate gk (i) produced in the i-th filter clock cycle can be considered as the correlation of the i-th Takens vector with the sequence of values ​​fk1 to fkn of the k-th singular vector. Mathematically it is the projection of the i-th Vector of Takens on the k-th singular vector.

The coordinates g1 (i) to gd (i) output in parallel and in synchronism from the FIR filter 34 on a filter bus 40. Successive sets of coordinates g1 (i + 1) to gd (i + 1), g1 (i + 2) to gd (i + 2) etc. they come out in successive clock cycles. The filter bus 40 has a first branch 40a connected to a multipole, two-position switch 42, and thus to an AI response input of a heuristic processor indicated generally by 44. Here, the legend "HEUR PROCESSOR^andISTICO "indicates elements AI, DI and 50 to 68b described later. The positions of switch 42 connect the response input AI, either to the first branch 40a of the bus, or to a zero signal (not shown). A second branch 40b of the filter bus 40 is connected to a delay unit 46, which delays each of the coefficients g1 (i) to gd (i) in equal time intervals of pT. Here p is an integer and τ is the length of a clock cycle, so the delay is the same whole number of clock cycles for each of the coordinates g1 (i) to gd (i).

The delay unit 46 is connected via a delay output bus 48 to a bus

ES 2 149 172 T3 data input DI of the heuristic processor 44. A second branch 38b of the second output bus of the computer 38 is connected via the switch 42 to a subsidiary response input bus SAI of the heuristic processor 44. An interface of Processor input 44 is indicated by a dotted line 50.

The heuristic processor 44 is of a known type. A detailed description of its construction and mode of operation is given in the international patent application No. PCT / GB90 / 00142 published under the Patent Cooperation Treaty (PCT) on August 23, 1990 with an international publication no.^° WO 90/09643. Its characteristics will be briefly discussed.

The heuristic processor 44 includes a digital arithmetic unit 52. This incorporates arithmetic circuits P11 aPdr shown as straight angles. These were arranged in rows P11 to P1r, ... Pd1 aPdr and columns P11 sPd1 ... P1r to Pdr. The parameter d is the dimension of the coordinates g1 (i) etc. like before. The parameter r is the dimension number of the space in which the input data to unit 52 is transformed. Unit 52 also incorporates memories M1 ... Mr shown as diamonds. These are each arranged to receive output signals from a respective column of arithmetic units, such as P11 aPd1 for memory M1, for example. Devices of class P or M that are not explicitly illustrated are denoted by dot-and-dash lines.

The output of the digital arithmetic unit goes to a triangular systolic matrix 54 which is also part of the heuristic processor 44. The triangular matrix 54 is shown vertically in foreshortening for the convenience of illustration. It incorporates border cells B11 aBrr indicated by circles, along with inner cells, I12 to Ir-1, r, indicated by squares. Structure not explicitly illustrated is denoted by dashed lines. Each boundary cell (generally referred to as B) of the triangular array 54 is arranged to compute rotational parameters of the data input from above, to update a respective stored array element, and to output parameters laterally to the right. Each inner cell (generally referred to as I) of the triangular array 54 is arranged to receive rotational parameters from the left, to apply them to data received from above, to generate an output below, to update an element. of the respective stored matrix, and pass in the rotation parameters to the right. These functions take place in each clock cycle, with the heuristic processor 44 receiving a clock signal from the computer 22 via a line 56. This clock signal is transmitted to each element of the heuristic processor 44.

Boundary and internal cells B and I are transputers programmed to perform arithmetic functions. They are designed to implement a rotation algorithm. A transformed data matrix exits the digital arithmetic unit 52 at a rate of one row matrix per clock cycle. This is fed to the triangular array 54 at the same speed, and the elements in the rows are progressively removed by a series of rotations according to the algorithm at the speed of one element per row. The result is that the matrix of data transformed from unit 52 undergoes rotations that make it a triangular matrix. There are various rotation algorithms that the triangular matrix 54 could use. One such uses pairs of rotation parameters consisting of sines and cosines of the respective angles of rotation. The boundary and inner cells B and I store and update the elements of the triangular matrix. Other more computationally convenient algorithms produce rotation parameters and stored items mathematically related to such sines, cosines, and matrix items. In particular, in an algorithm, border cells B store a matrix of a first class and internal cells I store a matrix of a second type intended to multiply the matrix of the first class to generate the triangular matrix.

The process of triangularization of a matrix by a series of rotations in the plane is known as QR decomposition, where Q is a matrix of rotation parameters and R is the triangular matrix. Triangular matrix processing functions for implementing rotation algorithms are well known in the prior art, such as the international application referenced above, and will not be described by me.

Rotation parameters pass along rows of cells such as the top row B11 aI1r. These then pass into a rectangular matrix 58 of internal cells I1, r + 1 aIr-1 and multiplier cells Mr + 1, r + 1 to Mr + 1, r + d + 1 indicated by hexaigons. Boundary cells B11 aBrr have the additional function of cumulatively multiplying their rotational cosine parameters (or cosine type parameters in other algorithms). For example, when the first boundary cell B11 generates sine and cosine parameters in one clock cycle and outputs them to the inner cell to the right, it also outputs the cosine parameter diagonally below and to the right to the second row boundary cell. , second column (not shown, it would be designated B22). This cell multiplies by its cosine parameter generated two cycles later, with a delay of two clock cycles between a diagonal output from a border cell and an equal output from its diagonal adjoining cell that incorporates the first such cell output. This diagonal boundary cell function generates products that arise from cosine parameters multiplied cumulatively, each product generated from a respective input row of the matrix that comes out of the digital arithmetic unit 52. Each product is fed to the row of multiplying cells Mr +1 to Mr + 1, r + d + 1, and passes through it at a rate of one multiplier cell per clock cycle. Similarly, the rotation parameters passing along rows of the triangular matrix 54 are fed into respective rows of the triangular matrix 58. These parameters are used by the inner cells I1, r + 1 etc. from the rectangular array to rotate received data from above, to update the respective stored items, output

ES 2 149 172 T3 below and pass the parameters to the right; that is, the internal cells of the rectangular matrix operate in precisely the same way as those of the triangular matrix 54.

Multiplying cells (generally referred to as M) have the processing function of receiving data from above, multiplying it by an input from the left, pulling the product from below, and passing the input from the left to the right. The product is the scalar least squared residual error obtained by fitting the transformed data from the arithmetic unit 52 to the data that passes down the column that exceeds the relevant multiplier cell in each case. The sets of residuals emerging from the rectangular matrix 58 form the respective least squares residual vectors.

The rectangular matrix 58 is imaginatively partitioned into a multi-column portion 58a and a one-column portion 58b.

The triangular matrix 54 and the rectangular matrix 58 receive different input data. The first receives data from the DI input that has undergone transformation in the digital arithmetic unit 52, as has been said. The rectangular array receives input data from the AI ​​response input via lock chains (not explicitly illustrated). These are indicated as an L-shaped delay device 60 incorporated in the heurastic processor 44; Device 60 is imaginatively partitioned into multichannel portion 60a and uan channel portion 60b respectively connected to portions 58a and 58b of rectangular array 58.

The data input bus DI is connected to the arithmetic unit 52 via a time staggering device indicated by a triangle 62. This device is a triangular array of timed locks (not shown). Introduces a time delay that varies in magnitude across the DI input bus. The delay is zero at the upper edge of the DI input bus, increasing by one clock interval τ per bus segment (not shown) to (d-1) r at the lower edge of the bus. Similarly, a second time stepping device 64 introduces a similar stepping delay variation at 64a on the AI ​​response input bus along with a dT delay at 64b on the subsidiary response input bus SAI. The outlet of the rectangular die 58 passes through a temporary de-stepping device 66 similar to device 62. The first introduces a rectangular matrix output delay that varies from dT for the column of cells I1, r + 1 to Mr + 1, r + 1 to zero in the column furthest to the right of cells I1, r + d + 1 a Mr + 1, r + d + 1. The delay is reduced by one clock cycle interval T per column.

The purpose of input delay devices 62 and 64 is to ensure correct timing of the arrival of data within elements 52, 54, 58a and 58b of the heuraistic processor regardless of the data path. The de-stepping delay device 66 provides synchronous input of data at DI, AI, and SAI to result in saynchronous outputs at the output interface indicated by a line 68.

Line 68 has portions 68a and 68b corresponding to outlets of rectangular matrix portions 58a and 58b, respectively. The use of timed latch triangular arrays, such as 62, 64, and 66 is well known in the art of systolic array processing and will not be described further.

The heuraistic processor 44 has two alternative modes of operation selectable by means of a signal applied to a Tr / Te input. A high signal (1 binary) implements Tr (training) mode, and a low signal (0 binary) implements Te (test) mode. This mode selection signal is passed to the triangular and rectangular matrices 54 and 58 via an appropriate delay (not shown) within the arithmetic unit 52 via a connection 70.

In the training mode, the training data is derived from the non-linear system 14 via the computer 22 and the FIR filter 30 as described above, and is fed to the DI data input bus after the delay at 46. Version no. The delay of the training data is fed to the AI ​​response input, with switch 42 being set to connect buses 40a and AI together. The lagged data undergoes a transformation in the digital arithmetic unit 52. This involves the computation of the Euclidean norm (geometrical distance) of each data vector from each of a set of centers (geometrical origins). Here, each data vector is a respective set of phase space coordinates. Each of these norms then undergoes a non-linear transformation. The transformation results then undergo a QR decomposition in the triangular matrix 54. This provides the rotation parameters for the application to non-delayed data (advanced by p clock cycles) in the rectangular matrix 58 received via the delay device 60. In training mode, internal and boundary cells B11 to Ir, r + d + 1 all operate adaptively; that is, when generating and applying rotation parameters, respectively, they update their stored matrix elements in each clock cycle according to each new set of inputs.

The result of operating in the training mode is that the data delayed by p clock cycles (interval pT) is transformed and least squares fitted to the non-delayed data. The precision of the fit is indicated by the residuals (errors) that appear at the output of the heuraistic processor 68. An acceptable fine fit is indicated by small fluctuating residuals.

At the end of the training, the internal and boundary cells B11 a Ir, r + d + 1 have respective matrix elements stored inside them. In the first portion of the rectangular matrix (several columns) 58a, the elements correspond to a point having coordinates gk (i) (k = 1 to d) on a trajectory in a phase space of dimension d one point p stages later in said trajectory that has coordinates gk (i + p). In the second portion of the rectangular matrix (a single column) 58b, the matrix elements correspond to the point gk (i) (k = 1ad) unva6

ES 2 lor of digitized time series V (ti + n-1 + p) also p advanced stages.

When the training mode is complete, the Tr / Te control mode input is set to a low voltage signal level (binary 0) to implement the test mode. This is transmitted to the triangular and rectangular matrices 54 and 58 via the connection 70. It connects all the boundary and internal cells of the classes B and I in these matrices to modes in which the update of the internal matrix mode is suppressed. The evaluation and application of the rotation parameters in cells B and I respectively continues, however. Consequently, the mathematical model expressed by the matrix elements stored during training is frozen.

There are two alternative test modes of operation depending on the position of the switch 42. In both modes, the "known" or "comparison" non-linear system 14 is replaced by an "unknown" non-linear system required to be compared with the First. The first mode implies that the switch 42 remains in the training position, and thus connects the branch 40a of the bus to the response input AI. The data derived from the unknown system are then transformed into coordinates gk '(i) etc. by the FIR filter 34 using singular vectors f1 afd derived from the known system 14.

The coordinates gk '(i) etc. they are fed directly to the AI ​​response input due to the position of switch 42. They are also fed via delay 46 to the DI data input. The associated time series values ​​V '(ti) are fed to the subsidiary response input SAI. The input in DI is nonlinearly transformed into an r-dimensional space by the arithmetic unit 52. The output of unit 52 is further transformed by the frozen mathematical model stored in the triangular and rectangular matrices 54 and 58. After this other transformation it is subtracted in the rectangular matrix 58 from the coordinates and time series values ​​received from the device. delay 60. This results in error values ​​at the output of processor 68. Each error value is the difference between a respective measured value of the unknown non-linear system and an associated predicted value calculated by the heurastic processor 44 of its frozen mathematical model. If the error values ​​are small enough, it is shown that the unknown non-linear system is equivalent to the known system 14, because its future behavior is predictable p stages earlier based on the behavior of the known system. The predictability will exist in phase space (rectangular matrix portion 58a) and in the corresponding time series (rectangular matrix portion 58b). This is a demonstration that signals from a test system can be used to detect and quantify their similarity to those from a previously observed "standard" system.

The first test mode of operation is also important for the purpose of detecting changes in the functioning of systems that may be chaoatic but are deterministic in the field.

172 T3 12 mathematical sense. For example, an oscillatory system such as a noisy mechanical bearing will affect the readings of a transducer arranged to monitor a parameter of the device of which the bearing is a part. In this case, the "known" nonlinear system 14 used in training will be the device with the bearing in good working order. The output of the transducer (containing the bearing noise) is also the "unknown" system. The “error values” that come out at 68 during test mode will be small as long as the bearing is operating normally. The error values ​​will increase as the bearing degenerates from wear. This is a demonstration that the system signals are separable according to the invention when they are made different.

The second test mode of operation requires that switch 42 be set to its second position. This disconnects branch 40a from the bus and applies a binary zero signal level to each segment of the AI ​​response input. In other respects, the operation in the second test mode is as in the first. Consequently, the input of data at DI and transforms at 52, 54, and 58 is compared to the input of zero signals to the rectangular matrix 58 from the delay device 60. Consequently, the error values ​​that appear at the output interface 68 are in fact generated from comparisons with zero signals. These values ​​are therefore estimates of the coordinates in phase space (output 68a) and the time series (output 68b) of the "unknown" system p stages (clock cycles) in the future in each case. The predictions are made based on the multidimensional matrix model formed during the training of the triangular and rectangular matrices 54 and 58 and stored inside. This shows that the prediction device 10 has the ability to generate predictions of the future states of a non-linear dynamic system.

The mode of operation of the computer 22 will now be described in more detail. Computer 22 is programmed to execute a sequence of operations; this produces a singular value decomposition of the time series data of the type V (ti) of the "known" system 14. The sequence is as follows:

(a) determine an acceptable value for n, the dimension of the phase space;

(b) carry out a singular value decomposition of the time series V (ti) etc. to define a basis for the phase space;

(c) determine singular noise-corrupted functions generated in (b) and omit them to produce a subspace of dimension d;

(d) load singular vectors into filter

FIR 34;

(e) activate the FIR filter 34 by the clock signal to read in V (ti) etc; Y

ES 2 149 172 T3 (f) activate the heuristic processor 44 n clock cycles after (e) to start the input of coordinates and signals of the class gk (i) and V (ti).

To carry out the above operations, the computer 22 reads in the digitized time series that has the general term V (t_i), where i = 1 to M. A typical value for M would be 50,000. Convert the terms of the time series into a set of normalized values ​​by calculating their statistical mean and variance, subtracting the mean of each term, and dividing the result by the variance. For convenience, a normalized time series value derived from and equivalent to a sampled ADC output signal V (ti) will hereinafter be referred to as V (ti).

For the purposes of testing a "known" non-linear system 14 was built. This consisted of a mechanical bearing in the form of a ball channel, an aluminum bar contacting the center of the channel, and a positional transducer contacting the bar. The transducer was a piezoelectric crystal needle of the type used to establish record contact in home turntables. It gave an electrical output indicating bearing noise, that is, unwanted non-linear bearing vibration and at least partially deterministic.

The transducer produced a time series of 50,000 (ie M) values. The first 2048 of them are shown in Fig. 2 in digitized form after normalization by the computer.

22.

To carry out operation (a) (as defined above) and determine an acceptable number n for the dimension of the phase space, the computer 22 calculates the autocorrelation function of the time series V (ti) (now normalized) etc. If Vi is written instead of V (ti), the general term for the normalized time series to reduce the complexity of equations, the autocorrelation function has successive values ​​of type Gk defined by:

Ti = [Vi, Vi + 1, ..., Vi + n-1 (5)

Equation (5) shows that the i-th Takens vector Ti (i = 1 to M-n + 1) has vector elements that are a sequence of n values ​​of the successive time series beginning with Vi. This was explained by F Takens in "Detecting Stange Attractor in Turbulence", Lectures Notes in Mathematics, Ed. DA Rand and L.-S. Young, Springer, Berlin, 1981 page 366. X is therefore given by

T1 ^t2

V1

V2

V2

V3

Vn

Vn + 1

X = (6) ^TM-n + 1 ^VM -n + 1 ^.^VM

The transpose of X is X^T given by:

X^T = [T1^T ... T^TM-n + 1] (7) that is,

The autocovariance matrix Ξ is defined

<td colspan="6">for:</td>

<td></td><td>_Ξ</td><td>= X^T X</td><td> (9)</td><td></td><td></td>

<td></td><td>V1</td><td>V2 ..</td><td>^..VM -n + 1</td><td></td><td></td>

<td>XT =</td><td></td><td></td><td></td><td></td><td> (8)</td>

<td></td><td>Vn</td><td>Vn + 1 ..</td><td>.. VM</td><td></td><td></td>

MK

Gk = (MK)^-1 ViVi + k (4) i = 1 where k = 0, 1, 2, 3, ... (K-1), where K is the maximum value of k, and K is less than M. 50

Figure 3 shows the autocorrelation function G for the normalized time series data shown in Fig. 2. Computer 22 calculated the function G for k values ​​from 1 to 128, and the Gk values ​​are shown plotted 55 in units arbitrary with respect to k in Fig. 3.

The autocorrelation function crosses the k axis, that is, Gk = 0, at approximately k =

28. It has been found empirically that an appropriate value of n, the dimension of the phase space, 60 is equal to the first value of k where Gk = 0. The value of n chosen is therefore 28.

Computer 22 now implements operation (b), the singular value decomposition of the time series. To do this, form 65 the matrix X in which the successive rows are the successive vectors of Takens Ti etc. These latter defined by

Ξ is then processed by computer 22 to provide the required singular value decomposition. This matrix is ​​subjected to the well-known mathematical diagonalization process. The process consists of the Householder transformation of the matrix to the tridiagonal form, followed by the application of the QR algorithm to obtain the eigenvalues. The eigenvalues ​​are the squares of the required singular values. This was described in WH Press et al., In "Numerical Recipes - The Art of Scientific Computing" Cambridge, page 350 et seq. Also, the decomposition produces a set of n singular vectors (n = 28 in the present case) expressed as sets of digital values ​​of defining mathematical functions f1 to fn. Each singular vector is associated with a respective singular value s1 asn, which could be loosely referred to as the force of the respective vector. The decomposition could be considered as the generation of the higher dimensional equivalent of an ellipsoid that has axes directed in the directions of the respective singular vectors and axial lengths given by the respective singular values. This is explained by WH Press and others (as referenced above) on pages 52-64.

The lengths or forces of the singular vectors, that is, their respective singular values,

ES 2 149 172 T3 are arranged to decrease monotonously as j increases, j being the index number of the vector at fj and j = 1an. Normally I will drop to the noise level of the nonlinear system before j reaches n. This indicates that the singular vectors associated with singular values ​​that are small are significantly corrupted by noise. Therefore, such vectors should be omitted. Consequently, the phase space dimension is reduced from nad, where the singular values ​​si asd are greater in magnitude than the noise level.

In a properly designed analyzer 10, the ADC 16 would be arranged to make the noise level approximately equal to but greater than the least significant bit of the digital word exiting the ADC.

Programming a computer to deal with equations (4) through (9) to truncate the singular value decomposition and properly supply timed data and clock signals is straightforward and would not be described. Singular value decomposition in the preceding context was described by Broomhead et al. In Physica 20D, 1986, page 217.

Referring now to Fig. 4, the results of a singular value decomposition of the 50,000 normalized digital time series values ​​obtained from the bearing ball channel referred to above are shown. A graph indicated generally by 80 provides the spectrum of the singular values ​​if as produced in the decomposition. The present example n = 28, so the values ​​are si as28, with a general value si, where i = 1 a

28. Graph 80 shows a decreasing monotonic curve 82 of logi0 si with respect to the singular value number i (or singular vector) with i varying from 1 to 28. At i = 1, curve 82 has a logarithm of approximately 0.6, corresponding at a relative magnitude of approximately 4. At i = 28, it also decreases to a logarithmic value of approximately -2.3, equivalent to approximately 5 x 10^-3. It does not cross the noise level of the ADC 16 (indicated by the dotted line 84). However, as indicated at a point 85 where horizontal and vertical dashed lines intersect and 87 kink appears in curve 82. To the left of vertical line 87, curve 82 varies rapidly with the singular value number i ; to the right of this line the variation of logi0 si with i is insignificant. The vertical line is at i = 7. For the same measurements, it can be inferred that the nineteen singular values ​​were corrupted by the noise produced in the system 14 and the transducer sensor 12. This inference can be deduced because the noise is homogeneously distributed, and the lack of variation in a spectrum is associated normally at a noise level. The horizontal dashed line 86 is therefore treated as the system noise level (as opposed to the ADC quantization noise level 84). Singular vectors for which i> 7 are therefore considered unacceptably corrupted by noise. Singular vectors or functions fi to f7 have an acceptable degree of noise. Therefore, these are retained and f8 to f28 are discarded. This creates a reduced noise subspace of the non-linear system, with said subspace having dimension 7 in the example. The value of the parameter d in equation (2.d) is therefore 7.

Fig. 4 also shows the twenty-eight singular vectors generated by the singular value decomposition in the present example, the vectors being shown as curves labeled fi to f28.

Comparing Fig. 4 with classical Fourier analysis, it is noted that this assumes infinitely long (and therefore impractical) time series. Fourier analysis characterizes linear systems with discrete frequencies (line spectrum). It is not useful for non-linear systems displaying broadband spectra.

If necessary, it is possible to combine the alternatives provided by switch 42 in Fig. 1. This requires switch 42 to be removed and a third portion 58 (not shown) added to rectangular die 58. Portion 58c would be of a construction similar to that of portions 58a and 58b combined. In training mode, portion 58c would receive AI input, which would be branched for this purpose. In test mode, portion 58c would receive zero inputs. A third output (say 68c) would then provide the predicted values ​​of the coordinates and time series while the errors between these values ​​and the measured parameters would appear in 68a and 68b.

Referring now to Fig. 5, results similar to those shown in Fig. 4 but obtained from a modified Van der Pol electroinic oscillator are shown. Such an oscillator has very little inherent noise. Any noise that appears in the course of signal processing is due to ADC quantization. Fig. 5 shows a singular value spectrum 90 decreasing below an estimate of a quantization level of ADC 91 at approximately i = 13. Level 91 is the standard deviation of the relevant time series divided by one in the least significant bit channel of the associated ADC. The spectrum 90 becomes constant at approximately i = 20. The ADC quantization noise therefore dominates the results, and the spectrum requires truncation at i = 13 (giving d = 13). With a more sensitive 16 ADC, the truncation could extend to i = 20.

Fig. 5 also shows the shapes of the singular vectors or functions obtained in a singular value decomposition. There are thirty-eight, that is i = 1 to 38.

Referring now to FIG. 6, a functional block diagram of a portion of another embodiment 110 of the invention is shown. This is a dynamical systems analyzer which is a modification of the one shown in Fig. 1, and the parts equivalent to those described above have similar reference symbols. In the case of numerical references, 100 has been added to each one.

In view of the similarities between analyzers 10 and 110, the description of the latter would focus on the different aspects.

The analyzer 110 has a heuristic processor 144 and a FIR filter 134 equivalent to the

ES 2 149 172 T3 described above. The FIR filter 134 is illustrated laterally inverted so that its latch chain 136 appears on the right. The bolt chain 136 is extended by three additional bolts 136n + 1, 136n + 2 and 136n + 3 compared to the bolt chain 36. The bolts 136n to 136n + 3 are connected to the respective inputs SAI0 to SAI3 collectively forming an input. subsidiary response SAI. The inputs SAI1, SAI2 and SAI3 are connected via the respective input delay units 1151 to 1153 to a delay device 160 incorporated in the heurastic processor 144. There is no equivalent of switch 42. The input delay units 1151, 1152 and 1153 implement the respective delays of τ, 2τ, and 3τ, where τ is one clock cycle as before.

The signals pass from the delay device 160 into the respective columns C0 to C3, in combination forming an output matrix 158b. Matrix 158b is equivalent to a multi-column version of column 58b from the previous embodiment, in which a least squares fit to the time series data was implemented. Columns C0, C1, and C2 are in series with output delay units 1170, 1171, and 1172, respectively; the latter provides delays of 3Υ, 2τ and τ, respectively, and consequently signals reaching an output interface 168 are temporarily de-staggered (compared to the equivalent output of columns C0 to C3).

System analyzer 110 operates as follows. A training procedure equivalent to that described above is used. A comparison system (not shown) provides a series of signal times from which singular vectors f1 to fd are derived by singular value decomposition and truncation. The vectors are loaded into the FIR filter 134 and the time series V1 Vi ... is timed along the lock chain 136. In the training mode, the phase space coordinates exiting the FIR filter 134 are non-linearly transformed, QR decomposed, and fitted to the time series data as described for embodiment 10 above. In the present embodiment, instead of a single number (p) of clock cycle delays in device 46, there are four independent delays of 0, 1, 2 and 3 clock cycles (intervals τ). To implement this, the delay circuitry has been rearranged somewhat making use of the latch chain 136 and adjusting for delays experienced by the rotation parameters that pass between successive columns C0 to C3. The input delay units 1151 etc provide such an adjustment.

Consequently, during training the i-th set of phase space coordinates g1 (i) to gd (i) that come out of the FIR filter 134 is associated by the heurastic processor 144 to each of the following points of the series of times: Vi + n-1, Vi + n, Vi + n + 1 and Vi + n + 2. This is carried out by cutting into respective columns C0 to C3; that is, the Cq + 1 column is responsible for the adjustment of the transformed version from g1 (i) to gd (i) to Vi + n + q, where q is -1, 0, 1 or 2 and q is equivalent to p-1.

Columns C0 to C3 provide error values ​​between the transformed versions from g1 (i) to gd (i) and Vi + n-1, Vi + n, Vi + n + 10, and Vi + n + 2. The prediction interval increases from 0 to 3 clock cycles between columns C0 to C3, and therefore the non-conformities or errors are likely to increase as well.

Fig. 7 is an indication of what might be expected as the output of an analyzer similar to the analyzer 110 but having eleven columns of output C0 to C10 (not shown). It is a graph 200 of the logarithm in base 10 of the RMS prediction error with respect to the prediction interval, being this the number p of intervals (τ) of clock cycles; here p varies from 0 to 10, and is the sub-index of the respective output column C0, C1, C2, ..., C10 in each case. A graph such as Graph 200 can be used to identify the degree of predictability of a nonlinear dynamic system. The RMS prediction error shown in Fig. 7 is seen to increase as p increases to a maximum level at which it saturates. The saturation level is caused by the fact that even if the analyzer made the worst possible guess the maximum prediction error will be limited to some maximum, determined by the "attractor sizes" of training and processing.

The result shown in Fig. 7 is calculated as follows. The sequence of values ​​xt, t = 0 to 4000 are generated by iterating the equation in differences

Xt = 1 - a (xt-i)² + ext-2 (10)

The parameters α and β are constant with values ​​α = 1.4 and β = 0.3. These values ​​are chosen so that the sequence is chaotic, and a small initial value is used for x0. The sequence of values ​​of xt, where t = 1001 to 2000, is used to train the dynamic systems analyzer using 36 centers uniformly distributed through the plane bounded by the points (2, 2), (2, -2), ( -2, -2) and (-2, 2). The sequence of xt values, where t = 3001 to 4000 is then processed, and the log10 RMS prediction error is calculated using the respective set of 1000 error values ​​that appear from each of the output columns C0 to C10.

The invention can be implemented in an alternative way to that described above. Analyzer 10, for example, was operated in successive training and testing modes. However, it is possible for a training procedure on one analyzer to produce results used by a different device. Referring to Figures 1 through 1a once again, the training procedure in the heurastic processor 44 generates a set of array elements. The values ​​of matrix elements are terminated when the training is finished, and are stored in individual internal or boundary cells I and B. In test mode these values ​​are frozen, and implement a mathematical transformation (matrix) on the data of output of the digital arithmetic unit 52, which in turn implements a non-linear expansion around a set of centers. The

ES 2 149 172 T3 transformation matrix is ​​equivalent to multiplication by a vector of weights or matrix of weights; the individual elements of the weight vector or matrix can be obtained from the trained heuristic processor. The procedure is established in the PCT patent application No.^°PCT / GB90 / 00142 cited above, and will not be described in detail. However, having such singular weight elements and vectors determined in a training process, together with the centers and a non-linear transformation (unit 52) ​​used in training, the operation in the test mode can be by means of a simple device. of weighing that works at the output of the digital arithmetic unit 52. This allows the apparatus requirements of the test device to be reduced.

It has been found that the dynamical systems analyzer of the invention could, in certain circumstances, prove to be too sensitive to minor deviations of a dynamic system with respect to the expected behavior. In other words, a small deviation could lead to a big mistake. In this case it is possible to use the analyzer of the invention in a different way, as will now be described.

As said, the heuristic processor 44 operating in the training mode constitutes stored elements that are referred to a weight matrix. The elements of the weight matrix can be obtained at the output 68 by inputting the unit vectors to the triangular systolic matrix 54 as described in document PCT / GB90 / 00142. This allows a dynamical system to be characterized by the elements of the weight matrix that it gives rise to in the heuristic processor 44. The data of a dynamic comparison system can, therefore, be entered in the analyzer and used in it to generate the elements of the weight matrix for later reading. These elements are then available for comparison with the elements of the weight matrix of any test system. If the weight matrices of the comparison system and the test system are similar, then the two systems are similar.

In fact, it is not essential to convert the elements stored in the inner or boundary cells I and B into weight matrices to compare the test and comparison systems. Since each set of stored elements is directly related to a weight matrix, and the relationship is invariant, the stored elements of the test and comparison systems can be compared with each other directly without intervening the transformation to weight matrices.

There are two approaches to obtaining the weight matrix of the test system. One approach is to limit the analyzer 10 to using data from the test system; that is, after the analyzer 10 has been trained using the data from the comparison system, the weight matrix of the comparison system is read and the entire training process described above is carried out one more time. The computer 22 consequently performs a singular value decomposition of the data derived from the test system, and loads the FIR filter 34 with singular vectors derived from said decomposition. The heuristic processor 44 subsequently processes the output of the FIR filter 34 as described above to form elements stored in the inner or boundary cells.

The second approach is that the heuristic processor operates as described above for the training mode, but the FIR filter 34 is loaded with singular functions derived from the singular value decomposition of the comparison system data. There is no such decomposition of data from the test system. Consequently, the test system data is transformed into phase space coordinates using the singular functions of the comparison system. Therefore, the phase space coordinates thus produced undergo processing in the heuristic processor as described above to form the set of stored elements of the test system relative to the weight matrix.