A Bayes Estimator of Parameters of Nonlinear Dynamic Systems

Abstract

A new multipolynomial approximations algorithm (the MPA algorithm) is proposed for estimating the state vector 𝜃 of virtually any dynamical (evolutionary) system. The input of the algorithm consists of discrete-time observations 𝑌. An adjustment of the algorithm is required to the generation of arrays of random sequences of state vectors and observations scalars corresponding to a given sequence of time instants. The distributions of the random factors (vectors of the initial states and random perturbations of the system, scalars of random observational errors) can be arbitrary but have to be prescribed beforehand. The output of the algorithm is a vector polynomial series with respect to products of nonnegative integer powers of the results of real observations or some functions of these results. The sum of the powers does not exceed some given integer 𝑑. The series is a vector polynomial approximation of the vector 𝐸(𝜃∣𝑌), which is the conditional expectation of the vector under evaluation (or given functions of the components of that vector). The vector coefficients of the polynomial series are constructed in such a way that the approximation errors uniformly tend to zero as the integer 𝑑 increases. These coefficients are found by the Monte-Carlo method and a process of recurrent calculations that do not require matrix inversion.

1. Introduction

Consider a dynamical system whose mathematical model has the form 𝑥(𝑘+1)=𝑓(𝑥(𝑘),𝜂_𝑘,𝑘). (1.1) Here, the integers 0,1,…,𝑘,… refer to time instants 𝑡₀,𝑡₁,…,𝑡_𝑘,…; 𝑥(𝑘)∈𝑅^𝑛 is the state vector of the dynamical system at the instant 𝑡_𝑘; 𝜂_𝑘∈𝑅^𝑛 is the vector of random perturbations; 𝑓(𝑥(𝑘),𝜂_𝑘,𝑘) is a given function (nonlinearly depending on its arguments, in general); 𝑥(0),𝜂_𝑘 are random vectors with given distributions.

Let the results of observations at the instants 𝑡₁,…,𝑡_𝑘,…,𝑡_𝑁 be described by a scalar sequence 𝑦₁,…,𝑦_𝑘,…,𝑦_𝑁 generated by the mathematical model 𝑦_𝑘=𝐻_𝑘(𝑥(𝑘),𝜉_𝑘), (1.2) where 𝐻(𝑘,𝑥(𝑘)) is a given function of its arguments, and the distribution is given for the sequence of random variables 𝜉₁,…,𝜉_𝑘,…,𝜉_𝑁.

In what follows, it is assumed that 𝑦₁,…,𝑦_𝑘 are components of the vector 𝑌_𝑘.

A parameter vector 𝜃 is defined as a vector whose components coincide with some components (or given functions of these) of state vectors of the dynamical system at given time instants. In the present paper, we describe a new algorithm for estimating the vector 𝜃.

The problem of estimating the vector 𝜃 is called [1, 2] the problem of smoothing, if 𝜃=𝑥(0); the problem of filtration, if 𝜃=𝑥(𝑁); the problem of extrapolation, if 𝜃=𝑥(𝑁∗),𝑁∗>𝑁.

It is well known that the mean-square optimal estimate of the random vector 𝜃 on the basis of observations of the vector 𝑌_𝑁 is the vector of conditional expectation 𝐸(𝜃∣𝑌_𝑁). Therefore, we consider the problem of creating the MPA algorithm for the construction of approximations converging to the vector 𝐸(𝜃∣𝑌_𝑁).

There is only one requirement with regard to computer power and the expressions (1.1) and (1.2): for given distributions of random variables 𝑥(0),𝜂_𝑘,𝜉_𝑖, it is possible to generate a set that consists of a considerable number (several thousands) of sequences of random state vectors and observation results satisfying (1.1) and (1.2).

Elements of the sequence 𝑦₁,…,𝑦_𝑘,…,𝑦_𝑁 are sent to the input of the algorithm at the instants 𝑡₁,…,𝑡_𝑘,…,𝑡_𝑁.

Let us briefly consider some known algorithms that give an approximate solution to the above estimation problems.

(1) The algorithm of the nonlinear least-squares method should be regarded as the most general method of solution. This algorithm is mentioned in a huge number of publications such as, for instance, [3]. Thus, in the smoothing problem for 𝜂(𝑘)≡0, the algorithm determines the estimate vector 𝜃=𝑧_𝑥(0)(𝑌_𝑁) that approximately minimizes the heuristic quality functional 𝑆(𝑌_𝑁,𝑧_𝑥(0)(𝑌_𝑁)) with the constraints (1.1). The estimation quality functional is the sum of squared residuals 𝑆(𝑌_𝑁,𝑧_𝑥(0)(𝑌_𝑁))=^𝑖=𝑁_𝑖=1(𝑦_𝑖−𝐻(𝑖,𝑥_𝑖(𝑧_𝑥(0)(𝑌_𝑁))))². (1.3) Here, the vectors 𝑥_𝑖(𝑧_𝑥(0)(𝑌_𝑁)) are obtained by consecutive calculations of state vectors on the basis of (1.1) with the initial data 𝑧_𝑥(0)(𝑌_𝑁). For the numerical solution of this problem one utilizes numerous versions of the gradient method or the Newton method. In order to apply the Newton method, one has to solve a system of algebraic equations, which is obtained if we equate to zero partial derivatives of the right-hand side of (1.3) with respect to the components of the vector 𝑧_𝑥0(𝑌_𝑁).

The basic drawbacks of the nonlinear least-squares method are the following:

(i)the quality functional in the form of the sum of squared residuals is theoretically adequate for an estimate quality criterion only in the unrealistic case of the residuals being independent random variables with Gaussian distribution;(ii)the versions of the gradient method or the Newton method are applicable if the functions in the mathematical models (1.1) and (1.2) are differentiable; these methods require the existence of a good first approximation and the construction of a complicated interactive computation process, if the function 𝑆(𝑌_𝑁,𝑧_𝑥0(𝑌_𝑁)) happens to have a set of local minima;(iii)to get an idea of the statistical characteristics of estimation errors is possible only by means of the Monte-Carlo method (there are no explicit formulas for the estimate error variance).

(2) The extended Kalman filter algorithm (EKF algorithm) yields an heuristic solution to the filtration problem (there are many publications dedicated to the exposition of the EKF algorithm; it suffices to mention [4]). The EKF algorithm is based on a sequence of linearizations of nonlinear functions involved in the mathematical model of a dynamical system, the linearizations being constructed in a neighborhood of a sequence of estimate vectors.

The recurrent scheme of the EKF algorithm has a two-step structure.

Step 1. Suppose that after the instant 𝑡_𝑘−1, approximations have been found for the first and the second statistical moments of the state vector at the instant 𝑡_𝑘−1 (the estimation vector for this vector coincides with the first statistical moment). On the time interval [𝑡_𝑘−1,𝑡_𝑘−0] (prior to observation results at the instant 𝑡_𝑘) prediction is made and its results are approximations for the first and the second statistical moments of the state vector at the instant 𝑡_𝑘−0. The vector and the matrix of the prediction are found by the calculation of a matrix formed by partial derivatives of the state vector (the Jacobian matrix). Therefore, realization of the prediction requires the differentiability of the right-hand side of (1.1) and the admissibility of its linearization with respect to the components of the state vector increments arising from 𝑡_𝑘−1 to 𝑡_𝑘.

Step 2. After actual observations at the instant 𝑡_𝑘 and linearization of the right-hand side of (1.2), mean-square optimal linear correction of the first statistical moment of the state vector at the instant 𝑡_𝑘 is implemented, as well as the corresponding correction of the second statistical moment. This step results in approximations for the first and the second statistical moments of the state vector components at the instant 𝑡_𝑘.

The two-step scheme of the recurrent EKF algorithm is quite credible and gives quite satisfactory results in many practically important cases.

The basic drawbacks of the EKF algorithm are similar to those of the nonlinear least-squares method:

(i)the estimate quality functional is defined and coincides with the estimation error variance vector only if the functions in the right-hand sides of (1.1) and (1.2) linearly depend on their arguments (in this situation, the EKF turns into the standard discrete Kalman filter);(ii)the EKF algorithm is applicable only if the functions involved in the mathematical models (1.1) and (1.2) are differentiable;(iii)if the functions in (1.1) and (1.2) are essentially nonlinear, it is necessary to have a fairly good first approximation in order to ensure a convergent calculation process (the notion of a “fairly good approximation’’ is heuristic and is determined on the basis of numerical experiments);(iv)there are no explicit formulas for the estimate error variance.

