Parameter Identification of the Discrete-Time Stochastic Systems with Multiplicative and Additive Noises Using the UD-Based State Sensitivity Evaluation

: The paper proposes a new method for solving the parameter identification problem for a class of discrete-time linear stochastic systems with multiplicative and additive noises using a numerical gradient-based optimization. The constructed method is based on the application of a covariance UD filter for the above systems and an original method for evaluating state sensitivities within the numerically stable, matrix-orthogonal MWGS transformation. In addition to the numerical stability of the proposed algorithm to machine roundoff errors due to the application of the MWGS-UD orthogonalization procedure at each step, the main advantage of the obtained results is the possibility of analytical calculation of derivatives at a given value of the identified parameter without the need to use finite-difference methods. Numerical experiments demonstrate how the obtained results can be applied to solve the parameter identification problem for the considered stochastic system model.


Introduction
This paper addresses a parameter identification problem for a class of discrete-time linear stochastic systems represented by state-space difference equations with additive and multiplicative noises.A distinctive feature of this class of dynamic systems is multiplicative noise, which can be included in both the state and measurement equations [1].The reasons for the appearance of multiplicative noise are different; for example, it can be due to linearization, quantization, and modeling errors, or physical phenomena such as fading in communication channels.Most often, systems with additive and multiplicative noises are considered when solving problems related to various kinds of measurement processing.
The parameter identification problem consists of determining the unknown parameters of a mathematical model of the system belonging to a selected class of models using known input and output measurement data [2].As noted in [3,4], the main approaches to solving parameter identification problems, which remain of current interest today, are subspace identification methods and minimum prediction error (MPE) methods.The first approach is based on the application of projections in Euclidean space, and the second one is based on minimizing the identification criterion depending on the system parameters.The foundation of these approaches was laid in the seminal works [5,6].A great contribution to the development of MPE methods was made by Lennart Ljung, who defined the basic concepts of MPE methodology [2,4].
Linear discrete-time stochastic systems with associated Kalman-type filtering algorithms have been extensively used in practice.As a rule, application of the filter equations assumes a complete a priori knowledge of the system model parameters, but it is a rare case.The classical way of solving the parameter identification problem is to use adaptive filters where the model parameters are estimated together with the dynamic state [7].This requires determination of the sensitivities of the system state to unknown parameters, i.e., partial derivatives of the state estimates.Straightforward differentiation of the filter equations is a direct approach to compute the state sensitivities.This leads to a set of filter sensitivity equations.It is well known that, for discrete-time linear stochastic systems with additive noises, a conventional Kalman filter may suffer from numerical instability caused by machine roundoff errors [8].This is also true for such systems with multiplicative noises [9].Nevertheless, there are currently many numerically stable modifications of the conventional Kalman filtering algorithm, which are based on various methods for factorizing the estimation error covariance matrices [8,10].It is worth noting that, for discrete-time stochastic systems with multiplicative noises, such modifications were developed relatively recently (see, for example, [9,11]).
In this paper, we chose a numerically stable UD-based modification of the Kalmantype filtering algorithm [9] for solving the parameter identification problem, since it has the following attractive numerical properties: numerical stability to machine roundoff errors, absence of square root and matrix inversion operations, and the compactness and simplicity of the block matrix array form for implementation on a computer, including parallel computing [10].
Thus, the purpose of our work is to develop a numerically stable UD-based method for solving the parameter identification problem for a class of discrete-time linear stochastic systems with multiplicative and additive noises using a numerical, gradient-based optimization of the identification criterion in the form of a negative logarithmic likelihood function.
We propose the following solution steps to achieve the goal: 1.
Replace the conventional Kalman-type filtering algorithm with a UD-based covariance filtering algorithm that is numerically stable to machine roundoff errors.

2.
Construct a new method for calculating the state sensitivity values in the adaptive UD-based filter.

3.
Apply this method to the gradient-based minimization of the identification criterion.
The paper is organized as follows.Section 2 provides the problem statement of parameter identification for a considered class of stochastic systems.Then, the UD-based covariance filtering algorithm for discrete-time stochastic systems with multiplicative and additive noises is described.Lemma 1 provides the method of obtaining the main results presented in the next section.Section 3 contains the new UD-based state sensitivity evaluation method, which is described in Proposition 1. Proposition 2 explains how the values of the identification criterion and its gradient can be calculated using the suggested algorithm.Section 4 demonstrates how these two methods can be applied for solving the parameter identification problem of the considered stochastic system model.Section 5 concludes the paper.

