# Volterra black-box models identification methods: direct collocation vs
least squares

Denis Sidorov
Applied Mathematics Department
Melentiev Energy Systems Institute SB RAS
Irkutsk, Russia
<EMAIL_ADDRESS>
&Aleksandr Tynda
Higher and Applied Mathematics Department
Penza State University
Penza, Russia
<EMAIL_ADDRESS>
&Vladislav Muratov
Institute of Mathematics and Information Technologies
Irkutsk State University
Irkutsk, Russia
<EMAIL_ADDRESS>
&Eugeny Yanitsky
Intermediate Radio Frequency Lab
Huawei Russian Research Institute
Moscow, Russia
<EMAIL_ADDRESS>

###### Abstract

The Volterra integral-functional series is the classic approach for nonlinear
black box dynamical systems modeling. It is widely employed in many domains
including radiophysics, aerodynamics, electronic and electrical engineering
and many other. 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 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.

_K_ eywords Volterra series $\cdot$ collocation method $\cdot$ kernels
identification $\cdot$ Chebyshev polynomials $\cdot$ memory effects

##  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
do 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 Fig. 1, where its shown in blue, and the main
signal is red color respectively.

Figure 1: Power spectrum density [1]

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 to meet the requirements of
standards and at the same time positively affect 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’s consider the general statement of the problem
of digital correction (DPD) with the following structure of the some model of
correction as shown in Fig. 2.

Figure 2: Digital correction scheme

Here $F_{DPD}(.)$ is a nonlinear operator reflecting the essence of nonlinear
correction and imagine it as some function dependent on parameters
$\vec{W}=[w_{1},...,w_{p}]^{T},\vec{W}\in\mathbf{C}^{p}$. $F_{PA}(.)$ is a
nonlinear operator identified with a nonlinear device and generating some
complex vector $\vec{Y}=[y_{1},...,y_{n}]^{T},\vec{Y}\in\mathbf{C}^{n}$ and
also define some vector from a complex field of numbers
$\vec{Y_{d}}=[y_{d,1},...,y_{d,n}]^{T},\vec{Y_{d}}\in\mathbf{C}^{n}$ on which
the operator $F_{DPD}(.)$ depends. Under the error $E\in\mathbf{C}^{n}$ we
will understand the difference between vectors $Y$ and $Y_{d}$

$\vec{E}=\vec{Y_{d}}-\vec{Y}.$

Then we can formulate the requirements for the definition of parameters
$\vec{W}$ as follows:
$\vec{\omega}=\operatorname*{arg\,min}_{W}\left\|E\right\|^{2}$, where
$\left\|.\right\|$ is Euclidean norm. Considering $Y=F_{PA}(F_{DPD}(Y_{d})),$
the above introduced expression can be rewritten as

$\vec{\omega}=\operatorname*{arg\,min}_{W}\left\|Y_{d}-F_{PA}(F_{DPD}(Y_{d}))\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 $\vec{W}$. How this function is defined will depend on the
quality of the correction of nonlinear distortions;

