Volterra Black-Box Models Identification Methods: Direct Collocation vs. Least Squares

Abstract

The Volterra integral-functional series is the classic approach for nonlinear black box dynamical system modeling. It is widely employed in many domains including radiophysics, aerodynamics, electronic and electrical engineering and many others. Identifying the time-varying functional parameters, also known as Volterra kernels, poses a difficulty due to the curse of dimensionality. This refers to the exponential growth in the number of model parameters as the complexity of the input-output response increases. The least squares method (LSM) is widely acknowledged as the standard approach for tackling the issue of identifying parameters. Unfortunately, the LSM suffers with many drawbacks such as the sensitivity to outliers causing biased estimation, multicollinearity, overfitting and inefficiency with large datasets. This paper presents an alternative approach based on direct estimation of the Volterra kernels using the collocation method. Two model examples are studied. It is found that the collocation method presents a promising alternative for optimization, surpassing the traditional least squares method when it comes to the Volterra kernels identification including the case when input and output signals suffer from considerable measurement errors.

1. Introduction

At the current stage of development of wireless technologies like 5G/6G communication system networks based on antenna arrays with digital beam forming (Massive Multiple Input Multiple Output system), it is impossible to manage without such digital signal processing algorithms as digital correction of the nonlinear distortion DPD (Digital Predistortion). Nonlinear distortions of the signal occurring inside the transceiver path strongly distort the spectrum of this signal, as shown in Figure 1, where it is shown in red, and the main signal is blue in color, respectively.

However, the international wireless standards like 3GPP, ETSI impose strict requirements on the spectral power of the radiated signal. The use of digital nonlinear distortion correction algorithms allows for meeting the requirements of standards and at the same time positively affecting the overall efficiency, that is, the energy consumption of the entire signal receiving and transmitting system. There are different approaches to the implementation of such algorithms, both purely digital and analog and mixed. One of them, a purely mathematical approach to the description of nonlinear distortions, we will describe below. However, let us consider the general statement of the problem of digital correction (DPD) with the following structure of the model of correction as shown in Figure 2.

Here, $F_{DPD}(.)$ is a nonlinear operator reflecting the essence of nonlinear correction—imagine it as some function dependent on parameters $\overrightarrow{W}={[w_1,\dots ,w_p]}^T,\overrightarrow{W}\in {\mathbb{C}}^p$. $F_{PA}(.)$ is a nonlinear operator identified with a nonlinear device which generates some complex vector $\overrightarrow{Y}={[y_1,\dots ,y_n]}^T,\overrightarrow{Y}\in {\mathbb{C}}^n$ and also defines some vector from a complex field of numbers $\overrightarrow{Y_d}={[y_{d,1},\dots ,y_{d,n}]}^T,\overrightarrow{Y_d}\in {\mathbb{C}}^n$ on which the operator $F_{DPD}(.)$ depends. Under the error $E\in {\mathbb{C}}^n$ we will understand the difference between vectors Y and $Y_d$

$$\overrightarrow{E}=\overrightarrow{Y_d}-\overrightarrow{Y}.$$

Then we can formulate the requirements for the definition of parameters $\overrightarrow{W}$ as follows: $\overrightarrow{\omega }={arg\; min}_W{∥E∥}^2$, where $∥.∥$ is Euclidean norm. Considering $Y=F_{PA}\left(F_{DPD}\left(Y_d\right)\right),$ the above introduced expression can be rewritten as

$$\overrightarrow{\omega }=\underset{W}{arg\; min}{∥Y_d-F_{PA}\left(F_{DPD}\left(Y_d\right)\right)∥}^2.$$

This equation will be task of DPD (Digital Predistortion). Here, we can highlight several important sub-tasks, which in themselves are quite complex both theoretically and computationally:

(a)

Since we have formulated, in fact, the problem of approximation of a function, we need to derive the analytical regression dependence $F_{DPD}(.)$ on the parameters $\overrightarrow{W}$. How this function is defined will depend on the quality of the correction of nonlinear distortions;

(b)

The procedure of searching for the parameters $\overrightarrow{W}$ is a classical optimization problem, which is a linear or nonlinear regression with respect to the parameters $\overrightarrow{W}$. Finding efficient methods of convex or non-convex optimization is one of the major difficulties in this problem;

(c)

Compression of a function $F_{DPD}(.)$, i.e., reducing its computational complexity.

