Identification of model parameters in Engee

Example description

Let's consider the implementation of the task of parametric identification of a model based on data. We have to identify several time-invariant parameters of some model, the structure of which we know in advance.

In this example we will analyse the signal generated by the dynamic model Engee. It is located in the folder start/examples/data_analysis/parametric_identification.

The model consists of two subsystems: one generates the analysed signal, the other is a reference model (without noise).

In the subsystem ModelSystem, five signals are combined by means of summing:

• constant component: the number 5,
• time component: a line with a slope of 0.01,
• first harmonic: a sinusoid with amplitude 1 and angular frequency rad/s,
• less prominent second harmonic: sinusoid with amplitude 0.6 and angular frequency rad/s,
• white noise with zero offset and dispersion 0.3.

In the subsystem ModelBlock we will be able to substitute the identified parameters. As long as reference parameters are set in it and it generates a reference signal. In it the model is described as a code:

function (c::ModelBlock)(t::Real)
    T0 = 5
    T1 = 0.01
    A1, f1, phi1 = 1, 0.02, 0
    A2, f2, phi2 = 0.6, 0.008, 0
    
    signal = T0 + T1 * t +
             A1 * sin( 2*pi * f1 * t + phi1 ) .+
             A2 * sin( 2*pi * f2 * t + phi2 );
    return signal
end

Let's save the reference values of the parameters:

As a result of running the model, we get 500 seconds of signal. Here is a graph of the last 200 seconds:

The block ModelBlock is now initialised with exactly the same parameters as the dynamic model in the subsystem ModelSystem, so the graph of from_code (dark blue) completely replicates the systematic part of the signal from_model (orange).

What data the model stores

Both signals are saved in a CSV table param_id_data.csv with the following columns:

• time - time column,
• m - column of data generated by the dynamic model (reference parameters and error),
• c - column of data generated by the reference model (executed as a code).

The order of actions in this example

We have the following steps to perform:

1. start the dynamic model to synthesise the model data
2. find the parameters of all the model links in the correct order (assuming that the signal allows us to guess the topology of the model, or that we already know it).
3. check that after subtracting systematic components from the input signal, the result is normally distributed noise.
4. fit the parameters into the Engee Function block code and check that the models match each other.

Synthesis of model data

Run the following cell to synthesise the data we are going to analyse:

Loading libraries

Connection of libraries for data loading, working with tables, analysing spectra and library for polynomial regression.

Downloading and connecting the library for hypothesis testing (we use normality testing):

Let's load the data

From here on we will only work with the signal generated from the dynamic model, but you can see that both methods are interchangeable, and if necessary, this example can be accommodated in a single script.

We will visualise the model signal ModelSystem, which consists of a vector of time and calculated values:

We will also keep the sampling rate of the time signal (a value inverse to the distance in time between two neighbouring measurements):

What we will identify

The reference model ModelBlock is initialised with reference parameter values (the same as in the model ModelSystem).

The model is defined by 8 parameters:

• linear trend parameters:
• amplitude, frequency and phase of the first sinusoid:
• amplitude, frequency and phase of the second sinusoid:

We will not be interested in noise parameters, although they are easily identified by the functions mean and var.

Let's identify the linear component

The first thing we will do is to identify and eliminate the linear trend from the signal.

Let's keep the polynomial parameters as separate variables

Let's display the found trend on the graph and show how our data looks like minus this trend:

Great, the signal has started to look like a combination of sinusoids and noise.

Identify the first sinusoid

To determine the parameters of the sinusoids included in the signal, we plot the spectrum of the signal using the FFTW library.

Its discretisation will be proportional to the number of points in the original signal. That is, if we leave the original discretisation, we can get low accuracy of the spectrum. To get rid of the high-frequency region, we reduce the number of points in the original signal by lowering the sampling frequency of the original signal. To do this, we set the variable downsample_factor. The value 40 means that we will take every 40th point from the original signal.

Since we are only interested in the positive part of the spectrum, we use the function rfft (we also don't have to do fftshift):

Let's find the maximum on this spectrum and determine which frequency this value corresponds to.

When we have found the maximum point on the spectrum, we can identify the amplitude frequency and phase of the most prominent harmonic component. Now we need to do a recombination - generate a harmonic signal from a complex number that describes its parameters.

Let's define the first sinusoid as a function of time and plot it:

OK, now the signal is a combination of a sine wave and noise.

Let's identify the second sine wave

Repeat the same operations as for finding the first sinusoid:

Let's find the parameters of the second sinusoid

Let's define the second sinusoid as a function of time and plot it:

How do we know if there are significant systematic components left in this signal? One way is to imagine that this signal is a realisation of a random variable and test it for normality of distribution.

Analysing the residual signal

Let's test the hypothesis of normal distribution of the residual signal component using Harker-Bera test (JB-test):

Note the result (outcome) of the test:

• fail to reject h_0 means that there is no good reason to reject the hypothesis : with 95% confidence we assert that the observed signal is generated by a random process with a normal distribution
• reject h_0 means that the confidence that the observed signal originates from a normal distribution does not exceed 5%

Try to exclude some component from the formula - most likely such a signal will fail Harkke-Bera test and you will get the result reject h_0.

Analysing the results

The result obtained serves as a basis for completing the search for additional systematic components.

Let us derive the final report and the parameter estimation error.

Create a function that calculates the model data with identified parameters:

Let's compare the result of running the identified model and the result of running ModelBlock:

In this example we:

• implemented a fairly simple parametric identification procedure,
• determined the identification error of all parameters,
• verified that the graph of the model with identified parameters is quite similar to the graph of the reference model.