The Problem Statement
Consider a discrete-time linear stochastic system with multiplicative and additive noises where x k ∈ R n is the system state vector; z k ∈ R m is the measurement vector; matrices F k , The considered systems find their application in wireless sensor networks, navigation, and other fields.Multiplicative noises are commonly used to account for stochastic uncertainties in system dynamics and measurements.In this paper, we assume that system (1) additionally contains parametric uncertainty, i.e., the matrices defining the equations of system (1) depend on the unknown parameter θ ∈ R p .Consequently, all elements of the system matrices F k , F k , G k , H k , H k , noise covariance matrices Q k and R k , and variances σ 2 ξ , σ 2  ζ , as well as the initial conditions x0 and Π 0 , may depend on the unknown parameter θ; i.e., Thus, the parameter identification problem of system (1) arises.It is known that the performance of the Kalman filter degrades in the presence of modeling uncertainties.Solving the problem of parameter identification makes it possible to cope with this problem.The well-known approaches involve numerical minimization by θ of the identification criterion [12] θmin = argmin θ∈D(θ) which implies solving an optimization problem with constraints, where D(θ) represents some compact set, such as a segment in the scalar case.
To solve the problem of parameter identification, we chose an identification criterion (2) in the form of the negative logarithmic likelihood function [13]: depending on the available measurement information, which includes both measurements Z M 1 = {z 1 , . . ., z k , . . ., z M } themselves and the observed measurement residuals ν k (θ) = z k − H k xk (θ), calculated by conventional Kalman-type filtering algorithm (Algorithm 1 in [9]).Here, In this paper, our first task is to replace the conventional Kalman-type filtering algorithm with a numerically stable, UD-based covariance filtering algorithm.

The UD-Based Covariance Filtering Algorithm for Discrete-Time Stochastic Systems with Multiplicative and Additive Noises
Originally, two algorithms based on the UDU T decomposition of covariance and information matrices for discrete-time stochastic systems with multiplicative and additive noises were proposed in [9].The first of these algorithms is the UD-based covariance filter, and the second one is the UD-based information filter.Both of these filters have an extended array form, and their computational schemes allow for updating all required filter quantities with the use of the numerically stable MWGS orthogonalization procedure.
The UD-based implementations imply the decomposition of the error covariance matrix in the form of P = U P D P U T P , where U P is an upper triangular matrix with 1's on the main diagonal, and D P is a diagonal matrix.To recursively update the resulting UDfactors U P and D P , we use the modified weighted Gram-Schmidt UD-based (MWGS-UD) orthogonalization [14] as follows: given a pair of the pre-arrays {A, D A }, compute a pair of the post-arrays {U , D U } by means of the MWGS-UD orthogonalization, i.e., where A ∈ R r×s , r > s, and M ∈ R r×s is the MWGS-UD transformation that produces the block upper triangular matrix U ∈ R s×s .The diagonal matrices D A ∈ R r×r and D U ∈ R s×s satisfy M T D A M = D U and D A > 0 (see Lemma VI.4.1 in [14] for an extended explanation).
In this paper, we consider a modification of the UD-based covariance filter [9], whose equations are presented in Algorithm 1.The proof of Algorithm 1 is similar to the proof presented in [9].

End For
Output: xk , Remark 1. Steps 3 and 7 of Algorithm 1 require the application of the UDU T modified Cholesky decomposition [15].Steps 4-6, 8, and 9 require the application of the MWGS-UD orthogonalization to a pair of block matrices ⟨A T i , D A i ⟩ that produces another pair of block matrices ⟨U i , D U i ⟩ (i = 1, . . ., 5) so that the equalities (4) are satisfied.
Let us rewrite the identification criterion (3) in terms of the UD-CF algorithm.Taking into account that det , we obtain where diagonal matrix D Σ k and normalized residual vector ēk are calculated in Steps 9 and 10 of Algorithm 1.