One of the methods to solve the problem (a) for the DPD task is the Volterra functional series. And it is also the conventional tool to characterize the complex nonlinear dynamics in various fields including the radiophysics, mechanical engineering, electronic and electrical engineering, energy sciences (here, readers may refer, e.g., to review [1] or [2]). Volterra series are widely employed to represent the input-output relationship of nonlinear dynamical systems with memory. Volterra power series are among the best-understood nonlinear system representations in signal processing. Such an integral functional series (also called Fréchet-Volterra series)

$$\begin{array}{c}y\left(t\right)=F\left(x\left(t\right)\right):=\underset{0}{\overset{t}{\int }}{K}_1\left(s\right)x(t-s)ds+\underset{0}{\overset{t}{\int }}\underset{0}{\overset{t}{\int }}{K}_2(s_1,s_2)x(t-s_1)x(t-s_2)d{s}_{1}d{s}_2+\dots \\ +\underset{0}{\overset{t}{\int }}\underset{0}{\overset{t}{\int }}\cdots \underset{0}{\overset{t}{\int }}{K}_n(s_1,s_2,\dots ,s_n)x(t-s_1)x(t-s_2)\dots x(t-s_n)d{s}_{1}d{s}_2\dots d{s}_n+\dots t\in [0,T]\end{array}$$

(1)

was proposed by Maurice Fréchet for a continuous nonlinear dynamical systems representation [3,4]. Here, readers may also refer to overview [5] and monograph [6] for more details on relevant Lyapunov–Liechtenstein operator and Lyapunov–Schmidt methods in the theory of non-linear equations.

The role of a reproducing kernel Hilbert space in the development of a unifying view of the Volterra theory and polynomial kernel regression is presented in [7].

In (1), $x\left(t\right)$ is the input signal and $y\left(t\right)$ is the output of a single-input-single-output (SISO) nonlinear system and $K_n(s_1,s_2,\dots ,s_n)$ are the multidimensional Volterra kernels (or transfer functions) to be identified based on nonlinear system’s response $y\left(t\right)$ as a reaction on input $x\left(t\right)$ (Figure 3). It is to be noted that for the basic case $n=1$, we have a conventional Finite Impulse Response (FIR) linear model which is optimal in the least-squares sense.

The Fréchet theorem [3] generalizes the famous Weierstrass approximation theorem which characterizes the set of continuous functions on a compact interval via uniform approximation by algebraic polynomials.

Power series (1) characterize the stationary dynamical systems. Stationarity here means that a transfer function does not vary during the transient process as $t\in [0,T]$. More general power series (2) models nonstationary dynamics when transfer functions depend explicitly on time t

$$\begin{array}{c}y\left(t\right)=\underset{0}{\overset{t}{\int }}{K}_1(t,s)x\left(s\right)ds+\underset{0}{\overset{t}{\int }}\underset{0}{\overset{t}{\int }}{K}_2(t,s_1,s_2)x\left(s_1\right)x\left(s_2\right)d{s}_{1}d{s}_2+\dots \\ +\underset{0}{\overset{t}{\int }}\underset{0}{\overset{t}{\int }}\cdots \underset{0}{\overset{t}{\int }}{K}_n(t,s_1,s_2,\dots ,s_n)x\left(s_1\right)x\left(s_2\right)\dots x\left(s_n\right)d{s}_{1}d{s}_2\dots d{s}_n+\dots t\in [0,T].\end{array}$$

(2)

The Volterra series is an essential tool for the mathematical modeling of the nonlinear dynamical systems appearing in the digital pre-distortion (DPD) iterative process [8]. DPD as we described before is an important part of the digital signal processing algorithms used in transmitters and receivers. Particularly for short-distance applications where the limitations of the transceiver are more significant. DPD is used to improve performance by compensating for the imperfect response of transmitter components, e.g., in [9] the frequency selective DPD was proposed. In [9] the Volterra series model structure consists of a basic linear part and a partial-band pre-compensation part, moreover, a generalized indirect learning architecture is employed to extract the coefficients. Several methods have been studied for DPD, with Volterra series-based methods being popular due to their ease of implementation and the straightforward interpretation of their nonlinear terms. The key issue with Volterra series is the curse of dimension: as the order of the series increases, the number of terms involved in the expansion grows exponentially, making it computationally demanding. On the other hand, estimating the functional coefficients (Volterra kernels) of the Volterra integral functional series can be challenging. It is often considered in its discrete form and requires a significant amount of data and complex optimization algorithms to find the best fit for the model coefficients. An alternative approach based on problem reduction to multi-dimensional integral equations solution [10,11] requires a special probe signal design.

