Symmetrical Augmented System of Equations for the Parameter Identiﬁcation of Discrete Fractional Systems by Generalized Total Least Squares

: This paper is devoted to the identiﬁcation of the parameters of discrete fractional systems with errors in variables. Estimates of the parameters of such systems can be obtained using generalized total least squares (GTLS). A GTLS problem can be reduced to a total least squares (TLS) problem. A total least squares problem is often ill-conditioned. To solve a TLS problem, a classical algorithm based on ﬁnding the right singular vector or an algorithm based on an augmented system of equations with complex coefﬁcients can be applied. In this paper, a new augmented system of equations with real coefﬁcients is proposed to solve TLS problems. A symmetrical augmented system of equations was applied to the parameter identiﬁcation of discrete fractional systems. The simulation results showed that the use of the proposed symmetrical augmented system of equations can shorten the time for solving such problems. It was also shown that the proposed system can have a smaller condition number.


Introduction
Equations with fractional-order derivatives and differences are widely used to describe various processes and phenomena. Models based on equations with fractional-order derivatives and differences find wide application in hydrology, economics, and forecasts of network data usage [1][2][3].
Most of the research in this field deals with the parameter identification of fractionalorder differential equations using the error in an equation or an output error. References [4,5] described time domain identification methods based on the least squares technique. References [6][7][8] presented an overview of different methods for identifying systems with fractional-order derivatives. The use of total least squares for estimating the parameters of fractional-order differential equations was proposed in [9,10].
Autoregressions with fractional differences are widely used in the analysis of time series with long memory [11,12]. There are a large number of different models with generalizations of fractional differences, such as Gegenbauer autoregressive moving average (GARMA) [13,14], fractional ARUMA [15], seasonal autoregressive fractionally integrated moving average (SARFIMA) [16,17], and autoregressive tempered fractionally integrated moving average (ARTFIMA) [18,19]. Various aspects of using fractional differences for time series analysis have been considered [20,21].
In [22], a discrete Kalman filter with fractional differences was proposed. Equations with fractional-order differences were also used to approximate differential equations with Grunwald-Letnikov derivatives [23]. Meanwhile, the modeling of various physical processes based on equations with fractional-order differences was considered in [24][25][26].

1.
The problem is sensitive: The minimum singularities have a small gap; 2.
The problem is large: The SVD of a large matrix is very expensive. Instead, we can approximate the secular equation in the large-scale problem by Gaussian quadrature rules.
When applying the secular equation approach, it is necessary to solve a biased normal system. A biased normal system has more conditionality than a solution based on the SVD approach.
The problem of estimating the parameters of fractional systems is often ill-conditioned. To reduce the condition number, augmented systems of equations can be used [37,38].
For TLS, [39] proposed an augmented system of equations with complex coefficients equivalent to a biased normal system.
In article [34], an algorithm was proposed for the parameter identification of discrete fractional systems based on an augmented system of equations.
This article proposes a symmetrical augmented system of equations with real coefficients for the parameter identification of discrete fractional systems by generalized total least squares. The proposed symmetrical augmented system of equations has a smaller dimension than the augmented system [39], as well as real coefficients.

Problem Statement
The dynamical system is described by linear stochastic equations with fractional differences: where b (m) and a (m) are constant coefficients; 0 < α 1 . .
The following assumptions are introduced: A1. The dynamic system (1) is asymptotically stable. Results on the asymptotic stability of discrete time fractional difference systems can be found in [40,41].
A2. The noise sequences {ξ i }, {ζ i } are independent sequences with E(ξ i ) = 0, where E is the expectation operator. A3. The input and output sequences { x i } , { z i } are independent of noise sequences A4. The input signal x i persistently excites the order r 1 [42]. A5. The noise variance ratio γ = σ 2 ζ /σ 2 ξ is known. It is required to estimate the unknown coefficients of the dynamic system described by Equation (1) from the observed sequences {y i }, {w i } with noise for the known orders r, r 1 , α m , and β m .
If r, r 1 , α m , and β m are unknown, it is necessary to apply algorithms based on global optimization, such as genetic algorithms [43].

The Objective Function for Parameter Estimation
In [13], the following objective function was proposed for estimating the parameters: where: Theorem 1. Let the dynamic system described by Equation (1) with initial zero conditions and Assumptions A1-A5 be introduced. Then, the estimate of the coefficients determined by expression (2) exists, is unique, and converges to the true value of the coefficients with probability 1, i.e., θ → a.s.

θ.
The proof of the theorem is similar to the proof given in [33]. The objective function (2) is a generalized total least squares objective function. In [32], an iterative procedure based on the minimization of the generalized Rayleigh ratio was proposed. The simulation results in [34] showed that for ill-conditioned systems, the algorithm has a low accuracy of parameter estimation.
We transformed the objective function (2) into the following form: Under Assumptions A1-A5, the matrix H ξζ is nonsingular. The generalized total least squares problem (3) can be reduced to the total least squares problem [44].
Then objective function (3) can be defined as: The generalized total least squares algorithm for a nonsingular matrix H ξζ has this form (Algorithm 1). Algorithm 1. Ref. [34].
Step 1. Change variables H ξζ = U T ξζ U ξζ and ϑ = U ξζ θ Step 2. Calculate the estimate of the vector of parameterθ by the objective function (5).