(3) In last ten years the significant number of researches with representation basic solution of a problem by definition of conditional expectation of a vector of parameters for mathematical model nonlinear dynamic system [1, 5–10] was published. By means of multiple application of formula Bayes to vectors from (1.1) and (1.2) and by using numerical quadratures, the recurrent equations are easyly discovered for the probability density function (pdf) of state vectors of dynamic system. The actual solutions of the previous recurrent equations are, however, unfeasible because of the unwieldy dimensions of the integrals involved [10]. Therefore several alternative strategies have been developed. One set of such alternatives consists of developing suboptimal filtering, such as those based on linearization or transformations, and the other consists of methods which employ Monte-Carlo simulation strategies (e.g., estimation by using the sequential importance sampling particle filter) to approximate evaluate the multidimensional integrals in a recursive manner.

The presented direction is perspective but is bulky and unsuitable for the solution of practical applied (instead of model!) problems of nonlinear identification.

2. Schematic Diagram of MPA Algorithm

Let us describe the basic structure of the MPA algorithm [11–14]. Denote by 𝜃 the vector of dimension 𝑟×1 to be estimated upon fixing the vector 𝑌_𝑁.

The MPA algorithm is a new recurrent algorithm, such as asymptotic accurately solve nonlinear problem to construction of the vectors conditional expectations of the vector of parameters for mathematical model nonlinear dynamic system by any random errors and perturbations with given distributions. Therefore MPA algorithm approximate solve problem the optimal mean-square estimation.

MPA algorithm in principle differs from the all above mentioned algorithms since

(i)MPA algorithm does not make linearizations and do not attempt to find the posterior probability density function of state vectors of dynamic system and never uses Bayes's formula, although essentially MPA algorithm is Bayes—it uses a priori information about a priori distributions. For given 𝑌_𝑁, MPA algorithm once build approximations to vectors posterior expectation of the estimated parameters vector. The approximation determined in the form segments of the vectorial power series relatively products the integer non negative powers of components of the vector 𝑌_𝑁. For every term of series sum of powers bounded above given integer 𝑑. For given 𝑌_𝑁, the approximations uniformly bounded converge to the vectors posterior expectation of the estimated parameters vector by 𝑑→∞;(ii)the vectorial coefficients of power series are defined by means of an adjustment of MPA algorithm. The adjustment is realized be means of define the first and second statistical moments of basis random vector, which components are estimated parameters and the mentioned powers of components of the vector 𝑌_𝑁. This definition is realized by means of Monte-Carlo's method from (1.1) and (1.2). However, after these vectorial coefficients have been calculated and kept in computer memory, the estimation vectors are determined by simple calculations after any new observations of the vectors 𝑌_𝑁;(iii)there is a sequence 𝑟 vectors such as each consist from 𝑘 leading to a component of vector 𝜃, 𝑘=1,…,𝑟. After adjustment by given 𝑑,𝑁, MPA algorithm determines the recurrent approximations to the conditional expectation of these vectors and matrices covariance of errors of these approximations. In this connection, formulas of evaluations are simple and do not require the inverse of matrices and evaluations of quadratures.

The new MPA algorithm proposed in this paper for estimating the vector of parameter (and in particular state vector) of a nonlinear dynamical system, for the most part, has no drawbacks mentioned above.

Step 1. Suppose that 𝑑 is a given positive integer number and the set of integer numbers 𝑎₁,…,𝑎_𝑁 consists of all nonnegative solutions of the integer inequality 𝑎₁+⋯+𝑎_𝑁≤𝑑, the number of which we denote by 𝑚(𝑑,𝑁). The value 𝑚(𝑑,𝑁) is given by the recurrent formula proved by induction, 𝑘𝑚(𝑑,𝑁)=𝑚(𝑑−1,𝑁)+(𝑁+𝑑−1)⋯𝑑!,𝑚(1,𝑁)=𝑁. (2.1) We obtain the vector 𝑊_𝑁(𝑑) of dimension 𝑚(𝑑,𝑁)×1, the components 𝑤1,…,𝑤_{𝑚(𝑑,𝑁)} of which are all possible values of the form that represent the powers of measurable values.

Next, we define a basic vector 𝑉(𝑑,𝑁) of dimension (𝑟+𝑚(𝑑,𝑁))×1,𝑉(𝑑,𝑁)=‖𝜃𝑊_𝑁(𝑑)‖^𝑇.

Step 2. We use a known statistical generator of random vectors 𝑥(0),𝜂_𝑘,𝜉_𝑖 to solve repeatedly the Cauchy problem for (1.1) for a given initial conditions 𝑥(0), a control law 𝑢 and various realizations of random vectors 𝑥(0),𝜂_𝑘,𝜉_𝑖 and 𝑘. As a result, the computer memory will contain a set of realizations of random basic vectors 𝑉(𝑑,𝑁) sufficient for the calculation of the statistical characteristics of the basic vector 𝑉(𝑑,𝑁). We apply the Monte-Carlo method to find the prior first and second statistical moments of the vector 𝑉(𝑑,𝑁), that is, the mathematical expectation 𝑉(𝑑,𝑁)=𝐸(𝑉(𝑑,𝑁)), and the covariance matrix 𝐶_𝑉(𝑑,𝑁)=𝐸((𝑉(𝑑,𝑁)−𝐸(𝑉(𝑑,𝑁)))(𝑉(𝑑,𝑁)−𝐸(𝑉(𝑑,𝑁)))^𝑇). (2.2)

Implementation of Step 2 is a learning process for the algorithm, adjusting it to solve the particular problem described by (1.1) and (2.3).

Step 3. For given and 𝑁 and a fixed vector 𝑌_𝑁, we assign the vector 𝑍_{𝐸(𝜃∣𝑌𝑁)}(𝑊_𝑁(𝑑)) to be the solution to the estimation problem. This vector gives an approximate estimate of the vector 𝐸(𝜃∣𝑌_𝑁) that is optimal in the root-mean-square sense on the set of vector linear combinations of components of the vector 𝑊_𝑁(𝑑). Then the vector 𝑍_{𝐸(𝜃∣𝑌𝑁)}(𝑊_𝑁(𝑑)) is an element of this set and is an approximate mean-square optimal estimate for the vector 𝐸(𝜃∣𝑌_𝑁) of the conditional expectation of the vector 𝜃.
The vector 𝑍_{𝐸(𝜃∣𝑌𝑁)}(𝑊_𝑁(𝑑)) can be presented in a following aspect: 𝑍_{𝐸(𝜃∣𝑌𝑁)}(𝑊_𝑁(𝑑))=_{𝑎1+⋯+𝑎𝑁≤𝑑}𝜆(𝑎₁,…,𝑎_𝑁)𝑦^𝑎1₁⋯𝑦^𝑎𝑁_𝑁, (2.3) where vectors 𝜆(𝑎₁,…,𝑎_𝑁) are some weight vectorial factors.

The vector 𝑉(𝑑,𝑁)=𝐸(𝑉(𝑑,𝑁)) and the matrix 𝐶_𝑉(𝑑,𝑁) are the initial conditions for the process of recurrent calculations that realizes the principle of observation decomposition [6] and consists of 𝑚(𝑑,𝑁) steps. Once the final step is performed, we obtain vector weight coefficients 𝜆(𝑎₁,…,𝑎_𝑁) for (1.1). Moreover, we determine the matrix 𝐶(𝑑,𝑁), which is the covariance matrix of the estimation errors for the vector of conditional mathematical expectation estimated by the vector 𝜃. Calculating the elements of the matrix 𝐶(𝑑,𝑁), we have the method of preliminary (prior to the actual flight) analysis of observability of identified parameters for the given control law, structure of measurements and their expected random errors. Recurrent calculations do not require matrix inversion and indicate the situations when the next component of the vector 𝑊_𝑁(𝑑) is close to linear combination of its previous components. To implement the recursion, we process the components of the vector 𝑊_𝑁(𝑑) one after another. However, the adjustment of the algorithm performed by applying the Monte-Carlo method to find the vector 𝑉(𝑑,𝑁)=𝐸(𝑉(𝑑,𝑁)) and the matrix 𝐶_𝑉(𝑑,𝑁) takes into account a priori ideas on stochastic structure of components of the whole set of possible vectors 𝑊_𝑁(𝑑) that can appear in any realizations of the random vectors 𝑥(0),𝜂_𝑘,𝜉_𝑖 and 𝑘 allowed by a priori conditions. This adjustment is the price we have to pay if we want the MPA algorithm to solve nonlinear identification problems efficiently. This is what makes the MPA algorithm differs fundamentally from, for instance, the standard Kalman filter designed to solve linear identification problems only or from multiple variations of algorithms resulted from attempts to extend the Kalman filter to nonlinear filtration problems.