In present paper, the alternative approach for the identification of Volterra kernels is proposed using the direct collocation method. The results are compared with the conventional least squares method (LSM) widely employed for the Volterra series identification problem in the telecommunication domain.

The rest of the paper is structured as follows: The subsequent section provides the problem statement. Section 3 focuses on the collocation method. Section 4 carries out computational experiments with LSM, while Section 5 discusses concluding remarks and future work.

2. Identification Problem Statement

Let us consider the following segment of the truncated Volterra series (1) for $n=2$

$$y\left(t\right)=\underset{0}{\overset{t}{\int }}{K}_1\left(s\right)x(t-s)ds+\underset{0}{\overset{t}{\int }}\underset{0}{\overset{t}{\int }}{K}_2(s_1,s_2)x(t-s_1)x(t-s_2)d{s}_{1}d{s}_2,t\in [0,T].$$

(3)

Our current problem in this section is to determine the kernels $K_1\left(s\right)$ and $K_2(s_1,s_2)$ by a known input and output pair $\left(x\left(t\right),y\left(t\right)\right)$.

In contrast to the linear case $n=1$, when it is sufficient to specify a single pair $\left(x\left(t\right),y\left(t\right)\right)$ to determine the kernel $K_1\left(s\right)$, in the nonlinear case $n=2$, for the unique identification of the two-dimensional kernel $K_2(s_1,s_2)$, it is necessary to specify a two-dimensional continuum of equalities. This means that problem (4) has an infinite set of solutions.

Remark 1.

It should be noted that if we consider this problem as an integral equation with two unknown functions $K_1\left(s\right)$ and $K_2(s_1,s_2)$, then this problem is essentially ill-posed. There are an infinite number of solutions and this problem is insufficiently defined. In this regard, no classical numerical methods designed for integral equations are applicable in this case. And as a result, there are no any attempts to solve the problem in this form in the literature.

Remark 2.

A fundamentally different situation takes place in the problem of determining an unknown input signal $x\left(t\right)$ with a known output signal $y\left(t\right)$ after kernels identification. It is to be noted that in this case we have the problem of nonlinear Volterra integral equations’ solution. Here, readers may refer to Section 9 in book [11], papers [12,13,14] and references therein regarding the Kantorovich principal solutions and the blow-up phenomenon.

Within the framework of this paper, from a practical point of view, we will be satisfied with any pair of approximately found kernels ${\tilde{K}}_1\left(s\right)$ and ${\tilde{K}}_2(s_1,s_2)$ that provides a sufficiently small residual norm

$$\epsilon =\underset{t\in [0,T]}{max}\left|y\left(t\right)-\underset{0}{\overset{t}{\int }}{\tilde{K}}_1\left(s\right)x(t-s)ds-\underset{0}{\overset{t}{\int }}\underset{0}{\overset{t}{\int }}{\tilde{K}}_2(s_1,s_2)x(t-s_1)x(t-s_2)d{s}_{1}d{s}_2\right|.$$

(4)

Denoted by $B_i\left(t\right),i=0,1,2,\dots ,$ the basis functions form a complete orthogonal system of functions on the segment $[0,T]$.

We look for an approximate solution of the problem (3) in the form of segments of series of expansions according to the selected system of basis functions

$${\tilde{K}}_{1,m}\left(s\right)=\sum _{i=0}^{m-1}{A}_i{B}_i\left(s\right),{\tilde{K}}_{2,m_1,m_2}(s_1,s_2)=\sum _{i=0}^{m_1-1}\sum _{j=0}^{m_2-1}{C}_{ij}{B}_i\left(s_1\right)B_j\left(s_2\right).$$

(5)

3. Collocation Method

Collocation-type methods are widely used in the discretization of various kinds of integro-functional equations [15]. With sufficiently good accuracy and stability, they are also computationally less expensive in comparison with projection methods of the Galerkin type requiring additional integration [16].

In order to determine the unknown coefficients $A_i$ and $C_{ij}$, we introduce a uniform grid of nodes

$$t_k\in [0,T],k=0,1,\dots ,N,$$

(6)

where $N+1$ is number of nodes.

Substitute (5) in (3) and then demand that the equalities be fulfilled at the points (6)

$$y\left(t_k\right)=\underset{0}{\overset{t_k}{\int }}{\tilde{K}}_{1,m}\left(s\right)x(t_k-s)ds+\underset{0}{\overset{t_k}{\int }}\underset{0}{\overset{t_k}{\int }}{\tilde{K}}_{2,m_1,m_2}(s_1,s_2)x(t_k-s_1)x(t_k-s_2)d{s}_{1}d{s}_2,k=\overline{0,N}.$$