Numerical Method for the Total Least Squares Problem
There are many algorithms for solving total least squares problems. The classical algorithm for solving the total least squares problem is based on SVD (singular value decomposition) [45]. The solution to a total least squares problem based on augmented systems of the equation was considered in [40]. To solve a large dimension or with a sparse matrix, iterative methods for solving total least squares problems are used, such as Newton's method [46,47], Rayleigh iterations [48], and Lanschoz iterations [49].
The problem of finding the singular vector is a nonlinear vector problem. Solving this problem numerically involves significant difficulties [50] related to convergence issues, high computational complexity, and the stability of search algorithms. Another approach is based on solving a biased normal system. Reference [51] showed that given the satisfaction of the condition:  (5) can be found as a solution to the biased normal system of equations When system (6) is solved, only the scalar problem of finding the minimal singular value σ remains nonlinear. This problem is always well-conditioned [48]. The vector problem is linear, but the system (6) is often ill-conditioned. The spectral condition number of the biased normal matrix is obtainable from the expression: The use of Cholesky decomposition can make the solution more stable. As this method applies to systems of linear equations with symmetrical, positively definite matrices, it can also be used for the biased normal system (6). However, Cholesky decomposition for ill-conditioned matrices leads to large errors.

Augmented System of Equations for Total Least Squares
In [39], an approach was proposed for solving normal systems of equations based on equivalent augmented systems of equations. Normal systems of equations corresponding to ordinary least squares: Can be represented as an equivalent augmented system of equations: Mathematics 2021, 9, 3250 5 of 13 (7) is defined as: In [38], a generalization of the augmented system (7) to the case of the biased normal system (6) is proposed. It was proposed to use a complex matrix G and vector g defined as: Using notation (8), the biased normal system (6) can be defined as: Or: where k is an arbitrary positive variable factor: The number of conditions (9) is defined as [17]: where µ max and µ min are the maximal and minimal eigenvalues of a matrix respectively. The problem of finding the minimum value of the condition number can be considered as a problem of choosing the optimal factor k: Problem (9) has no analytical solution but can be solved by numerical methods. In [38], an estimate was proposed for factork opt : Mathematics 2021, 9, 3250 6 of 13 Equation (7) can be solved by standard methods for solving systems of equations, such as using LU decomposition [52].
Generalized total least squares with the augmented system of Equation (9) in the second step algorithm were proposed for estimating the parameters of the system (1) [15].
The simulation results in [15,39] showed that the use of the augmented system of equations (7) is superior in accuracy to methods based on minimization of the generalized Rayleigh ratio or Cholesky decomposition for the biased normal system of equations.
The limitations of the augmented system (9) include the following: 1.
The augmented system of equations has complex coefficients, even if the original system has real coefficients.

2.
The matrix A(k) is non-Hermitian. It is impossible to apply LDL T decomposition for the solution of the system of Equation (9) [52]. 3.
The matrix A(k) has the dimensions 2(r + r 1 ) + N.

4.
It is necessary to determine the optimal factor k opt .

Symmetry Augmented System of Equations for Total Least Squares
Let us obtain an alternative augmented system of equations free of these limitations. We transformed the system of Equation (6): We combined Equation (13) with the equation Let us represent this in matrix form: Or: The resulting augmented system has a real symmetric matrix. The dimension of the matrix A is r + r 1 + N.

Theorem 2.
The spectral condition number of the system of Equation (15) is defined as: Proof of Theorem 2. We introduced a new variable, ρ = σ −1 ε, and considered the eigenvalue problem: We transformed the system of Equation (17): From Equation (18), it follows that if υ = 0, then the eigenvector υ and λ 2 − 2λσ are the corresponding eigenvalue of the matrix where the eigenvalues of the matrix by solving the quadratic equation are: Solving Equation (19), we can obtain λ i = σ ± σ i . Then, the spectral condition number is:

Results
Test cases were compared to the normalized root mean square error (NRMSE) of parameter estimation, defined as: The results were based on 50 independent Monte-Carlo simulations. The number of data points N in each simulation was 200. The input signal was a step function. The dynamic system is described by the equation: Figure 1 shows the step response of the system (20).

Example 1.
The noise-to-signal ratios were σ ξ /σ z = 10 −12 and σ ζ /σ x = 10 −12 . The gap between the minimal singular values of the matrix CU −1 ξζ d is g = 4.5199 ± 0.0048. Table 1 shows the mean values of the condition numbers of the matrices and their standard deviations.
8246.5 ± 886.6 90.684 ± 4.839 215.410 ± 13.147 Figure 2 shows the histogram of the decision time. Table 2 shows the mean values of the decision time and their standard deviations. Figure 3 shows the histogram of the normalized root mean square error of the parameters. The results were based on 50 independent Monte-Carlo simulations. The number of data points N in each simulation was 200. The input signal was a step function. The dynamic system is described by the equation:

Example 1.
The noise-to-signal ratios were 12 10 z ξ σ σ − = and 12 10 .  Table 1 shows the mean values of the condition numbers of the matrices and their standard deviations.       Figure 3 shows the histogram of the normalized root mean square error of the parameters.   Table 3 shows the mean values of the normalized root mean square error of the parameter estimation and their standard deviations.  Figure 3 shows the histogram of the normalized root mean square error of the parameters.  Table 3 shows the mean values of the normalized root mean square error of the parameter estimation and their standard deviations.
The gap between the minimal singular values of the matrix is g = (7.628 ± 0.4828)10 −12 . Table 4 shows the mean values of the condition numbers of the matrices and their standard deviations.
Mathematics 2021, 9, x FOR PEER REVIEW 10 of 14  . Table 4 shows the mean values of the condition numbers of the matrices and their standard deviations. ( ) (1.6530 ± 0.1216)10 38 (4.0269 ± 0.1828)10 18 (6.4410 ± 0.1510)10 9 Figure 4 shows the histogram of the decision time. Table 5 shows the mean values of the decision time and their standard deviations.      Table 5 shows the mean values of the decision time and their standard deviations.  Figure 5 shows the histogram of the normalized root mean square error of the parameters.  Table 6 shows the mean values of the normalized root mean square error of parameter estimation and their standard deviations. In Table 3, the standard deviations and mean values of the normalized root mean square error of the parameters are nearly equal. This is due to the fact that the condition numbers were large enough and the resulting estimates were sensitive to the initial data.
In Table 6, the standard deviations for the classical algorithm significantly exceed the mean values of the normalized root mean square error of the parameters. This indicates that in some experiments the results are very inaccurate, suggesting that the algorithm is numerically unstable and cannot be applied without regularization methods.
As can be seen from Examples 1,2, conditionality is highly dependent on the signalto-noise ratios upon input and output. Let us investigate how the condition number changes with a change in the signal-to-noise ratio upon input and output. Upon output, there is a signal-noise ratio of As can be seen from Figure 6, the conditionality of the proposed system is lower than the conditionality of the system (8).  Table 6 shows the mean values of the normalized root mean square error of parameter estimation and their standard deviations. In Table 3, the standard deviations and mean values of the normalized root mean square error of the parameters are nearly equal. This is due to the fact that the condition numbers were large enough and the resulting estimates were sensitive to the initial data.
In Table 6, the standard deviations for the classical algorithm significantly exceed the mean values of the normalized root mean square error of the parameters. This indicates that in some experiments the results are very inaccurate, suggesting that the algorithm is numerically unstable and cannot be applied without regularization methods.
As can be seen from Examples 1,2, conditionality is highly dependent on the signal-tonoise ratios upon input and output. Let us investigate how the condition number changes with a change in the signal-to-noise ratio upon input and output. Upon output, there is a signal-noise ratio of σ ξ /σ z = 10 −12 , while there is a signal-to-noise ratio upon input of σ (m) ζ /σ x = 10 −(12−0.2m) . Figure 2 shows lg κ 2 A and lg(κ 2 (A(k))). As can be seen from Figure 6, the conditionality of the proposed system is lower than the conditionality of the system (8).

Discussion
This paper proposed an estimation method of the parameters of discrete fractional systems by generalized total least squares. The simulation results showed that the parameter estimates obtained with the augmented systems of the equation can exceed the accuracy of the SVD-based estimates. This is due to the fact that methods based on augmented systems are less sensitive to the size of the gap between the smallest singular numbers.
In Table 6, showing the standard deviations for the classical algorithm, the mean values of normalized root mean square error of parameters is significantly exceeded. This shows that the method based on the SVD approach in example 2 is numerically unstable.
The solution time for the proposed method is less than that used in [34]. It is also not required to determine the optimal parameter opt k . The LDL T decomposition-based implementation turned out to be faster than the LU decomposition-based implementation. The application of the proposed augmented system is especially justified in cases when it is necessary to repeatedly estimate the parameters of discrete fractional systems, such as in algorithms that estimate both linear and nonlinear parameters [43]. In the future, it is planned to investigate the use of regularized versions of these systems.

Discussion
This paper proposed an estimation method of the parameters of discrete fractional systems by generalized total least squares. The simulation results showed that the parameter estimates obtained with the augmented systems of the equation can exceed the accuracy of the SVD-based estimates. This is due to the fact that methods based on augmented systems are less sensitive to the size of the gap between the smallest singular numbers.
In Table 6, showing the standard deviations for the classical algorithm, the mean values of normalized root mean square error of parameters is significantly exceeded. This shows that the method based on the SVD approach in example 2 is numerically unstable.
The solution time for the proposed method is less than that used in [34]. It is also not required to determine the optimal parameter k opt . The LDL T decomposition-based implementation turned out to be faster than the LU decomposition-based implementation.
The application of the proposed augmented system is especially justified in cases when it is necessary to repeatedly estimate the parameters of discrete fractional systems, such as in algorithms that estimate both linear and nonlinear parameters [43]. In the future, it is planned to investigate the use of regularized versions of these systems.