Using the multidimensional analog of the Weierstrass theorem (a corollary of the Stone theorem [15]), we prove that with the increase of the integer 𝑑 the error estimate vector 𝑍_{𝐸(𝜃∣𝑌𝑁)}(𝑊_𝑁(𝑑))−𝐸(𝜃∣𝑌_𝑁) tends to zero uniformly in some domain.

The basic structure of the MPA algorithm described above allows us to use a strict quality criterion for the estimates obtained, since on each step the algorithm tends to ensure the minimal mean-square error of the estimate for the conditional expectation vector with a given volume of observations. Its closeness to the minimum increases with the increase of the integer 𝑑 and the increase of the number of realizations, if the Monte-Carlo method is used.

If the number 𝑚(𝑑,𝑁) of the approximating polynomial series is not small, a large volume of calculations has to be performed (after the initial choice of the integers 𝑛,𝑁,𝑑) for the determination of the vector coefficients 𝜆(𝑎₁,…,𝑎_𝑁).

3. Mean-Square Optimal Estimate for the Conditional Expectation Vector 𝐸(𝜃∣𝑌_𝑁)

Consider the principal algorithm (PA) for solving the problem of the mean-square optimal estimate for the vector 𝐸(𝜃∣𝑌_𝑁). It is known that the vector 𝐸(𝜃∣𝑌_𝑁) is the mean-square optimal estimate for the vector 𝜃 after the vector 𝑌_𝑁 has been fixed. Therefore, it makes sense that the PA should estimate the conditional expectation vector.

We are going to construct a PA that gives a polynomial approximation of the vector 𝐸(𝜃∣𝑌_𝑁). To that end, we obtain an estimate for the vector 𝐸(𝜃∣𝑌_𝑁) which is linear with respect to the components of the vector 𝑊_𝑁(𝑑) and is mean-square optimal.

In what follows the vector of that estimate is denoted by 𝑍_{𝐸(𝜃∣𝑌𝑁)}(𝑊_𝑁(𝑑)).

An explicit expression for the estimate vector is obtained after the calculation of the elements of the vector 𝑉(𝑑,𝑁)=𝐸(𝑉(𝑑,𝑁)) and the covariance matrix 𝐶_𝑉(𝑑,𝑁), which coincide with the first and the second (centered) statistical moments for the vector 𝑉(𝑑,𝑁). This vector and this matrix can be divided into blocks whose structure can be represented as follows: 𝐸(𝐸(𝜃∣𝑌_𝑁))=𝐸(𝜃),𝐸(𝐸(𝑊(𝑑,𝑁)∣𝑌_𝑁𝐶))=𝐸(𝑊(𝑑,𝑁)),_𝜃(0)=𝐸((𝐸(𝜃∣𝑌_𝑁)−𝐸(𝜃))(𝐸(𝜃∣𝑌_𝑁)−𝐸(𝜃))^𝑇)=𝐸((𝜃−𝐸(𝜃))(𝜃−𝐸(𝜃))^𝑇𝐿),_𝑁(𝑑)=𝐸((𝐸(𝜃𝑌_𝑁)−𝐸(𝐸(𝜃∣𝑌_𝑁(𝑑)))(𝑊_𝑁(𝑑)−𝐸(𝑊_𝑁(𝑑)))^𝑇)=𝐸(𝜃)𝑊_𝑁(𝑑)^𝑇)−𝐸(𝜃)𝐸(𝑊_𝑁(𝑑))^𝑇,𝑄_𝑁(𝑑)=𝐸((𝑊_𝑁(𝑑)−𝐸(𝑊_𝑁(𝑑)))(𝑊_𝑁(𝑑)−𝐸(𝑊_𝑁(𝑑)))^𝑇). (3.1) The right-hand sides of the above blocks coincide with the first and the second (centered) statistical moments determined by the Monte-Carlo method. On the other hand, the left-hand sides of the same blocks coincide with the first and the second (centered) statistical moments of the components of the conditional expectation vector. Therefore, the said statistical moments can be found experimentally on the basis of mathematical models (1.1) and (1.2) also for the conditional expectation vectors. This obvious statement is crucial for the practical numerical procedure of estimating the conditional expectation vector 𝐸(𝜃∣𝑌_𝑁).

Let 𝑍_{𝐸(𝜃∣𝑌𝑁)}(𝑊_𝑁(𝑑))=𝐸(𝜃)+Λ_𝑁(𝑑)(𝑊_𝑁(𝑑)−𝐸(𝑊_𝑁(𝑑))), (3.2) where Λ_𝑁(𝑑) is an 𝑛×𝑚(𝑑,𝑁)-matrix satisfying the equation Λ_𝑁(𝑑)𝑄_𝑁(𝑑)=𝐿_𝑁(𝑑). (3.3)

Let 𝑍_𝜃(𝑊_𝑁(𝑑))=𝑧+Λ_𝑁(𝑑)(𝑊_𝑁(𝑑)−𝐸(𝑊_𝑁(𝑑))), (3.4) where 𝑧 is an arbitrary 𝑛×1-vector and Λ_𝑁(𝑑) is an arbitrary 𝑛×𝑚(𝑑,𝑁)-matrix. Let 𝐶_𝑍(𝑑,𝑁) and 𝐶(𝑑,𝑁) be error estimate covariance matrices for the vector 𝐸(𝜃∣𝑌_𝑁) in terms of the estimate vectors 𝑍_{𝐸(𝜃∣𝑌𝑁)}(𝑊_𝑁(𝑑)) and 𝑍_𝜃(𝑊_𝑁(𝑑)).

Lemma 3.1. A matrix 𝐶(𝑑,𝑁)−𝐶_𝑍(𝑑,𝑁) is the nonnegative definite matrix: 𝐶_𝑍(𝑑,𝑁)≤𝐶(𝑑,𝑁).

Lemma 3.1 follows from the identity 𝐶(𝑑,𝑁)=𝐶_𝑍(𝑑,𝑁)+(Λ_𝑁(𝑑)−Λ_𝑁(𝑑))(Λ_𝑁(𝑑)−Λ_𝑁(𝑑))^𝑇+(Λ_𝑁(𝑑)𝑄_𝑁(𝑑)−𝐿_𝑁(𝑑))(Λ_𝑁(𝑑)−Λ_𝑁(𝑑))^𝑇+((Λ_𝑁(𝑑)𝑄_𝑁(𝑑)−𝐿_𝑁(𝑑))(Λ_𝑁(𝑑)−Λ_𝑁(𝑑))^𝑇)^𝑇+(𝑧−𝐸(𝜃)(𝑧−𝐸𝜃)^𝑇). (3.5)

Corollary from Lemma 3.1
The vector 𝑍_{𝐸(𝜃|𝑌𝑁)}(𝑊_𝑁(𝑑)) is a mean-square optimal estimate for the vector 𝐸(𝜃∣𝑌_𝑁) on the set of estimates linear with respect to the components of the vector 𝑊_𝑁(𝑑).

If 𝑄_𝑁(𝑑)>0, then the estimate vector is unique and 𝑍_{𝐸(𝜃∣𝑌𝑁)}(𝑊_𝑁(𝑑))=𝐸(𝜃)+𝐿_𝑁(𝑑)𝑄_𝑁(𝑑)⁻¹(𝑊_𝑁(𝑑)−𝐸(𝑊_𝑁(𝑑))). (3.6)

The covariance matrix 𝐶_𝑍(𝑑,𝑁) for the error of the estimate of vector 𝐸(𝜃∣𝑌_𝑁) is defined by 𝐶_𝑍(𝑑,𝑁)=𝐶_𝜃(0)−Λ_𝑁(𝑑)𝐿_𝑁(𝑑). (3.7)