(7)

Denote for a simplicity $y\left(t_k\right)=y_k$, and transform the last equalities as follows

$$y_k=\sum _{i=0}^{m-1}{A}_i\underset{0}{\overset{t_k}{\int }}{B}_i\left(s\right)x(t_k-s)ds+\sum _{i=0}^{m_1-1}\sum _{j=0}^{m_2-1}{C}_{ij}\underset{0}{\overset{t_k}{\int }}\underset{0}{\overset{t_k}{\int }}{B}_i\left(s_1\right)B_j\left(s_2\right)x(t_k-s_1)x(t_k-s_2)d{s}_{1}d{s}_2.$$

(8)

As a system of basis functions $B_i\left(t\right),i=0,1,\dots$, we choose Chebyshev polynomials of the first kind

$$T_0\left(t\right)=1,T_1\left(t\right)=t,T_{i+1}\left(t\right)=2t{T}_i\left(t\right)-T_{i-1}\left(t\right),i=1,2,\dots .$$

(9)

Sufficient conditions for the applicability of Chebyshev polynomial expansions of the form (5) are the limitation of the first derivatives of the approximated kernels. For more detailed information about convergence, we refer, for example, to the book [17].

Since these polynomials are orthogonal on the segment $[-1,1]$, we apply a linear mapping to the segment $[0,T]$.

The controlled norm of the residual corresponding to the selected values of $m,m_1$ and $m_2$ takes the form

$$\begin{array}{c}\hfill {\epsilon }_N=\underset{t\in [0,T]}{max}\left|y\left(t\right)-\sum _{i=0}^{m-1}{A}_i\underset{0}{\overset{t}{\int }}{B}_i\left(s\right)x(t-s)ds-\right.\\ \hfill \left.\sum _{i=0}^{m_1-1}\sum _{j=0}^{m_2-1}{C}_{ij}\underset{0}{\overset{t}{\int }}\underset{0}{\overset{t}{\int }}{B}_i\left(s_1\right)B_j\left(s_2\right)x(t-s_1)x(t-s_2)d{s}_{1}d{s}_2\right|.\end{array}$$

(10)

Let us denote $N=m+m_1{m}_2-1$. The number of equalities (number of nodes in the grid) equals the number of unknown coefficients.

Thus, we have the following system of linear algebraic equations

$$y_k=\sum _{i=0}^{m-1}{A}_i{\beta }_{ik}+\sum _{i=0}^{m_1-1}\sum _{j=0}^{m_2-1}{C}_{ij}{\gamma }_{ijk},k=\overline{0,m+m_1{m}_2-1},$$

(11)

with respect to the unknown coefficients $A_i,i=0,1,\dots ,m-1$ and $C_{ij},i=0,1,\dots ,m_1-1,j=0,1,\dots ,m_2-1$. Here,

$${\beta }_{ik}=\underset{0}{\overset{t_k}{\int }}{B}_i\left(s\right)x(t_k-s)ds,{\gamma }_{ijk}=\underset{0}{\overset{t_k}{\int }}\underset{0}{\overset{t_k}{\int }}{B}_i\left(s_1\right)B_j\left(s_2\right)x(t_k-s_1)x(t_k-s_2)d{s}_{1}d{s}_2.$$

(12)

4. Least–Square Method

Let us denote $N>m+m_1{m}_2-1$. We have the situation where number of equalities is larger than number of unknown coefficients $A_i$ and $C_{ij}$. Thus we have the overdetermined system of linear equations with respect to the unknown coefficients $A_i,i=0,1,\dots ,m-1$ and $C_{ij},i=0,1,\dots ,m_1-1,j=0,1,\dots ,m_2-1$:

$$y_k=\sum _{i=0}^{m-1}{A}_i{\beta }_{ik}+\sum _{i=0}^{m_1-1}\sum _{j=0}^{m_2-1}{C}_{ij}{\gamma }_{ijk},$$

(13)

where

$${\beta }_{ik}=\underset{0}{\overset{t_k}{\int }}{T}_i\left(s\right)x(t_k-s)ds,{\gamma }_{ijk}=\underset{0}{\overset{t_k}{\int }}\underset{0}{\overset{t_k}{\int }}{T}_i\left(s_1\right)T_j\left(s_2\right)x(t_k-s_1)x(t_k-s_2)d{s}_{1}d{s}_2.$$