b) A procedure of the searching the parameters $\vec{W}$ is a classical
optimization problem, which is a linear or nonlinear regression with respect
to the parameters $\vec{W}$. Finding efficient methods of convex or non-convex
optimization is one of the major problems 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 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 [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 integral
functional series (also called Fréchet-Volterra series) (1)

$y(t)=F(x(t)):=\int\limits_{0}^{t}K_{1}(s)x(t-s)\,ds+\int\limits_{0}^{t}\int\limits_{0}^{t}K_{2}(s_{1},s_{2})x(t-s_{1})x(t-s_{2})\,ds_{1}ds_{2}+\dots$
(1)
$\dots+\int\limits_{0}^{t}\int\limits_{0}^{t}\dots\int\limits_{0}^{t}K_{n}(s_{1},s_{2},\dots,s_{n})x(t-s_{1})x(t-s_{2})\dots
x(t-s_{n})\,ds_{1}ds_{2}\dots ds_{n}+\dots\;t\in[0,T]$

were first 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 development of a unifying
view of the Volterra theory and polynomial kernel regression is presented in
[7].

In (1) $x(t)$ is input signal and $y(t)$ is 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(t)$ as a reaction on input $x(t)$
(Fig. 3). It is to be noted, that for the basic case $n=1$ we have a
conventional Finite Impulse Response (FIR) linear model which optimal in the
least-squares sense.

$x(t)$$F(x(t))$$y(t)$Figure 3. Behavioral modeling of the black box system

Fréchet theorem [3] generalises 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 functions do 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$

$y(t)=\int\limits_{0}^{t}K_{1}(t,s)x(s)\,ds+\int\limits_{0}^{t}\int\limits_{0}^{t}K_{2}(t,s_{1},s_{2})x(s_{1})x(s_{2})\,ds_{1}ds_{2}+\dots$
(2)
$\dots+\int\limits_{0}^{t}\int\limits_{0}^{t}\dots\int\limits_{0}^{t}K_{n}(t,s_{1},s_{2},\dots,s_{n})x(s_{1})x(s_{2})\dots
x(s_{n})\,ds_{1}ds_{2}\dots ds_{n}+\dots\;t\in[0,T].$

The Volterra series is essential tool of mathematical modeling the nonlinear
dynamical systems appearing in 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. 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. From other
hand, estimating the functional coefficients (Volterra kernels) of the
Volterra integral functional series can be challenging. It 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.
Alternative approach based on problem reduction to multi-dimensional integral
equations solution [10, 11] needs special probe signals design.

In present paper the alternative approach for the Volterra kernels
identification is proposed using the direct collocation method. The results
are compared with the conventional least squares method (LSM) widely employed
for Volterra series identification problem in 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(t)=\int\limits_{0}^{t}K_{1}(s)x(t-s)\,ds+\int\limits_{0}^{t}\int\limits_{0}^{t}K_{2}(s_{1},s_{2})x(t-s_{1})x(t-s_{2})\,ds_{1}ds_{2},\;t\in[0,T].$
(3)

Our current problem in this section is to determine the kernels $K_{1}(s)$ and
$K_{2}(s_{1},s_{2})$ by a known input and output pair
$\bigl{(}x(t),y(t)\bigr{)}$.

In contrast to the linear case $n=1$, when it is sufficient to specify a
single pair $\bigl{(}x(t),y(t)\bigr{)}$ to determine the kernel $K_{1}(s)$, 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 A1

It should be noted that if we consider this problem as an integral equation
with two unknown functions $K_{1}(s)$ and $K_{2}(s_{1},s_{2})$, then this
problem is essentially ill-posed. There are an infinite number of solutions,
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 A2

A fundamentally different situation takes place in the problem of determining
an unknown input signal $x(t)$ with a known output signal $y(t)$ 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
sec. 9 in book [11], papers [12], [13] and references therein regarding the
Kantorovich main solutions and 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 $\widetilde{K}_{1}(s)$
and $\widetilde{K}_{2}(s_{1},s_{2})$ that provides a sufficiently small
residual norm

$\varepsilon=\max\limits_{t\in[0,T]}\left|y(t)-\int\limits_{0}^{t}\widetilde{K}_{1}(s)x(t-s)\,ds-\int\limits_{0}^{t}\int\limits_{0}^{t}\widetilde{K}_{2}(s_{1},s_{2})x(t-s_{1})x(t-s_{2})\,ds_{1}ds_{2}\right|.$
(4)

Denote by $B_{i}(t),\;i=0,1,2,\ldots,$ the basis functions forming 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

$\widetilde{K}_{1,m}(s)=\sum\limits_{i=0}^{m-1}A_{i}B_{i}(s),\quad\widetilde{K}_{2,m_{1},m_{2}}(s_{1},s_{2})=\sum\limits_{i=0}^{m_{1}-1}\sum\limits_{j=0}^{m_{2}-1}C_{ij}B_{i}(s_{1})B_{j}(s_{2}).$
(5)

##  3 Collocation method

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

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,\ldots,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(t_{k})=\int\limits_{0}^{t_{k}}\widetilde{K}_{1,m}(s)x(t_{k}-s)\,ds+\int\limits_{0}^{t_{k}}\int\limits_{0}^{t_{k}}\widetilde{K}_{2,m_{1},m_{2}}(s_{1},s_{2})x(t_{k}-s_{1})x(t_{k}-s_{2})\,ds_{1}ds_{2},\;k=\overline{0,N}.$
(7)

Denote for a simplicity $y(t_{k})=y_{k}$, and transform the last equalities as
follows

$y_{k}=\sum\limits_{i=0}^{m-1}A_{i}\int\limits_{0}^{t_{k}}B_{i}(s)x(t_{k}-s)\,ds+\sum\limits_{i=0}^{m_{1}-1}\sum\limits_{j=0}^{m_{2}-1}C_{ij}\int\limits_{0}^{t_{k}}\int\limits_{0}^{t_{k}}B_{i}(s_{1})B_{j}(s_{2})x(t_{k}-s_{1})x(t_{k}-s_{2})\,ds_{1}ds_{2}.$
(8)

As a system of basis functions $B_{i}(t),\;i=0,1,\ldots$, we choose Chebyshev
polynomials of the first kind

$T_{0}(t)=1,\;T_{1}(t)=t,\;T_{i+1}(t)=2tT_{i}(t)-T_{i-1}(t),\;i=1,2,\ldots.$
(9)

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

The controlled norm of residual corresponding to the selected values of
$m,m_{1}$ and $m_{2}$ takes the form

$\varepsilon_{N}=\max\limits_{t\in[0,T]}\left|y(t)-\sum\limits_{i=0}^{m-1}A_{i}\int\limits_{0}^{t}B_{i}(s)x(t-s)\,ds-\right.\\\
\left.-\sum\limits_{i=0}^{m_{1}-1}\sum\limits_{j=0}^{m_{2}-1}C_{ij}\int\limits_{0}^{t}\int\limits_{0}^{t}B_{i}(s_{1})B_{j}(s_{2})x(t-s_{1})x(t-s_{2})\,ds_{1}ds_{2}\right|$
(10)

Let us denote $N=m+m_{1}m_{2}-1$. Number of equalities (number of nodes in the
grid) equals to the number of unknown coefficients.

Thus, we have the following system of linear algebraic equations

$y_{k}=\sum\limits_{i=0}^{m-1}A_{i}\beta_{ik}+\sum\limits_{i=0}^{m_{1}-1}\sum\limits_{j=0}^{m_{2}-1}C_{ij}\gamma_{ijk},\quad
k=\overline{0,m+m_{1}m_{2}-1},$ (11)

with respect to the unknown coefficients $A_{i},\;i=0,1,\ldots,m-1$ and
$C_{ij},\;i=0,1,\ldots,m_{1}-1,\;j=0,1,\ldots,m_{2}-1$. Here

$\beta_{ik}=\int\limits_{0}^{t_{k}}B_{i}(s)\,x(t_{k}-s)\,ds,\quad\gamma_{ijk}=\int\limits_{0}^{t_{k}}\int\limits_{0}^{t_{k}}B_{i}(s_{1})B_{j}(s_{2})x(t_{k}-s_{1})x(t_{k}-s_{2})\,ds_{1}ds_{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,\ldots,m-1$ and
$C_{ij},\;i=0,1,\ldots,m_{1}-1,\;j=0,1,\ldots,m_{2}-1$:

$y_{k}=\sum\limits_{i=0}^{m-1}A_{i}\beta_{ik}+\sum\limits_{i=0}^{m_{1}-1}\sum\limits_{j=0}^{m_{2}-1}C_{ij}\gamma_{ijk},$
(13)

where

$\beta_{ik}=\int\limits_{0}^{t_{k}}T_{i}(s)x(t_{k}-s)\,ds,\quad\gamma_{ijk}=\int\limits_{0}^{t_{k}}\int\limits_{0}^{t_{k}}T_{i}(s_{1})T_{j}(s_{2})x(t_{k}-s_{1})x(t_{k}-s_{2})\,ds_{1}ds_{2}.$
(14)

The system is inconsistent. 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\limits_{k=0}^{N-1}\left(y_{k}-\sum\limits_{i=0}^{m-1}A_{i}\beta_{ik}-\sum\limits_{i=0}^{m_{1}-1}\sum\limits_{j=0}^{m_{2}-1}C_{ij}\gamma_{ijk}\right)^{2}\longrightarrow\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

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

The Figure 4 shows the graphs of the input and output signal (16).

Figure 4: Input and output functions (16)

####  5.1.1 Collocation method results for the model (16)

The Table 1 demonstrates the dependence of the residual $\varepsilon_{N}$ on
the values $m=m_{1}=m_{2}$ for the uniform mesh
$t_{k}=\frac{k}{N},\;k=0,1,\ldots,N,$ covering the segment $[0,1]$.

Table 1: Dependence of the residual $\varepsilon_{N}$ on the values $m,m_{1},m_{2}$. $\mathbf{m}$ |  $\varepsilon_{N}$
---|---
3 |  $1,41\cdot 10^{-2}$
4 |  $1,14\cdot 10^{-6}$
5 |  $4,72\cdot 10^{-9}$
6 |  $1,77\cdot 10^{-12}$
7 |  $1,83\cdot 10^{-14}$
8 |  $1,53\cdot 10^{-18}$
10 |  $2,84\cdot 10^{-26}$

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). 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(t)$ in most cases allows 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(t)$.

Figure 5: Residual for $m=3$

Figure 6: Residual for $m=5$

Figure 7: Residual for $m=7$

Figure 8: Residual for $m=10$

####  5.1.2 Least–square method results for the model (16)

For the simplicity we assume that $m=m_{1}=m_{2}$. The Table 2 demonstrates
the dependence of the residual $\varepsilon_{N}$ on the parameters.

Table 2: Dependence of the residual $\varepsilon_{N}$ on the values $m$ and $k$. |  m = 3 |  m = 5 |  m = 7
---|---|---|---
$\mathbf{k=(m+m^{2})\cdot 2}$ |  $8.07\cdot 10^{-4}$ |  $4.92\cdot 10^{-10}$ |  $2.50\cdot 10^{-16}$
$\mathbf{k=(m+m^{2})\cdot 5}$ |  $8.07\cdot 10^{-4}$ |  $3.90\cdot 10^{-10}$ |  $1.50\cdot 10^{-16}$
$\mathbf{k=(m+m^{2})\cdot 10}$ |  $8.07\cdot 10^{-4}$ |  $4.90\cdot 10^{-10}$ |  $2.87\cdot 10^{-15}$

All calculations for least–square method were performed in MATLAB.
Overdetermined matrix is solved using lsqminnorm function. It also should be
noted that all the integrations during calculation were carried out
analytically and didn’t introduce additional error in the results.

Figure 9: Residual for $m=5$ and $k=(m+m^{2})\cdot 5$ Figure 10: Residual for
$m=7$ and $k=(m+m^{2})\cdot 5$

###  5.2 Model 2. Fading input signal

$\begin{split}x(t)=e^{-3t}\sin(10t),\\\
y(t)=\int\limits_{0}^{t}\cos\left(\frac{s}{2}\right)x(t-s)\,ds+\\\
+\int\limits_{0}^{t}\int\limits_{0}^{t}\sin(s_{1}+2s_{2})x(t-s_{1})x(t-s_{2})\,ds_{1}ds_{2}.\end{split}$
(17)

The Figure 11 shows the graphs of the input and output signals (17).

Figure 11: Input and output functions (17)

####  5.2.1 Collocation method results for the model (17)

The Table 3 demonstrates the dependence of the residual $\varepsilon_{N}$ on
the values $m=m_{1}=m_{2}$ for the uniform mesh
$t_{k}=\frac{k}{N},\;k=0,1,\ldots,N,$ covering the segment $[0,1]$.

Table 3: Dependence of the residual $\varepsilon_{N}$ on the values $m,m_{1},m_{2}$. $\mathbf{m}$ |  $\varepsilon_{N}$
---|---
3 |  $3,16\cdot 10^{-5}$
4 |  $9,85\cdot 10^{-9}$
5 |  $8,58\cdot 10^{-12}$
6 |  $2,17\cdot 10^{-16}$
7 |  $5,37\cdot 10^{-20}$

Figure 12: Residual for $m=3$

Figure 13: Residual for $m=7$

Let us also discuss the stability of suggested numerical technique. Let the
input data of the problem (17) be determined with some random error
$\varepsilon_{rand}$ varying within the $\delta$ value, namely
$\left|\varepsilon_{rand}\right|\leqslant\delta$. Table 4 shows the dependence
of the averaged residual $\varepsilon_{N}$ on the $\delta$ value at a fixed
$m=3$ based on the results of $10$ measurements.

Table 4: Stability results for collocation $\delta$ |  $\varepsilon_{N}$
---|---
$10^{-2}$ |  $0.01729$
$10^{-3}$ |  $2.71\cdot 10^{-3}$
$10^{-4}$ |  $2.56\cdot 10^{-4}$
$10^{-5}$ |  $7,54\cdot 10^{-5}$
$10^{-6}$ |  $1,66\cdot 10^{-5}$

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 the model (17)

Table 5: Dependence of the residual $\varepsilon_{N}$ on the values $m$ and $k$. |  m = 3 |  m = 5 |  m = 7
---|---|---|---
$\mathbf{k=(m+m^{2})\cdot 2}$ |  $2.38\cdot 10^{-6}$ |  $7.77\cdot 10^{-14}$ |  $2.93\cdot 10^{-16}$
$\mathbf{k=(m+m^{2})\cdot 5}$ |  $2.63\cdot 10^{-6}$ |  $7.46\cdot 10^{-14}$ |  $3.05\cdot 10^{-16}$
$\mathbf{k=(m+m^{2})\cdot 10}$ |  $3.48\cdot 10^{-6}$ |  $7.41\cdot 10^{-14}$ |  $3.80\cdot 10^{-16}$
Figure 14: Residual for $m=7$ and $k=(m+m^{2})\cdot 7$

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

Table 6: Stability results for LSM $\delta$ |  $\varepsilon_{N}$
---|---
$10^{-2}$ |  $0.00628$
$10^{-3}$ |  $5.11\cdot 10^{-4}$
$10^{-4}$ |  $6.02\cdot 10^{-5}$
$10^{-5}$ |  $5.36\cdot 10^{-6}$
$10^{-6}$ |  $2.64\cdot 10^{-6}$

##  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 (in the sense of the tendency
of the residual to zero). However, 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.