If 𝑄_𝑁(𝑑)≥0, then the vectors giving the linear mean-square optimal estimate are nonunique, but the variance of the components of the difference of these vectors is equal to zero.

Formula (3.2) gives explicit expressions for the vector coefficients 𝜆(𝑎₁,…,𝑎_𝑁) in (2.3). These expressions are obtained by equating the right-hand side of (3.2) to the right-hand side of (2.3) and writing out explicit expressions for the components of the vector 𝑊_𝑁(𝑑).

Let us examine asymptotic errors of the estimates obtained when using formula (3.2).

For a given vector 𝑌_𝑁, suppose that the vector 𝐸(𝜃∣𝑌_𝑁) is defined as a function of the argument 𝑌_𝑁 in some a priori given compact domain Ω_𝑌𝑁⊂𝑅^𝑁 and is continuous in that domain. Then the following result holds.

Theorem 3.2. sup_{𝑌𝑁∈Ω𝑌𝑁}|𝑍_{𝐸(𝜃∣𝑌𝑁)}(𝑊_𝑁(𝑑))−𝐸(𝜃∣𝑌_𝑁)|⇒0,𝑑⇒∞. (3.8)

Proof. According to the multidimensional version of the Weierstrass theorem [7], for any 𝜀>0 there is a multidimensional polynomial 𝑃(𝑊_𝑁(𝑑_𝜀)) such that sup_{𝑌𝑁∈Ω𝑌𝑁}|𝑃(𝑊_𝑁(𝑑_𝜀))−𝐸(𝜃∣𝑌_𝑁)|<𝜀. (3.9) This relation can be rewritten as sup_{𝑌𝑁∈Ω𝑌𝑁}|𝑃(𝑊_𝑁(𝑑))−𝐸(𝜃∣𝑌_𝑁)|⇒0,𝑑⇒∞. (3.10) Let 𝐶 be the covariance matrix for the random vector 𝑃(𝑊_𝑁(𝑑))−𝐸(𝜃∣𝑌_𝑁): 𝐶=𝐸(𝑃(𝑊_𝑁(𝑑))−𝐸(𝜃∣𝑌_𝑁))(𝑃(𝑊_𝑁(𝑑))−𝐸(𝜃∣𝑌_𝑁))^𝑇. (3.11) From (3.10), it follows that 𝐶⇒0_𝑛,𝑑⇒∞. (3.12) By construction, the vector 𝑍_{𝐸(𝜃∣𝑌𝑁)}(𝑊_𝑁(𝑑)) gives a mean-square optimal estimate for the vector 𝜃 and this estimate is linear with respect to the components of the vector 𝑊_𝑁(𝑑). But Lemma 3.1 implies that for any other nonoptimal linear estimate, in particular that like the vector 𝑃(𝑊_𝑁(𝑑)), we have 𝐶≥𝐶_𝑍(𝑑,𝑁). Hence, taking into account (3.12), we get 𝐶_𝑍(𝑑,𝑁)⇒0_𝑛,𝑑⇒∞. (3.13) The statement (3.13) is equivalent to (3.8), if we take into account that 𝐶_𝑍(𝑑,𝑁)=(𝑍_{𝐸(𝜃∣𝑌𝑁)}(𝑊_𝑁(𝑑))−𝐸(𝜃∣𝑌_𝑁))×(𝑍_{𝐸(𝜃∣𝑌𝑁)}(𝑊_𝑁(𝑑))−𝐸(𝜃∣𝑌_𝑁))^𝑇𝑝(𝜃,𝑌_𝑁)𝑑𝜃𝑑𝑌_𝑁, (3.14) where 𝑝(𝜃,𝑌_𝑁) is the joint distribution density for the random vectors 𝜃,𝑌_𝑁. Theorem 3.2 is proved.

Thus, the PA, by virtue of (2.3), determines a vector series which, with the increase of the number of its terms 𝑚(𝑑,𝑁), approximates the conditional expectation for the vector 𝜃 of the estimated parameters with arbitrarily small uniform mean-square error.

4. Recurrently MPA Algorithm

Suppose that for a given observation vector 𝑌_𝑁 and an integer 𝑑 we have constructed the vector 𝑊_𝑁(𝑑) with the components 𝑤₁,𝑤₂,…,𝑤_{𝑚(𝑑,𝑁)}.

The implementation of the algorithm starts with approximate calculation (by the Monte-Carlo method) of 𝑛+𝑚(𝑑,𝑛)-dimensional integrals that determine the components of the vector 𝑉(𝑑,𝑁)=𝐸(𝑉(𝑑,𝑁)) and the matrix 𝐶_𝑉(𝑑,𝑁)=𝐸((𝑉(𝑑,𝑁)−𝑉(𝑑,𝑁))(𝑉(𝑑,𝑁)−𝑉(𝑑,𝑁))^𝑇). Further, splitting this vector and this matrix into suitable blocks and using (2.3), we find the desired estimate vector 𝑍_{𝐸(𝜃∣𝑌𝑁𝑘)}(𝑊_𝑁(𝑑)) and the estimate error covariation matrix 𝐶_𝑍(𝑑,𝑁).

However, this approach is unreasonable, since it requires the inversion of the matrix 𝑄_𝑁(𝑑), which is a difficult task for 𝑄_𝑁(𝑑) of a large dimension 𝑚(𝑑,𝑁)×𝑚(𝑑,𝑁) or 𝑄_𝑁(𝑑) being close to a singular matrix.

Below we describe a recurrent calculation process based on the principle of decomposition of observations expounded in [11].

We are going to construct a recurrently MPA algorithm that does not involve matrix inversion and consists of 𝑚(𝑑,𝑁) steps of calculations of the first and the second statistical moments for a sequence of special vectors 𝑉₁,…,𝑉_𝑖,…,𝑉_{𝑚(𝑑,𝑁)} after a priori moments 𝑉(𝑑,𝑁),𝐶_𝑉(𝑑,𝑁) for the basic vector 𝑉(𝑑,𝑁) have been found.

Denote by 𝑉₁ the vector formed by the components of the basic vector 𝑉(𝑑,𝑁) remaining after the elimination of the component 𝑤₁;…; by 𝑉_𝑖+1 the vector formed by the components of 𝑉_𝑖 remaining after the elimination of 𝑤_𝑖; and so forth.

The quantity 𝑤_{𝑚(𝑑,𝑁)} is the last component of the vector 𝑊_𝑁(𝑑), and after its elimination, the last vector 𝑉_{𝑚(𝑑,𝑁)} turns out equal to the estimate vector 𝑍_{𝐸(𝜃∣𝑌𝑁)}(𝑊_𝑁(𝑑)).

On Step 1, formulas (3.2) and (3.6) are used for the calculation of the following quantities:

(i)the vector 𝑉₁ that yields an estimate for the vector 𝑉₁ which is mean-square optimal and linear with respect to 𝑤₁;(ii)the estimate error covariance matrix 𝐶(𝑉₁).

The estimate vector consists of the conditional expectation estimate vector 𝐸(𝜃∣𝑤₁) and a vector of dimension (𝑚(𝑑,𝑁)−1)×1. For 𝑤₁ fixed, the latter coincides with the vector of statistical prediction of mean future values of the quantities 𝑤₂,…,𝑤_{𝑚(𝑑,𝑁)}.

Note that the calculations of Step 1 are based on a priori values 𝑉(𝑑,𝑁),𝐶(𝑉(𝑑,𝑁)) found beforehand.

In a similar way, suppose that on steps 1,…,𝑖 of the calculation process, the vector 𝑉_𝑖 and the matrix 𝐶(𝑉_𝑖) have been found after fixing the quantities 𝑤₁,…,𝑤_𝑖.

On step 𝑖+1, the following quantities are calculated by means of (3.2) and (3.6):

(i)the vector 𝑉_(𝑖+1) of the estimate of 𝑉_𝑖+1, which is mean-square optimal and linear with respect to 𝑤₁,…,𝑤_𝑖+1;(ii)the estimate error covariance matrix 𝐶(𝑉_(𝑖+1)).

The vector 𝑉_(𝑖+1) is composed of the conditional expectation estimate 𝐸(𝜃∣𝑤₁,…,𝑤_𝑖+1) and a vector which, after the quantities 𝑤₁,…,𝑤_𝑖+1 have been fixed, is the vector of statistical prediction for the mean values of “future” quantities 𝑤_𝑖+2,…,𝑤_{𝑚(𝑑,𝑁)}.