(14)

The system is inconsistent. The least–square method is used to find the approximate solution of the system. The point of the method is to find such coefficients $A_i$ and $C_{ij}$ such that the following criteria is minimized:

$$\sum _{k=0}^{N-1}{\left(y_k-\sum _{i=0}^{m-1}{A}_i{\beta }_{ik}-\sum _{i=0}^{m_1-1}\sum _{j=0}^{m_2-1}{C}_{ij}{\gamma }_{ijk}\right)}^2⟶min.$$

(15)

5. Numerical Experiments

Let us illustrate the operation of the proposed identification methods on two pairs of model signals.

5.1. Model 1. Periodic Signal

Let us consider the case of periodic input signal, where

$$\begin{array}{c}\hfill x\left(t\right)=sin\left(20t\right),y\left(t\right)=\frac{1}{81002}(199{cos}^2\left(20t\right)-15sin\left(40t\right)-200cos\left(20t\right)e^{-2t}+1+\\ \hfill 10sin\left(20t\right)e^{-2t}+20sin\left(20t\right)e^{-t})+\frac{1}{409}\left(3sin\left(20t\right)-20cos\left(20t\right)+\frac{850920}{40501}{e}^{-3t}\right).\end{array}$$

(16)

The Figure 4 shows the graphs of the input $x\left(t\right)$ and output signal $y\left(t\right).$

5.1.1. Collocation Method Results for the Model 1

Table 1. Dependence of the residual ${\epsilon }_N$ on the values $m,m_1,m_2$ for Model 1.

m	${\mathit{\epsilon }}_{\mathit{N}}$
3	$1.41\times {10}^{-2}$
4	$1.14\times {10}^{-6}$
5	$4.72\times {10}^{-9}$
6	$1.77\times {10}^{-12}$
7	$1.83\times {10}^{-14}$
8	$1.53\times {10}^{-18}$
10	$2.84\times {10}^{-26}$

demonstrates the dependence of the residual ${\epsilon }_N$ on the values $m=m_1=m_2$ for the uniform mesh $t_k=\frac{k}{N},k=0,1,\dots ,N,$ covering the segment $[0,1]$.

All calculations were performed in the Maple system with parameter Digits:=30 (the number of digits that Maple uses when making calculations with software floating-point numbers). It should be noted that when using other values of parameter Digits, the order of the residual changes on average in direct proportion to this value. Also note that the integration during the formation of the system (11) was carried out analytically and did not introduce additional error in the calculation results. This is due to the fact that the input signal $x\left(t\right)$ in most cases allows for the analytical calculation of the values (12). In the case of using input signals of a more complex structure, special approximation methods should be applied to the integrals (12), taking into account the possible fast oscillation of $x\left(t\right)$. Figure 5, Figure 6, Figure 7 and Figure 8 demonstrate residual error for fixed parameter m.

5.1.2. Least–Square Method Results for the Model (16)

For simplicity, we assume that $m=m_1=m_2$.

Table 2. Dependence of the residual ${\epsilon }_N$ on the values m and k.

	m = 3	m = 5	m = 7
$k=(m+m^2)\times 2$	$8.07\times {10}^{-4}$	$4.92\times {10}^{-10}$	$2.50\times {10}^{-16}$
$k=(m+m^2)\times 5$	$8.07\times {10}^{-4}$	$3.90\times {10}^{-10}$	$1.50\times {10}^{-16}$
$k=(m+m^2)\times 10$	$8.07\times {10}^{-4}$	$4.90\times {10}^{-10}$	$2.87\times {10}^{-15}$

demonstrates the dependence of the residual ${\epsilon }_N$ on the parameters.

All calculations for least–square method were performed in MATLAB. Overdetermined matrix is solved using the lsqminnorm function. lsqminnorm solves the linear equation $AX=B$ and minimizes $\parallel AX-B\parallel$. The function uses the complete orthogonal decomposition to find a low-rank approximation of a matrix. It also should be noted that all the integrations during calculation were carried out analytically and did not introduce additional error in the results. Residual error for fixed $m=5$ and $k=(m+m^2)\times 5$ is presented in Figure 9. Residual error for fixed $m=7$ and $k=(m+m^2)\times 5$ is presented in Figure 10.

5.2. Model 2. Fading Input Signal

Let us consider the case of fading input signal, where