Derivative Evaluation of the MWGS-Based Array of Block Matrices
Solving the parameter identification problem requires minimization of the identification criterion (5) with respect to unknown system parameters.It is often done by using the gradient approach [7], where the computation of ∇ θ J UD (θ; Z M 1 ) is necessary.For the discrete-time stochastic system (1), the J UD (θ; Z M 1 ) and ∇ θ J UD (θ; Z M 1 ) evaluation demands an implementation of the UD-CF and of the so-called "differentiated" UD-CF to determine the state sensitivities of the system state to the unknown system parameters [7,16,17].The ∇ θ J UD (θ; Z M 1 ) computation will lead to a set of p filter sensitivity equations for computing ∂ xk /∂θ and a set of p matrix Riccati-type sensitivity equations for computing ∂P k /∂θ.Such a method of state sensitivity evaluation within the UD-based array covariance filter was proposed in [17] for a class of discrete-time linear stochastic systems with only additive noises.In this paper, we extend this approach to systems with multiplicative and additive noises.Solving the problem of state sensitivity evaluation, we augment the numerical scheme of the UD-CF (Algorithm 1) with a procedure for numerically efficient evaluation of the derivatives of the UD-based filter variables with respect to unknown system parameters.
We use the following basic result of [17] that gives a simple and convenient technique which naturally augments any MWGS-based array of block matrices for computing derivatives of its elements.

Lemma 1 ([17]
).Let entries of the pre-arrays A, D A in (4) be known differentiable functions of a parameter θ.Consider the transformation in (4).Given the derivatives of the pre-arrays A ′ θ and (D U ) ′ θ , the following formulas calculate the corresponding derivatives of the post-arrays: where the quantities L0 , D 0 , and Ū0 are, respectively, the strictly lower triangular, diagonal, and strictly upper triangular parts of the matrix product M T D A A ′ θ U −T .Additionally, D 2 and Ū2 are the diagonal and strictly upper triangular parts of the product M T (D A ) ′ θ M, respectively.

Main Results
One can see from Algorithm 1 that elements xk and the UD factors U P k and D P k are readily available from this UD-based filter.Hence, our aim is to augment equations of the UD-CF so that the derivatives ∂ xk (θ)/∂θ i and ∂U P k (θ) /∂θ i , ∂D P k (θ) /∂θ i , i = 1, . . ., p can be computed using quantities available from this UD-CF algorithm.

The New UD-Based State Sensitivity Evaluation Method
Now, we are ready to present our new result-the UD-based state sensitivity evaluation method for a class of discrete-time linear stochastic systems with multiplicative and additive noises.Proposition 1.Let the elements of matrices defining system (1) be known differentiable functions of a parameter θ.Then, for a given value of parameter θ, the estimates of state vector x k , their sensitivity values ∂ xk /∂θ, and the UD factors U P k , D P k of the error covariance matrices and their sensitivity values ∂U P k /∂θ, ∂D P k /∂θ can be evaluated simultaneously using the subsequent Algorithm 2.
Proof.Steps 1, 3, and 7 of Algorithm 2 are obtained from the corresponding steps of Algorithm 1 by direct differentiation of the vector x0 and the UD factors of the covariance matrices P 0 , Q k−1 , and R k with respect to parameter θ.
Steps 2, 10, and 11 of Algorithm 2 are obtained from the corresponding steps of Algorithm 1 by direct differentiation of the equations with respect to parameter θ.
Let us consider the equations that define Steps 4-6, 8, and 9 in Algorithm 2. They have the same form and are obtained using Lemma 1 as follows.At each step, the MWGS-UD procedure orthogonalizes the columns of the block matrix A i with respect to the weight matrix D A i so that the equalities in (4) are satisfied.To calculate sensitivity values (U i ) ′ θ and (D U i ) ′ θ for a given value of parameter θ, it is necessary and sufficient to find the values of the partial derivatives of the elements of the block matrices A i and D A i , i.e., calculate the matrices (A i ) ′ θ and (D A i ) ′ θ .Lemma 1 with A := A i , D A := D A i , and U := U i , D U := D U i (i = 1, . . ., 5) allows for achieving this goal.As a result, we find block matrices (U i ) ′ θ and (D U i ) ′ θ , from which we obtain the required sensitivity values.Thus, Proposition 1 gives the method for state sensitivity evaluation using the UDbased covariance filter for discrete-time linear stochastic systems with multiplicative and additive noises.