Note that the calculations of step 𝑖+1 are based on the already determined quantities 𝑉_𝑖,𝐶(𝑉_𝑖), which it is natural to call a priori first and second statistical moments of the “future’’ quantities 𝑤_𝑖+1,…,𝑤_{𝑚(𝑑,𝑁)}. The vector 𝑉_𝑖 and the matrix 𝐶(𝑉_𝑖) represent a priori data regarding the statistical moments of the components of the vector 𝑉_𝑖+1 before the quantity 𝑤_𝑖+1 is sent to the input of the algorithm.

The formulas of the computation precess described above are the following: 𝑉_(𝑖+1)=𝑉¹_𝑖+𝑞_𝑖⁻¹𝑏_𝑖(𝑤_𝑖+1−𝑧_𝑤𝑖+1),𝐶(𝑉_𝑖+1)=𝐶(𝑉_𝑖)¹−𝑞_𝑖⁻¹𝑏_𝑖𝑏^𝑇_𝑖. (4.1) Here, the scalar 𝑧_𝑤𝑖+1 is the (𝑛+1)th component of the vector 𝑉_𝑖 (this scalar is a linear mean-square optimal estimate for the component 𝑤_𝑖+1 after the algorithm has used the components 𝑤₁,…,𝑤_𝑖); the vector 𝑉¹_𝑖 is obtained from the vector 𝑉_𝑖 by eliminating its component 𝑧_𝑤𝑖+1; the scalar 𝑞_𝑖 is the (𝑛+1)th diagonal element of the matrix 𝐶(𝑉_𝑖) (this scalar is the variance of the estimate error for the component 𝑤_𝑖+1 after the components 𝑤₁,…,𝑤_𝑖 have been used); the matrix 𝐶(𝑉_𝑖)¹ is obtained from the matrix 𝐶(𝑉_𝑖) by elimination of its (𝑛+1)th row and its (𝑛+1)th column; the vector 𝑏_𝑖 coincides with the (𝑛+1)th column of the matrix 𝐶(𝑉_𝑖) with the (𝑛+1)th component excluded.

If the scalar 𝑞_𝑖 happens to be close to zero, then the component 𝑤_𝑖+1 becomes close to a linear combination of the components 𝑤₁,…,𝑤_𝑖. In this case, the component 𝑤_𝑖+1 gives no new information about 𝜃 and should be dropped from the computation process..

Note that the sequence of random variables (𝑤_𝑖+1−𝑧_𝑤𝑖+1) is a renewable sequence.

The top left block of the (𝑛×𝑛)-matrix 𝐶(𝑉_𝑖) coincides with the estimate error covariance matrix 𝐶_𝑍(𝑑,𝑖) for the vector 𝐸(𝜃∣𝑤₁,…,𝑤_𝑖) after the vector 𝑊_𝑖(𝑑) has been used by the algorithm.

Denote by 𝑙(𝑖) the vector formed by the first 𝑛 components of the vector 𝑏_𝑖. The formula representing the evolution of the covariance matrix 𝐶_𝑍(𝑑,𝑖) as a function of 𝑖 (the number of the components of the vector 𝑊_𝑖(𝑑)) has the form 𝐶_𝑍(𝑑,𝑖)=𝐶_𝜃(0)−𝑞₁⁻¹𝑙(1)𝑙(1)^𝑇−⋯−𝑞⁻¹_𝑖−1𝑙(𝑖−1)𝑙(𝑖−1)^𝑇. (4.2)

In what follows, the above MPA algorithm is tested in several numerical problems of estimating components of state vectors of essentially nonlinear dynamical systems. The components to be estimated are unknown random constant parameters 𝜃₁,…,𝜃_𝑟 of a dynamical system. The specific applied problems considered below are those of smoothing and filtration.

The Monte-Carlo method is used if the number of random realizations is from 5000 to 10000. This number has no strong effect on estimate errors arising in the MPA algorithm. It is assumed that the estimated random parameters are statistically independent and have uniform distribution a priori.

The quantity 𝜎(𝑖,cal) is determined theoretically by the calculation of variances, that is, the diagonal elements of the covariance matrix 𝐶_𝑍(𝑑,𝑁).

The quantity 𝜎(𝑖,exp) is determined experimentally by the Monte-Carlo method with the number of realizations being 10000.

The fact that the experimental and the theoretical mean-square deviations are approximately the same proves the correctness of the MPA algorithm formulas specified above.

5. Estimates for the Principal Moments of Inertia and Orientation Angles of the Principal Axes of Inertia of a Solid

Consider a solid in space (e.g., a satellite) with a fixed reference trihedron, that is, orthogonal axes 𝑋, 𝑌, 𝑍 and their origin 𝑂. The corresponding projections 𝜔_𝑥,𝜔_𝑦,𝜔_𝑧 of the angular velocity of the body, 𝜔, are measured by angular velocity sensors. The projections 𝑀_𝑥,𝑀_𝑦,𝑀_𝑧 of the moment 𝑀 of external forces are assumed known.

In the general situation, a solid has 6 moments of inertia relative to a rectangular reference frame 𝑋,𝑌,𝑍, among which three are called centrifugal moments. It is known that for a point 𝑂 there is a rectangular coordinate system 𝑋_𝑜,𝑌_𝑜,𝑍_𝑜 whose axes coincide with the principal axes of inertia. In this coordinate system, the centrifugal moments of inertia are equal to zero and the remaining ones, 𝐽_𝑥,𝐽_𝑦,𝐽_𝑧, are called the principal moments of inertia relative to the point 𝑂.

Let 𝛼,𝛽,𝛾 be the Euler angles that determine the angular position of 𝑋_𝑜,𝑌_𝑜,𝑍_𝑜 with respect to 𝑋,𝑌,𝑍, let 𝐴 be the orthogonal matrix of the direction cosines of these angles.

Generally, a solid in outer space is a complex nonsymmetric structure whose parts play different roles and are constructed of different materials. Therefore, it is difficult to calculate the quantities 𝛼,𝛽,𝛾,𝐽_𝑥,𝐽_𝑦,𝐽_𝑧 on the basis of design diagrams of these parts and their junctures. However, for the analysis of the dynamics of a rotating solid it is convenient to use the coordinate axes that coincide with the principal axes of inertia, since the equations of motion (dynamical Euler equations) have a simple form in these axes.

Thus, it is reasonable to state the problem as follows: using the MPA algorithm, find estimates for the quantities 𝜃₁=𝛼,𝜃₂=𝛽,𝜃₃=𝛾,𝜃₄=𝐽_𝑥,𝜃₅=𝐽_𝑦,𝜃₆=𝐽_𝑧, provided that the solid is rotating in space and in the reference frame the quantities 𝜔_𝑥,𝜔_𝑦,𝜔_𝑧,𝑀_𝑥,𝑀_𝑌,𝑀_𝑧 are observed as discrete-time functions.

In the continuous-time model, the equations of motion of a nonlinear dynamical system have the form 𝐽_𝑥𝑑𝜔^𝑜_𝑥𝑑𝑡+(𝐽_𝑧−𝐽_𝑦)𝜔^𝑜_𝑦𝜔^𝑜_𝑧=𝑀^𝑜_𝑥,𝐽_𝑦𝑑𝜔^𝑜_𝑦𝑑𝑡+(𝐽_𝑥−𝐽_𝑧)𝜔_𝑧𝜔_𝑥=𝑀^𝑜_𝑦,𝐽_𝑧𝑑𝜔^𝑜_𝑧𝑑𝑡+(𝐽_𝑦−𝐽_𝑥)𝜔^𝑜_𝑥𝜔^𝑜_𝑦=𝑀^𝑜_𝑧,𝑑𝐽_𝑥𝑑𝑡=0;𝑑𝐽_𝑦𝑑𝑡=0;𝑑𝐽_𝑧𝑑𝑡=0,𝑑𝛼𝑑𝑡=0;𝑑𝛽𝑑𝑡=0;𝑑𝛾𝑑𝑡=0. (5.1)(5.2) Here, 𝜔^𝑜_𝑥,𝜔^𝑜_𝑦,𝜔^𝑜_𝑧,𝑀^𝑜,𝑀^𝑜_𝑦,𝑀^𝑜_𝑧 are projections of the vectors 𝜔 and 𝑀 to the axes 𝑋_𝑜,𝑌_𝑜,𝑍_𝑜: 𝜔^𝑜𝑀=𝐴𝜔,_𝑜=𝐴𝑀. (5.3) Further, assume that observations of the time-functions 𝜔_𝑥,𝜔_𝑦,𝜔_𝑧,𝑀_𝑥,𝑀_𝑌,𝑀_𝑧 take place with intervals 0.1 second and have duration 10 seconds. In order to create a mathematical model (1.1) with discrete time, we use numerical integration by the fourth-order Runge-Kutta method with constant step 0.2 second between observation instants and the initial conditions 𝜔_𝑥(0)=1/20,𝜔_𝑦(0)=−1/20,𝜔_𝑧(0)=1/20and𝑀_𝑥=100kg,𝑀_𝑌=−100kg,𝑀_𝑧=100kg.