$$\begin{array}{c}\hfill x\left(t\right)=e^{-3t}sin\left(10t\right),\\ \hfill y\left(t\right)=\underset{0}{\overset{t}{\int }}cos\left(\frac{s}{2}\right)x(t-s)ds+\\ \hfill \underset{0}{\overset{t}{\int }}\underset{0}{\overset{t}{\int }}sin(s_1+2s_2)x(t-s_1)x(t-s_2)d{s}_{1}d{s}_2.\end{array}$$

(17)

The Figure 11 shows the graphs of the input signal $x\left(t\right)$ and output signal $y\left(t\right)$.

5.2.1. Collocation Method Results for Model 2

Table 3. Dependence of the residual ${\epsilon }_N$ on the values $m,m_1,m_2$.

m	${\mathit{\epsilon }}_{\mathit{N}}$
3	$3.16\times {10}^{-5}$
4	$9.85\times {10}^{-9}$
5	$8.58\times {10}^{-12}$
6	$2.17\times {10}^{-16}$
7	$5.37\times {10}^{-20}$

demonstrates the dependence of the residual ${\epsilon }_N$ on the values $m=m_1=m_2$ for the uniform mesh $t_k=\frac{k}{N},k=0,1,\dots ,N,$ covering the segment $[0,1]$.

Let us also discuss the stability of suggested numerical technique. Let the input data of the problem (17) be determined with some random error ${\epsilon }_{rand}$ varying within the $\delta$ value, namely $\left|{\epsilon }_{rand}\right|⩽\delta$.

Table 4. Stability results for collocation.

$\mathit{\delta }$	${\mathit{\epsilon }}_{\mathit{N}}$
${10}^{-2}$	$0.01729$
${10}^{-3}$	$2.71\times {10}^{-3}$
${10}^{-4}$	$2.56\times {10}^{-4}$
${10}^{-5}$	$7.54\times {10}^{-5}$
${10}^{-6}$	$1.66\times {10}^{-5}$

shows the dependence of the averaged residual ${\epsilon }_N$ on the $\delta$ value at a fixed $m=3$ based on the results of 10 measurements. Figure 12 and Figure 13 demonstrate the dependence of residual error with fixed parameter m.

It can be seen from the results of the Table 4 that residual continuously depends on the limits of random measurement errors of the input and output signals. Thus, we can conclude about the stability of the suggested method.

5.2.2. Least–Square Method Results for Model 2

Table 5. Dependence of the residual ${\epsilon }_N$ on the values m and k.

	m = 3	m = 5	m = 7
$k=(m+m^2)\times 2$	$2.38\times {10}^{-6}$	$7.77\times {10}^{-14}$	$2.93\times {10}^{-16}$
$k=(m+m^2)\times 5$	$2.63\times {10}^{-6}$	$7.46\times {10}^{-14}$	$3.05\times {10}^{-16}$
$k=(m+m^2)\times 10$	$3.48\times {10}^{-6}$	$7.41\times {10}^{-14}$	$3.80\times {10}^{-16}$

demonstrates the dependence of the residual ${\epsilon }_N$ on the parameters. Figure 14 demonstrates the residual error for fixed m and k.

As for collocation method, let us check the stability of the least-square method on this model. For testing stability, 10 rounds of experiments were performed and the average residual ${\epsilon }_N$ was calculated. Also, $m=3,k=(m+m^2)\times 5$ were fixed.

Table 6. Stability results for LSM.

$\mathit{\delta }$	${\mathit{\epsilon }}_{\mathit{N}}$
${10}^{-2}$	$0.00628$
${10}^{-3}$	$5.11\times {10}^{-4}$
${10}^{-4}$	$6.02\times {10}^{-5}$
${10}^{-5}$	$5.36\times {10}^{-6}$
${10}^{-6}$	$2.64\times {10}^{-6}$

shows the stability results.

6. Conclusions

Two numerical approaches to solving the problem of identification of the Volterra model were proposed in the paper. As can be seen from the presented results, both methods showed stable convergence. Convergence here can be interpreted only as the dependence of the residual on the increase in the number of terms in the expansions of kernels by Chebyshev polynomials (5). And this dependence is presented in numerical results.

It is to be noted, from the point of view of the arithmetic complexity of calculations, the collocation method turns out to be less expensive. And this factor is more pronounced the more parameters of the model are to be determined. This is due to the need to calculate a significantly larger number of integrals proportional to the square of the number of measurements being processed.

Further development of research suggests an increase in the number of terms n in the model (1) to identify a more accurate functional relationship between the input and output signals. It is also planned to develop special methods for approximating integrals (12) for the case of using input signals of a more complex structure, including fast oscillating signals.