## References

  * [1] Guan, L., Kearney, R., Yu, C., and Zhu, A. High-performance digital predistortion test platform development or wideband RF power amplifiers. International Journal of Microwave and Wireless Technologies, 5(02), 2013, pp. 149–162.
  * [2] Cheng, C.M., Peng, Z.K., Zhang, W.M. and Meng, G. Volterra-series-based nonlinear system modeling and its engineering applications: A state-of-the-art review. Mechanical Systems and Signal Processing, Volume 87, Part A, 2017, pp. 340–364.
  * [3] Fréchet, M. Sur les fonctionnels continues. Ann. Ec. Norm. Sup. 1910\.
  * [4] Volterra, V. Theory of functionals and of integral and integrodifferential equations. Dover Publications: New York, 2005.
  * [5] Vainberg, M.M., Trenogin, V.A. The methods of Lyapunov and Schmidt in the theory of non-linear equations and their further development. Russian Math. Surveys, 17:2, 1962, pp. 1-–60.
  * [6] Sidorov, N.A., Loginov B.V., Sinitsin A.V. and Falaleev M.V. Lyapunov-Schmidt methods in nonlinear analysis and applications. Dordrecht: Springer Science+Business Media, 2013.
  * [7] Franz, M.O., Schölkopf, B. A Unifying View of Wiener and Volterra Theory and Polynomial Kernel Regression. Neural Computation, vol. 18, no. 12, 2006, pp. 3097-3118.
  * [8] Schreurs, D., O’Droma, M., Goacher, A.A., Gadringer, M. RF power amplifier amplifier behavioral modeling. The Cambridge RF and Microwave Engineering Series. Cambridge Univ. Press, 2009.
  * [9] Xia, X., Quan, X., Liu, Y., Shao, S., Tang, Y. A Frequency-Selective Digital Predistortion Method Based on a Generalized Indirect Learning Architecture. IEEE Transactions on Signal Processing, vol. 70, 2022, pp. 2334-2348.
  * [10] Apartsyn, A. Nonclassical linear Volterra equations of the first kind. Walter de Gruyter: Utrecht, 2003.
  * [11] Sidorov, D. Integral Dynamical Models: Singularities, Signals & Control, World Scientific: S’pore, 2014.
  * [12] Apartsin, A.S. Studying the polynomial Volterra equation of the first kind for solution stability. Automation and Remote Control. 72, 2011, pp. 1229–1236.
  * [13] Sidorov, D.N., Sidorov, N.A. Convex majorants method in the theory of nonlinear Volterra equations. Banach J. Math. Anal. 6(1), 2012, pp. 1–10.
  * [14] Brunner, H., Collocation Methods for Volterra Integral and Related Functional Differential Equations. Cambridge University Press. Cambridge. 2004. 597 p.
  * [15] Polyanin, A.D., Manzhirov, A.V. Handbook of Integral Equations: Second Edition. Chapman and Hall/CRC. New York. 2008.