Setting 𝑦_𝑥(𝑡_𝑘)=𝜔_𝑥(𝑡_𝑘),𝑦_𝑦(𝑡_𝑘)=𝜔_𝑦(𝑡_𝑘),𝑦_𝑧(𝑡_𝑘)=𝜔_𝑧(𝑡_𝑘), we find that at the instant 𝑡_𝑘 the components of the initial observation vector 𝑦(𝑡_𝑘) are related to nonobservable quantities 𝛼,𝛽,𝛾 and nonobservable functions 𝜔^𝑜_𝑥,𝜔^𝑜_𝑦,𝜔^𝑜_𝑧 by 𝑦(𝑡_𝑘)=𝐴^𝑇𝜔^𝑜(𝑡_𝑘). (5.4)

Using the sequence 𝑦_𝑥(𝑡_𝑘),𝑦_𝑦(𝑡_𝑘),𝑦_𝑧(𝑡_𝑘) of primary observation vectors and (5.1), it is required to estimate the parameters 𝐽_𝑥,𝐽_𝑦,𝐽_𝑧,𝛼,𝛽,𝛾.

The secondary observations 𝑦_𝑥(𝑖),𝑦_𝑦(𝑖),𝑦_𝑧(𝑖),𝑖=1÷4, which serve as the input data of the estimation algorithm, are found by summing up the components of each 𝑙 vectors of consecutive primary observations: 𝑦_𝑥(𝑖)=^𝑗=𝑙_𝑗=1𝑦_𝑥(𝑡_{𝑙(𝑖−1)+𝑗}𝑦),_𝑦(𝑖)=^𝑗=𝑙_𝑗=1𝑦_𝑦(𝑡_{𝑙(𝑖−1)+𝑗}𝑦),_𝑧(𝑖)=^𝑗=𝑙_𝑗=1𝑦_𝑧(𝑡_{𝑙(𝑖−1)+𝑗}). (5.5) In our model, it is assumed that 𝑙=25 and the array of 300 primary observations has been converted to an array of 𝑁=12 scalar secondary observations.

For the operation of the MPA algorithm it is necessary to find a priori first and second statistical moments of the basic 18-dimensional vector formed by the estimated quantities 𝛼,𝛽,𝛾,𝐽_𝑥,𝐽_𝑦,𝐽_𝑧 and 𝑦_𝑥(𝑖),𝑦_𝑦(𝑖),𝑦_𝑧(𝑖),𝑖=1÷4. These data were found by the Monte-Carlo method under the assumption that the random parameters to be estimated have uniform distribution on a priori given segments: 𝐽−𝜌≤𝛼,𝛽,𝛾≤𝜌,_𝑜𝑥(1−𝜌)≤𝐽_𝑥≤𝐽^𝑜_𝑥𝐽(1+𝜌),_𝑜𝑦(1−𝜌)≤𝐽_𝑦≤𝐽^𝑜_𝑦𝐽(1+𝜌),_𝑜𝑧(1−𝜌)≤𝐽_𝑧≤𝐽^𝑜_𝑧(1+𝜌), (5.6) where 𝐽^𝑜_𝑥=4000,𝐽^𝑜_𝑦=20000,and𝐽^𝑜_𝑧=15000.

The quantity 𝜌 describes the range in which a priori errors of the estimated parameters may vary. In our model, the values of 𝜌 were the following: 𝜌=0.2 (large a priori errors), 𝜌=0.1 (medium a priori errors), and 𝜌=0.05 (small a priori errors).

Let 𝜎_𝑎(𝑖) be a priori mean-square deviation (MSD) of the scattering of the unknown parameter 𝜃_𝑖, and let 𝜎(𝑖) be an a posteriori MSD of the error of estimating this parameter found after the utilization of the MPA. In what follows, the quality of the estimate of the parameter 𝜃_𝑖 is characterized by the ratios 𝜎(𝑖)/𝜎_𝑎(𝑖),𝑖=1,…,6, which show how many times a posteriori (after observation) dispersal of the estimated parameter is smaller than a priori dispersal.

Below, for different values of 𝜌 we give the values of the said ratios for the MPA algorithm operating with 𝑑=2 (the estimate of 𝜃_𝑖 is realized by second-order polynomials with respect to 12 secondary scalar observations; the number of terms of the polynomials is equal to 𝑚(2,12)=90):

𝜌=0.2𝜎(1)/𝜎_𝑎(1)𝜎(2)/𝜎_𝑎(2)𝜎(3)/𝜎_𝑎(3)𝜎(4)/𝜎_𝑎(4)𝜎(5)/𝜎_𝑎(5)𝜎(6)/𝜎_𝑎(6)0.55680.04400.56140.10050.08620.0886 (5.7)

𝜌=0.1𝜎(1)/𝜎_𝑎(1)𝜎(2)/𝜎_𝑎(2)𝜎(3)/𝜎_𝑎(3)𝜎(4)/𝜎_𝑎(4)𝜎(5)/𝜎_𝑎(5)𝜎(6)/𝜎_𝑎(6)0.50690.00630.50790.01230.01640.0161 (5.8)

𝜌=0.05𝜎(1)/𝜎_𝑎(1)𝜎(2)/𝜎_𝑎(2)𝜎(3)/𝜎_𝑎(3)𝜎(4)/𝜎_𝑎(4)𝜎(5)/𝜎_𝑎(5)𝜎(6)/𝜎_𝑎(6)0.49830.00120.49850.00170.00280.0026 (5.9) These data show that for 𝑑=2 and the chosen characteristics for the vectors of primary and secondary observations, a posteriori dispersal of estimate errors is about 100 or more times less than a priori dispersal of the unknown parameters 𝛽,𝐽_𝑥,𝐽_𝑦,𝐽_𝑧, provided that a priori errors of the parameters belong to the classes of medium or small errors. A posteriori dispersal of the parameters 𝛼,𝛾 is only about 2 times less than a priori dispersal. To increase the accuracy of the estimates of these parameters, it seems necessary to use larger values of 𝑑, as well as other characteristics of the vectors of primary and secondary observations. It is also possible to use iterations that on each step decrease the intervals of possible a priori errors.

For the general reduction of estimate errors it seems reasonable to split a priori error intervals into smaller intervals. Thus, suppose that we have large a priori errors. In order to obtain estimate errors corresponding to the class of medium a priori errors it suffices to divide every large a priori interval into two parts and apply the MPA algorithm to each of the 64 resulting families of a priori intervals from the class of medium error intervals. The MPA algorithm calculation of the MSDs of a posteriori dispersal allows us to find a family of a priori intervals with the minimal values of these MSDs.

6. Estimates of the State Vector of a Satellite by Means of Handling the Information on Angles of Vising

Many results of astronomical observations depend on measurements and subsequent processing of digital data with regard to a sequence of viewing angles of space objects obtained in the optical frequency band. The viewing angles are measured by an optical device or on the basis of images on a photographic plate.

These angles are often the only available data regarding high-orbit sun-exposed satellites with onboard radio communication system failure (the energy of terrestrial radio band viewing facilities is insufficient because of the long distances to a satellite). Therefore, it makes sense to develop an algorithm that would give fairly precise estimates for state vectors of objects in space on the basis of processing the data about viewing angles.

For high-orbit satellites such data is obtained by specially constructed optical-electronic systems (e.g., the optical-electronic system Window [16] which monitors near-the-earth-space in optical frequency band and gives viewing information about space bodies at distances up to 40000 km).