The UD-Based Computation of the Identification Criterion and Its Gradient
Algorithm 2 can be considered as an adaptive filter, in which parameter θ is adjusted according to the minimum of criterion J UD (θ; Z M 1 ).When solving the parameter identification problem using gradient-based algorithms, adaptive filters are used, supplemented by a sensitivity model equations [18] to calculate the gradient of the identification criterion.
Implementation of gradient-based numerical methods for minimizing identification criterion (5) requires calculating the values of its gradient ∇ θ J UD (θ; Z M 1 ).The conventional gradient-based method has the following iterative form: where θ j is the parameter vector at the jth iteration.In (7), ∇ θ denotes the gradient operator , which is applied here to (5) at point θ = θ j .Scalar step size parameter β j is designed to ensure that J UD (θ j+1 ; Z M 1 ) ≤ J UD (θ j ; Z M 1 ) + e, where e is a positive number that can be chosen in a variety of ways [19].
Let us write the equation we could use to evaluate the gradient of the identification criterion (3) in terms of the UD-CF algorithm.Let θ ∈ R p .Then, from (5), we have Taking into account matrix differentiation rules, we rewrite (8) to obtain the expression for ∇ θ J UD (θ; Z M 1 ) evaluation (i = 1, . . ., p): Proposition 2. Let the elements of matrices defining system (1) be known differentiable functions of a parameter θ ∈ R p .Then, for a given value of parameter θ, the values of identification criterion J UD (θ; Z M 1 ) and its gradient ∇ θ J UD (θ; Z M 1 ) can be evaluated according to Equations ( 5) and ( 9) using Algorithm 2.
Thus, for a given value of parameter θ, Algorithm 2 allows us to obtain all the quantities necessary to calculate the values of the identification criterion and its gradient that are used in the parameter identification gradient-based method.

Discussion
As a practical application of the proposed method, let us consider the parameter identification problem for a nearly constant velocity model of the uniform motion augmented with multiplicative noises [11]: where ), and θ is the model parameter to be identified.Let us put the "true" value of the parameter equal to θ * = 0.3.
To demonstrate the validity of the proposed approach, we have conducted numerical experiments in MATLAB, which is a common tool for simulation, optimization, control, and filtering [8,20,21].We have implemented all necessary functions for simulating system dynamics and measurements, as well as functions for calculating the identification criterion (5) and its gradient (9). Figure 1 shows the identification criterion and its gradient for the considered problem, with noise covariance matrix R = 0.5 2 I 2 averaged over 100 runs.The following MATLAB functions were used for numerical minimization of the identification criterion: simulannealbnd, ga, and fmincon.The first two functions implement gradient-free metaheuristic algorithms SA (Simulated Annealing) and GA (Genetic Algorithm), respectively.The third function was configured to use two different gradient-based algorithms: interior-point (IP) and trust-region-reflective (TRR).The first algorithm estimates the gradient using finite differences, and the second algorithm uses a user-provided gradient of the objective function.
A series of 100 numerical experiments was conducted for each value of the noise level σ.In each experiment, numerical identification of parameter θ was performed based on

Conclusions
In this paper, we have proposed a new method for solving the parameter identification problem for discrete-time linear stochastic systems with multiplicative and additive noises based on the application of the covariance UD-filter and the original method for the state sensitivity evaluation within the numerically stable, matrix-orthogonal MWGS-UD transformation.
The main theoretical results of the paper are the UD-based state sensitivity evaluation method and the method of calculating the values of identification criterion J UD (θ; Z M 1 ) and its gradient, which are formulated in Propositions 1 and 2, respectively.Both methods use Algorithm 2 to calculate the state vector estimates and their sensitivity values.In addition to the numerical stability of the MWGS orthogonalization procedure to machine roundoff errors, the main advantage of the proposed method is the possibility of analytical calculation of derivative values at a given value of the identified parameter without the need to use finite-difference methods.
Numerical experiments confirm that the obtained results can be used for solving the parameter identification problems of the considered stochastic systems.The gradient-based minimization of the identification criterion that uses the proposed method outperforms both non-gradient algorithms and the algorithm that estimates gradients using finite differences.
is the number of measurements;x 0 ∼N ( x0 , Π 0 ) is the initial state; ξ k ∈ R∼N (0, σ 2 ξ ) and w k ∈ R q ∼N (0, Q k ) aremultiplicative and additive noises in the state equation, respectively; ζ k ∈ R∼N (0, σ 2 ζ ) and v k ∈ R m ∼N (0, R k ) are multiplicative and additive noises in the measurement equation, respectively; covariance matrices Q k and R k of noises w k and v k , respectively, are positive definite; and all noises and the initial state are mutually independent.