In an inertial geocentric coordinate system, the equations of motion of a body whose mass is not very large have the form 𝑑𝑋_𝑖𝑑𝑡=𝑉_𝑖,𝑑𝑉_𝑖𝑑𝑡=−𝜇𝑋_𝑖𝑅. (6.1) Here, 𝑅=(𝑋₁2+𝑋₂2+𝑋₃2)^1/2, 𝜇 is the gravity constant, 𝑋_𝑖,𝑖=1÷3 are rectangular coordinates of the body.

The viewing angles are the two angles that determine (with respect to the earth) the direction of the viewing vector, that is, the direction of the straight line from the observer to an object in space.

We assume that at the current time instant 𝑡 the center of viewing is located at a point of the Earth surface with geographical coordinates 𝜆=𝜔𝑡 (longitude), 𝜑 (latitude), ℎ (altitude), where 𝜔 is the angular speed of the Earth rotation. In the rotating right geocentric coordinate system 𝑥_𝑖,𝑖=1÷3, the coordinates of the center of viewing are 𝑋^′₁=𝑟_𝑒cos(𝜑)−ℎsin(𝜑),𝑋^′₂=0,𝑋^′₃=𝑟_𝑒sin(𝜑)+ℎcos(𝜑), where 𝑟_𝑒 is the Earth radius.

The two angles of orientation of the viewing vector, in the local horizon plane of the viewing center, are represented by an angle 𝑦₁ and, in the local vertical plane, by an angle 𝑦₂. These angles are defined by 𝑦₁𝑏=arctan(𝑎),𝑦₂𝑐=arctan((𝑎2+𝑏2)^1/2), (6.2) where 𝑎=cos(𝜆)cos(𝜑)𝑋₁+sin(𝜆)cos(𝜑)𝑋₂+sin(𝜑)𝑋₃−𝑟_𝑒,𝑏=−sin(𝜆)𝑋₁+cos(𝜆)𝑋₂,𝑐=−sin(𝜑)cos(𝜆)𝑋₁+sin(𝜑)sin(𝜆)𝑋₂+sin(𝜑)𝑋₃−ℎ. (6.3)

In what follows, we assume that the algorithm processing the data on viewing angles should solve a filtration problem: in an inertial geocentric coordinate system at certain current time instants, the algorithm should estimate three rectangular coordinates of an object in space and three rectangular components of its velocity vector.

The possibility of creating such an algorithm is based on the existence or correlation between the varying viewing angles and varying rectangular coordinates of a space object. This correlation is due to the law of gravitation responsible for the gravity acceleration vector that changes the trajectory of its motion and the viewing angles and is a known function of rectangular coordinates of the object.

In the absence of gravity acceleration, it would be impossible to estimate the parameters of motion by measuring the viewing angles. Indeed, the law of variation of the viewing angles is the same for all bodies in uniform rectilinear motion with parallel velocity vectors and the same constant ratio of the distance and the velocity modulus. Therefore, in this situation the parameters of motion are unobservable and cannot be estimated on the basis of information about the viewing angles.

The problem of the algorithm for estimating current components of the state vector satisfying the system of six nonlinear differential equations (6.1) is a problem of nonlinear filtration. The input of the algorithm consists of a sequence of observation results (6.2) nonlinearly depending on the current values 𝑋_𝑖,𝑖=1÷3.

Let 𝜏 be the time interval between two consecutive observations of pairs of viewing angles of the space object. After the first 𝑘 observations, we construct integral observation results 𝑦_int,1,𝑦_int,2. Here, 𝑦_int,1 is the sum of the first 𝑘 fixed quantities 𝑦₁; similarly, 𝑦_int,2 is the sum of the first 𝑘 fixed quantities 𝑦₂. In like manner, we construct integral pairs 𝑦_int,3,𝑦_int,4,…,𝑦_{int,2𝑠−1},𝑦_int,2𝑠,…, and so forth. The sequence of these random variables is a sequence of inputs to the recurrent algorithm at the instants 𝑇,2𝑇,…,𝑟𝑇,…, where 𝑇=𝑘𝜏. The summation operator reduces the effect of random tracking errors on the accuracy of the estimates.

Choosing an integer 𝑠, we introduce the following definitions. Let us divide the sequence of pairs of integral observation results into a sequence of segments, each containing 𝑠 pairs.

The iteration vector 1 is the estimate vector obtained after the 1st pair segment has been sent to the input of the recurrent algorithm. The iteration vector 1 is determined at the instant 𝑠𝑇 after 𝑠𝑘 observations of viewing angle pairs.

The iteration vector 2 is the estimate vector obtained after sending to the input the second segment of the pairs of integral observation results. The iteration vector 2 is determined at the instant 2𝑠𝑇 after 2𝑠𝑘 observations.

Iteration vectors 3, 4, and so forth, are determined in a similar manner.

For the construction of iteration vector 1, it is necessary to find the vector and the matrix of a priori data: the first and the second statistical moments of the basic random vector 1 whose first six components coincide with the components of the actual (unknown) initial state vector of the space object; its other components have the form of products of integer powers 2𝑠 of integral observation results; the sum of these powers does not exceed the given integer 𝑑. For 𝑑=1, the dimension of the basic vector 1 is equal to 6+2𝑠; for 𝑑=2, the dimension of the basic vector 1 is equal to 6+2𝑠+(2𝑠+1)𝑠, and so forth.

A priori scattering of the vectors 𝑋(0),𝑉(0) determines a priori dispersal of the other components of the basic vector 1 by means of (6.1) and (6.2). Next, we set 𝑋(0)=𝑋(0)_𝑁+𝛿𝑋(0),𝑉(0)=𝑉(0)_𝑁+𝛿𝑉(0), where 𝑋(0)_𝑁 is the nominal initial radius-vector of the observed space object and 𝑉(0)_𝑁 is the nominal initial vector of its velocity. It is assumed that 𝑉(0)_𝑁 is the velocity vector of an object moving on a circular orbit of radius |𝑋(0)_𝑁|. Uniform scattering of the components of random vectors 𝛿𝑋(0),𝛿𝑉(0) is determined by the inequalities −𝜌|𝑋(0)|≤𝛿𝑋₁,𝛿𝑋₂,𝛿𝑋₃≤𝜌|𝑋(0)|,−𝜌|𝑉(0)|≤𝛿𝑉₁,𝛿𝑉₂,𝛿𝑉₃≤𝜌|𝑉(0)|, (6.4) where 𝜌=0.2.

As the modeling process goes on, the errors in the determination of the initial coordinates reach several thousands km, and the errors in the determination of the initial velocity vector reach several hundreds m/sec. In this situation, the extended Kalman filter based on linearization is inapplicable for estimation, since it diverges at a fast rate.

A priori first and second statistical moments of the components of the basic vector 1 are determined by integrals over domains Ω_𝑋,Ω_𝑉, which are calculated by the Monte-Carlo method.

After a priori statistical characteristics of the basic vector 1 have been defined, 2𝑠 random variables are sent to the input of the recurrent algorithm, namely, 2𝑠 integral observation results for the viewing angles on the time interval 𝑠𝑇. After 2𝑠 steps of calculations, the recurrent algorithm determines the iteration vector 1, which is the estimate vector 𝑍_{𝑋(𝑠𝑇)},𝑍_{𝑉(𝑠𝑇)} for the state vector of the space object at the instant 𝑠𝑇.

The accuracy of the estimate is described by the estimate error covariance matrix calculated on each step of the recurrent algorithm.

The determination of iteration vector 2 is similar to the above. Set 𝑋(0)=𝑍_{𝑋(𝑠𝑇)}+𝛿𝑋(0),𝑉(0)=𝑍_{𝑉(𝑠𝑇)}+𝛿𝑉(0) for −𝜌3𝜎_𝑥,𝑖≤𝛿𝑋_𝑖≤𝜌3𝜎_𝑥,𝑖,−𝜌3𝜎_𝑣,𝑖≤𝛿𝑉_𝑖≤𝜌3𝜎_𝑣,𝑖, (6.5) where 𝑖=1÷3.

The components of a priori initial state vector correspond to a satellite on a circular orbit: 𝑋₁(0)=26×106m,𝑋₂(0)=𝑋₃𝑉(0)=0,₁(0)=0,𝑉₂(0)=3.685×103m/sec,𝑉₃(0)=0. (6.6)

The viewing system is defined by the relations 𝑡=0∶𝜆(0)=0,𝜙=0,ℎ=0(𝑦₁(0)=𝑦₂(0)=0).

The components of the actual initial state vector are determined by the inequalities (6.4) for 𝜌=0.2, |𝑋(0)|=26×106,|𝑉(0)|=3.685×103.

In the process of simulation, the recurrent algorithm was operating under the conditions 𝜏=10 seconds, 𝑘=20,𝑠=6. These conditions mean that each 20 pairs of viewing angles 𝑦1,𝑦2 observed during the sequence of time-intervals of length 𝑇=200seconds are replaced by pairs 𝑦_int,1,𝑦_int,2 of integral observation results. Iteration vectors 1,2,3,… form the output of the algorithm at the instants 6𝑇,12𝑇,18𝑇,….

In the process of simulation, the method of polynomial approximation for the components of the estimate vector was used for 𝑑=1 and 𝑑=2. For 𝑑=1, on each step of the recurrent algorithm, the components of the estimate vector are linear combinations of the last 12 components of the basic vector and consist of 12 integral observation results. For 𝑑=2, on each step of the recurrent algorithm, the components of the estimate vectors are linear combinations of 90 terms, which coincide with the last 90 components of the basic vector and consist of products of powers 0,1,2 of 12 integral observation results.

Below, we give the values of 𝜎_𝑥,𝑖 in meters and 𝜎_𝑣,𝑖 in meters per second attained in the process of several iterations.

A priori MSDs 𝜎^𝑎_𝑥,𝑖,𝜎^𝑎_𝑣,𝑖 of the components of the state vector of the space object at the instant 6𝑇 have the form 𝜎^𝑎_𝑥,1𝜎^𝑎_𝑥,2𝜎^𝑎_𝑥,3𝜎^𝑎_𝑣,1𝜎^𝑎_𝑣,2𝜎^𝑎_𝑣,33.76×10⁶3.73×10⁶3.67×10⁶5.45×10²5.39×10²5.23×10² (6.7) Iteration 1, 𝑑=2: 𝜎_𝑥,1𝜎_𝑥,2𝜎_𝑥,3𝜎_𝑣,1𝜎_𝑣,2𝜎_𝑣,35.00×10⁵1.60×10⁵7.89×10⁴3.44×106.55×101.62×10 (6.8) Iteration 2, 𝑑=2: 𝜎_𝑥,1𝜎_𝑥,2𝜎_𝑥,3𝜎_𝑣,1𝜎_𝑣,2𝜎_𝑣,36.57×10⁴3.18×10⁴9.84×10¹6.26×10⁻¹8.348.94×10⁻³ (6.9) Iteration 3, 𝑑=2: 𝜎_𝑥,1𝜎_𝑥,2𝜎_𝑥,3𝜎_𝑣,1𝜎_𝑣,2𝜎_𝑣,32.03×10²1.39×10²1.65×10⁻³8.71×10⁻³3.25×10⁻²1.53×10⁻⁷ (6.10) Iteration 4, 𝑑=1: 𝜎_𝑥,1𝜎_𝑥,2𝜎_𝑥,3𝜎_𝑣,1𝜎_𝑣,2𝜎_𝑣,35.804.228.78×10⁻⁹4.43×10⁻⁴9.40×10⁻⁴3.88×10⁻¹²1 (6.11) The above simulation data show that each iteration decreases, about 100 times, the MSD of the dispersal of the components of the state vector of the space object. By the end of iteration 4 (𝑡=24𝑇), its coordinates and the components of the velocity vector are practically precise.

7. Identification of Parameters of a Non Linear Servomechanism

We consider a task of identification of parameters of the nonlinear servo-mechanism, which consists sequentially of an inertial link from a constant of time 𝑇, an amplifier with factor 𝐾₁, an integrating link 1, a nonlinear link a gap 1 (magnitude of the gap is 𝑐₁), an amplifier with factor 𝐾₂, an integrating link 2, a nonlinear link a gap 2 (magnitude of the gap is 𝑐₂), a transmitter of a feedback, which exit a variable 𝑥₅. A variable 𝑈 is an input of the servo-mechanism which should be repeated by an exit of the servo-mechanism—a variable 𝑥₄.

The block diagram of the servo-mechanism is similar to the hydraulic booster which is used for turn steering surfaces of an aircraft [17].

The equations of driving of the servo-mechanism look like 𝑑𝑥₁1𝑑𝑡=(𝑇)(−𝑥1+Δ),𝑑𝑥₂𝑑𝑡=𝐾₁𝑥₁,𝑥₃=𝐹(𝑥₁,𝑥₂,𝑐₁),𝑑𝑥₄𝑑𝑡=𝐾₂𝑥₃,𝑥₅=𝐹(𝑥₃,𝑥₄,𝑐₂),Δ=𝑈−𝑥₁−𝑥₂−𝛼𝑥₅, (7.1) where 𝐹(𝑎,𝑏,𝑐)=𝑏−𝑐, if 𝑎>0, 𝐹(𝑎,𝑏,𝑐)=𝑏+𝑐, if 𝑎<0, and 𝛼=0, if the system is open loop, 𝛼=1 if the system is closed loop.

We suppose that 𝑈(𝑡)=sin(2𝜋𝑡)—to the servo-mechanism set harmonic motions with frequency 1 hertz and amplitude 1, a priori data about gaps: 0<𝑐₁,𝑐₂<0.1, remaining a priori data: 𝑇=𝑇(𝑛)(1+𝜃₃),𝐾₁=𝐾₁(𝑛)(1+𝜃₄),𝐾₂=𝐾₂(𝑛)(1+𝜃₅), where 𝑇(𝑛),𝐾₁(𝑛),𝐾₂(𝑛) are nominal (a priori) magnitudes, 𝜃_𝑖=𝜌𝜈_𝑖,0<𝜌<1,−1<𝜈_𝑖<1,𝑖=3,…,5. Further we suppose 𝜌=0.2 (we know parameters with errors of not surpassing 20 percents).

Outcomes of observations are variables 𝑥₅ measure with frequency 1000 hertz: 𝑦₁=𝑥₅(𝑡₁)+𝜀₁,…,𝑦_𝑘=𝑥₅(𝑡_𝑘)+𝜀_𝑘,…,𝑦_𝑁=𝑥₅(𝑡_𝑁)+𝜀_𝑛 They are used for identification of parameters. Random errors of observations 𝜀_𝑘 are uniformly distributed on a segment −[𝜀₀,𝜀₀]. The sequence of observations decomposes on segments in length 𝐿. The sums of outcomes of observations, belonging each segment, are sequence of inputs in MPA algorithm, which estimates 5 parameters. Parameters 𝑐₁,𝑐₂,𝑇,𝐾₁,𝐾₂ have numbers 1,…,5.

The exactitude of identification is presented by magnitudes 𝜎(𝑖)/𝜎_𝑎(𝑖),𝑖=1,…,5, and the rationes are certain in Section 4.

At 𝑇(𝑛)=0.01seconds,𝐾₁(𝑛)=60,𝐾₂(𝑛)=30,𝑑=1,𝐿=20,𝛼=0,𝜀₀=0.0001 and various times of observations 𝑡_𝑁 the mentioned rationes are presented as follows: 𝜎(𝑖)/𝜎_𝑎𝑡(𝑖)12345_𝑁𝑡=2sec0.110730.000300.142570.129810.00091_𝑁=4sec0.056320.000190.078600.067380.00012 (7.2)

8. Conclusion

The algorithm described above gives a sequence of convergent estimates for the vector of conditional expectation and the corresponding estimate error covariance matrices practically in any situation, provided that one knows the mathematical models of the system and observations necessary for leaning of the algorithm on the basis of the Monte-Carlo method.

It is convenient to compress data by replacing the high-dimensional vector of primary observations with a vector of a substantially smaller dimension, whose components are functions of primary observations. Such a replacement is likely to cause an insignificant increase of theoretical estimate error variance.

Thus, in Sections 5 and 6 data compression was implemented by passing from a 200-dimensional vector of primary observations to a 12-dimensional vector whose components are formed by sums of 20 consecutive components of the primary vector.

The development of rational methods for the compression of primary data is necessary for practical utilization of the proposed general algorithm, which is a unique example of a situation in which the Mote-Carlo method provides synthesis of a new efficient method of nonlinear estimation.

Acknowledgment

This work has been supported by the Russian Foundation for Basic Research.