# A Bayes Estimator of Parameters of Nonlinear Dynamic Systems

Authors: I. A. Boguslavsky
Journal: Mathematical Problems in Engineering
 (2009)

---

## 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 Y. 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 d. The series is a vector polynomial approximation of the vector  E(θ∣Y), 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 d increases. These coefficients are found by the Monte-Carlo method and a process
of recurrent calculations that do not require
matrix inversion.


---

## Body

## 1. Introduction

Consider a dynamical system whose mathematical model has the form(1.1)x(k+1)=f(x(k),ηk,k).
Here, the integers 0,1,…,k,… refer to time instants t0,t1,…,tk,…; x(k)∈Rn is the state vector of the dynamical system at the instant tk;  ηk∈Rn is the vector of random perturbations; f(x(k),ηk,k) is a given function (nonlinearly depending on its arguments, in general);  x(0),ηk are random vectors with given distributions.Let the results of observations at the instantst1,…,tk,…,tN be described by a scalar sequence y1,…,yk,…,yN generated by the mathematical model
(1.2)yk=Hk(x(k),ξk),
where H(k,x(k)) is a given function of its arguments, and the distribution is given for the sequence of random variables ξ1,…,ξk,…,ξN.In what follows, it is assumed thaty1,…,yk are components of the vector Yk.Aparameter 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 θ=x(0); the problem of filtration, if θ=x(N); the problem of extrapolation, if θ=x(N*),N*>N.It is well known that the mean-square optimal estimate of the random vectorθ on the basis of observations of the vector YN is the vector of conditional expectation E(θ∣YN). Therefore, we consider the problem of creating the MPA algorithm for the construction of approximations converging to the vector E(θ∣YN).There is only one requirement with regard to computer power and the expressions (1.1) and (1.2): for given distributions of random variables x(0),ηk,ξi, 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 sequencey1,…,yk,…,yN are sent to the input of the algorithm at the instants t1,…,tk,…,tN.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 η(k)≡0, the algorithm determines the estimate vector θ̂=zx(0)(YN) that approximately minimizes the heuristic quality functional S(YN,zx(0)(YN)) with the constraints (1.1). The estimation quality functional is the sum of squared residuals
(1.3)S(YN,zx(0)(YN))=∑i=1i=N(yi-H(i,xi(zx(0)(YN))))2.
Here, the vectors xi(zx(0)(YN)) are obtained by consecutive calculations of state vectors on the basis of (1.1) with the initial data zx(0)(YN). 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 zx0(YN).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 S(YN,zx0(YN)) 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 instanttk-1, approximations have been found for the first and the second statistical moments of the state vector at the instant tk-1 (the estimation vector for this vector coincides with the first statistical moment). On the time interval [tk-1,tk-0] (prior to observation results at the instant tk) prediction is made and its results are approximations for the first and the second statistical moments of the state vector at the instant tk-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 tk-1 to tk.Step 2.
After actual observations at the instanttk 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 tk 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 tk.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 r×1 to be estimated upon fixing the vector YN.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 givenYN, 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 YN. For every term of series sum of powers bounded above given integer d. For given YN, the approximations uniformly bounded converge to the vectors posterior expectation of the estimated parameters vector by d→∞;(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 vectorYN. 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 YN;(iii)
there is a sequencer vectors such as each consist from k leading to a component of vector θ, k=1,…,r. After adjustment by given d,N, 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 thatd is a given positive integer number and the set of integer numbers a1,…,aN consists of all nonnegative solutions of the integer inequality a1+⋯+aN≤d, the number of which we denote by m(d,N). The value m(d,N) is given by the recurrent formula proved by induction,
(2.1)m(d,N)=m(d-1,N)+(N+d-1)⋯kd!,m(1,N)=N.
We obtain the vector WN(d) of dimension m(d,N)×1, the components w1,…,wm(d,N) of which are all possible values of the form that represent the powers of measurable values.Next, we define a basic vectorV(d,N) of dimension (r+m(d,N))×1,V(d,N)=∥θWN(d)∥T.Step 2.
We use a known statistical generator of random vectorsx(0),ηk,ξi to solve repeatedly the Cauchy problem for (1.1) for a given initial conditions x(0), a control law u and various realizations of random vectors x(0),ηk,ξi and k. As a result, the computer memory will contain a set of realizations of random basic vectors V(d,N) sufficient for the calculation of the statistical characteristics of the basic vector V(d,N). We apply the Monte-Carlo method to find the prior first and second statistical moments of the vector V(d,N), that is, the mathematical expectation V¯(d,N)=E(V(d,N)), and the covariance matrix
(2.2)CV(d,N)=E((V(d,N)-E(V(d,N)))(V(d,N)-E(V(d,N)))T).Implementation of  Step2 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  andN and a fixed vector YN, we assign the vector ZE(θ∣YN)(WN(d)) to be the solution to the estimation problem. This vector gives an approximate estimate of the vector E(θ∣YN) that is optimal in the root-mean-square sense on the set of vector linear combinations of components of the vector WN(d). Then the vector ZE(θ∣YN)(WN(d)) is an element of this set and is an approximate mean-square optimal estimate for the vector E(θ∣YN) of the conditional expectation of the vector θ.
The vectorZE(θ∣YN)(WN(d)) can be presented in a following aspect:
(2.3)ZE(θ∣YN)(WN(d))=∑a1+⋯+aN≤dλ(a1,…,aN)y1a1⋯yNaN,
where vectors λ(a1,…,aN) are some weight vectorial factors.The vectorV¯(d,N)=E(V(d,N)) and the matrix CV(d,N) are the initial conditions for the process of recurrent calculations that realizes the principle of observation decomposition [6] and consists of m(d,N) steps. Once the final step is performed, we obtain vector weight coefficients λ(a1,…,aN) for (1.1). Moreover, we determine the matrix C(d,N), 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 C(d,N), 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 WN(d) is close to linear combination of its previous components. To implement the recursion, we process the components of the vector WN(d) one after another. However, the adjustment of the algorithm performed by applying the Monte-Carlo method to find the vector V¯(d,N)=E(V(d,N)) and the matrix CV(d,N) takes into account a priori ideas on stochastic structure of components of the whole set of possible vectors WN(d) that can appear in any realizations of the random vectors x(0),ηk,ξi and k 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 d the error estimate vector ZE(θ∣YN)(WN(d))-E(θ∣YN) 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 integerd and the increase of the number of realizations, if the Monte-Carlo method is used.If the numberm(d,N) of the approximating polynomial series is not small, a large volume of calculations has to be performed (after the initial choice of the integers n,N,d) for the determination of the vector coefficients λ(a1,…,aN).


## 3. Mean-Square Optimal Estimate for the Conditional Expectation VectorE(θ∣YN)

Consider the principal algorithm (PA) for solving the problem of the mean-square optimal estimate for the vectorE(θ∣YN). It is known that the vector E(θ∣YN) is the mean-square optimal estimate for the vector θ after the vector YN 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 vectorE(θ∣YN). To that end, we obtain an estimate for the vector E(θ∣YN) which is linear with respect to the components of the vector WN(d) and is mean-square optimal.In what follows the vector of that estimate is denoted byZE(θ∣YN)(WN(d)).An explicit expression for the estimate vector is obtained after the calculation of the elements of the vectorV¯(d,N)=E(V(d,N)) and the covariance matrix CV(d,N), which coincide with the first and the second (centered) statistical moments for the vector V(d,N). This vector and this matrix can be divided into blocks whose structure can be represented as follows:
(3.1)E(E(θ∣YN))=E(θ),E(E(W(d,N)∣YN))=E(W(d,N)),Cθ(0)=E((E(θ∣YN)-E(θ))(E(θ∣YN)-E(θ))T)=E((θ-E(θ))(θ-E(θ))T),LN(d)=E((E(θYN)-E(E(θ∣YN(d)))(WN(d)-E(WN(d)))T)=E(θ)WN(d)T)-E(θ)E(WN(d))T,QN(d)=E((WN(d)-E(WN(d)))(WN(d)-E(WN(d)))T).
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 E(θ∣YN).Let(3.2)ZE(θ∣YN)(WN(d))=E(θ)+ΛN(d)(WN(d)-E(WN(d))),
where ΛN(d) is an n×m(d,N)-matrix satisfying the equation
(3.3)ΛN(d)QN(d)=LN(d).Let(3.4)Z̃θ(WN(d))=z+Λ̃N(d)(WN(d)-E(WN(d))),
where z is an arbitrary n×1-vector and Λ̃N(d) is an arbitrary n×m(d,N)-matrix. Let CZ(d,N) and C̃(d,N) be error estimate covariance matrices for the vector E(θ∣YN) in terms of the estimate vectors ZE(θ∣YN)(WN(d)) and Z̃θ(WN(d)).Lemma 3.1.
A matrixC̃(d,N)-CZ(d,N) is the nonnegative definite matrix: CZ(d,N)≤C̃(d,N).Lemma3.1 follows from the identity
(3.5)C̃(d,N)=CZ(d,N)+(ΛN(d)-Λ̃N(d))(ΛN(d)-Λ̃N(d))T+(ΛN(d)QN(d)-LN(d))(Λ̃N(d)-ΛN(d))T+((ΛN(d)QN(d)-LN(d))(Λ̃N(d)-ΛN(d))T)T+(z-E(θ)(z-Eθ)T).Corollary from Lemma3.1
The vectorZE(θ|YN)(WN(d)) is a mean-square optimal estimate for the vector E(θ∣YN) on the set of estimates linear with respect to the components of the vector WN(d).IfQN(d)>0, then the estimate vector is unique and
(3.6)ZE(θ∣YN)(WN(d))=E(θ)+LN(d)QN(d)-1(WN(d)-E(WN(d))).The covariance matrixCZ(d,N) for the error of the estimate of vector E(θ∣YN) is defined by
(3.7)CZ(d,N)=Cθ(0)-ΛN(d)LN(d).IfQN(d)≥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 λ(a1,…,aN) 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 WN(d).Let us examine asymptotic errors of the estimates obtained when using formula (3.2).For a given vectorYN, suppose that the vector E(θ∣YN) is defined as a function of the argument YN in some a priori given compact domain ΩYN⊂RN and is continuous in that domain. Then the following result holds.Theorem 3.2.
(3.8)supYN∈ΩYN|ZE(θ∣YN)(WN(d))-E(θ∣YN)|⇒0,d⇒∞.Proof.
According to the multidimensional version of the Weierstrass theorem [7], for any ε>0 there is a multidimensional polynomial P(WN(dε)) such that
(3.9)supYN∈ΩYN|P(WN(dε))-E(θ∣YN)|<ε.
This relation can be rewritten as
(3.10)supYN∈ΩYN|P(WN(d))-E(θ∣YN)|⇒0,d⇒∞.
Let C be the covariance matrix for the random vector P(WN(d))-E(θ∣YN):
(3.11)C=E(P(WN(d))-E(θ∣YN))(P(WN(d))-E(θ∣YN))T.
From (3.10), it follows that
(3.12)C⇒0n,d⇒∞.
By construction, the vector ZE(θ∣YN)(WN(d)) gives a mean-square optimal estimate for the vector θ and this estimate is linear with respect to the components of the vector WN(d). But Lemma 3.1 implies that for any other nonoptimal linear estimate, in particular that like the vector P(WN(d)), we have C≥CZ(d,N). Hence, taking into account (3.12), we get
(3.13)CZ(d,N)⇒0n,d⇒∞.
The statement (3.13) is equivalent to (3.8), if we take into account that
(3.14)CZ(d,N)=∫(ZE(θ∣YN)(WN(d))-E(θ∣YN))×(ZE(θ∣YN)(WN(d))-E(θ∣YN))Tp(θ,YN)dθdYN,
where p(θ,YN) is the joint distribution density for the random vectors θ,YN. 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 m(d,N), 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 vectorYN and an integer d we have constructed the vector WN(d) with the components w1,w2,…,wm(d,N).The implementation of the algorithm starts with approximate calculation (by the Monte-Carlo method) ofn+m(d,n)-dimensional integrals that determine the components of the vector V¯(d,N)=E(V(d,N)) and the matrix CV(d,N)=E((V(d,N)-V¯(d,N))(V(d,N)-V¯(d,N))T).  Further, splitting this vector and this matrix into suitable blocks and using (2.3), we find the desired estimate vector ZE(θ∣YNk)(WN(d)) and the estimate error covariation matrix CZ(d,N).However, this approach is unreasonable, since it requires the inversion of the matrixQN(d), which is a difficult task for QN(d) of a large dimension m(d,N)×m(d,N) or QN(d) 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 ofm(d,N) steps of calculations of the first and the second statistical moments for a sequence of special vectors V1,…,Vi,…,Vm(d,N) after a priori moments V¯(d,N),CV(d,N) for the basic vector V(d,N) have been found.Denote byV1 the vector formed by the components of the basic vector V(d,N) remaining after the elimination of the component w1;…; by Vi+1 the vector formed by the components of Vi remaining after the elimination of wi; and so forth.The quantitywm(d,N) is the last component of the vector WN(d), and after its elimination, the last vector Vm(d,N) turns out equal to the estimate vector ZE(θ∣YN)(WN(d)).On Step1, formulas (3.2) and (3.6) are used for the calculation of the following quantities:(i)
the vectorV¯1 that yields an estimate for the vector V1 which is mean-square optimal and linear with respect to w1;(ii)
the estimate error covariance matrixC(V1).The estimate vector consists of the conditional expectation estimate vectorE(θ∣w1) and a vector of dimension (m(d,N)-1)×1. For w1 fixed, the latter coincides with the vector of statistical prediction of mean future  values of the quantities w2,…,wm(d,N).Note that the calculations of  Step1 are based on a priori values V¯(d,N),C(V(d,N)) found beforehand.In a similar way, suppose that on steps1,…,i of the calculation process, the vector V¯i and the matrix C(Vi) have been found after fixing the quantities w1,…,wi.On stepi+1, the following quantities are calculated by means of (3.2) and (3.6):(i)
the vectorV¯(i+1) of the estimate of Vi+1, which is mean-square optimal and linear with respect to w1,…,wi+1;(ii)
the estimate error covariance matrixC(V(i+1)).The vectorV¯(i+1) is composed of the conditional expectation estimate E(θ∣w1,…,wi+1) and a vector which, after the quantities w1,…,wi+1 have been fixed, is the vector of statistical prediction for the mean values of “future” quantities wi+2,…,wm(d,N).Note that the calculations of stepi+1 are based on the already determined quantities V¯i,C(Vi), which it is natural to call a priori first and second statistical moments  of the “future’’ quantities wi+1,…,wm(d,N). The vector V¯i and the matrix C(Vi) represent a priori data regarding the statistical moments of the components of the vector Vi+1 before the quantity wi+1 is sent to the input of the algorithm.The formulas of the computation precess described above are the following:(4.1)V¯(i+1)=V¯i1+qi-1bi(wi+1-zwi+1),C(Vi+1)=C(Vi)1-qi-1bibiT.
Here, the scalar zwi+1 is the (n+1)th component of the vector V¯i (this scalar is a linear mean-square optimal estimate for the component wi+1 after the algorithm has used the components w1,…,wi); the vector V¯i1 is obtained from the vector V¯i by eliminating its component zwi+1; the scalar qi is the (n+1)th diagonal element of the matrix C(Vi) (this scalar is the variance of the estimate error for the component wi+1 after the components w1,…,wi have been used); the matrix C(Vi)1 is obtained from the matrix C(Vi) by elimination of its (n+1)th row and its (n+1)th column; the vector bi coincides with the (n+1)th column of the matrix C(Vi) with the (n+1)th component excluded.If the scalarqi happens to be close to zero, then the component wi+1 becomes close to a linear combination of the components w1,…,wi. In this case, the component wi+1 gives no new information about θ and should be dropped from the computation process..Note that the sequence of random variables(wi+1-zwi+1) is a renewable sequence.The top left block of the(n×n)-matrix C(Vi) coincides with the estimate error covariance matrix CZ(d,i) for the vector E(θ∣w1,…,wi) after the vector Wi(d) has been used by the algorithm.Denote byl(i) the vector formed by the first n components of the vector bi. The formula representing the evolution of the covariance matrix CZ(d,i) as a function of i (the number of the components of the vector Wi(d)) has the form
(4.2)CZ(d,i)=Cθ(0)-q1-1l(1)l(1)T-⋯-qi-1-1l(i-1)l(i-1)T.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θ1,…,θr 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σ(i,cal) is determined theoretically by the calculation of variances, that is, the diagonal elements of the covariance matrix CZ(d,N).The quantityσ(i,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 axesX, Y, Z and their origin O. The corresponding projections ωx,ωy,ωz of the angular velocity of the body, ω, are measured by angular velocity sensors. The projections Mx,My,Mz of the moment M of external forces are assumed known.In the general situation, a solid has 6 moments of inertia relative to a rectangular reference frameX,Y,Z, among which three are called centrifugal moments. It is known that for a point O there is a rectangular coordinate system Xo,Yo,Zo 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, Jx,Jy,Jz, are called the principal moments of inertia relative to the point O.Letα,β,γ be the Euler angles that determine the angular position of Xo,Yo,Zo with respect to X,Y,Z, let A 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α,β,γ,Jx,Jy,Jz 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θ1=α,θ2=β,θ3=γ,θ4=Jx,θ5=Jy,θ6=Jz, provided that the solid is rotating in space and in the reference frame the quantities ωx,ωy,ωz,Mx,MY,Mz are observed as discrete-time functions.In the continuous-time model, the equations of motion of a nonlinear dynamical system have the form(5.1)Jxdωxodt+(Jz-Jy)ωyoωzo=Mxo,Jydωyodt+(Jx-Jz)ωzωx=Myo,Jzdωzodt+(Jy-Jx)ωxoωyo=Mzo,(5.2)dJxdt=0;dJydt=0;dJzdt=0,dαdt=0;dβdt=0;dγdt=0.
Here, ωxo,ωyo,ωzo,Mo,Myo,Mzo are projections of the vectors ω and M to the axes Xo,Yo,Zo:
(5.3)ωo=Aω,Mo=AM.
Further, assume that observations of the time-functions ωx,ωy,ωz,Mx,MY,Mz 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 ωx(0)=1/20,ωy(0)=-1/20,ωz(0)=1/20andMx=100kg,MY=-100kg,Mz=100kg.Settingyx(tk)=ωx(tk),yy(tk)=ωy(tk),yz(tk)=ωz(tk), we find that at the instant tk the components of the initial observation vector y(tk) are related to nonobservable quantities α,β,γ and nonobservable functions ωxo,ωyo,ωzo by
(5.4)y(tk)=ATωo(tk).Using the sequenceyx(tk),yy(tk),yz(tk) of primary observation vectors and (5.1), it is required to estimate the parameters Jx,Jy,Jz,α,β,γ.The secondary observationsyx(i),yy(i),yz(i),i=1÷4, which serve as the input data of the estimation algorithm, are found by summing up the components of each l vectors of consecutive primary observations:
(5.5)yx(i)=∑j=1j=lyx(tl(i-1)+j),yy(i)=∑j=1j=lyy(tl(i-1)+j),yz(i)=∑j=1j=lyz(tl(i-1)+j).
In our model, it is assumed that l=25 and the array of 300 primary observations has been converted to an array of N=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α,β,γ,Jx,Jy,Jz and yx(i),yy(i),yz(i),i=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:
(5.6)-ρ≤α,β,γ≤ρ,Jxo(1-ρ)≤Jx≤Jxo(1+ρ),Jyo(1-ρ)≤Jy≤Jyo(1+ρ),Jzo(1-ρ)≤Jz≤Jzo(1+ρ),
where Jxo=4000,Jyo=20000,andJzo=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σa(i) be a priori mean-square deviation (MSD) of the scattering of the unknown parameter θi, and let σ(i) 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 θi is characterized by the ratios σ(i)/σa(i),i=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 d=2 (the estimate of θi is realized by second-order polynomials with respect to 12 secondary scalar observations; the number of terms of the polynomials is equal to m(2,12)=90):ρ=0.2(5.7)σ(1)/σa(1)σ(2)/σa(2)σ(3)/σa(3)σ(4)/σa(4)σ(5)/σa(5)σ(6)/σa(6)0.55680.04400.56140.10050.08620.0886ρ=0.1(5.8)σ(1)/σa(1)σ(2)/σa(2)σ(3)/σa(3)σ(4)/σa(4)σ(5)/σa(5)σ(6)/σa(6)0.50690.00630.50790.01230.01640.0161ρ=0.05(5.9)σ(1)/σa(1)σ(2)/σa(2)σ(3)/σa(3)σ(4)/σa(4)σ(5)/σa(5)σ(6)/σa(6)0.49830.00120.49850.00170.00280.0026
These data show that for d=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 β,Jx,Jy,Jz, 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 d, 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 systemWindow [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)dXidt=Vi,dVidt=-μXiR.
Here, R=(X12+X22+X32)1/2, μ is the gravity constant, Xi,i=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 instantt the center of viewing is located at a point of the Earth surface with geographical coordinates λ=ωt (longitude), φ (latitude), h (altitude), where ω is the angular speed of the Earth rotation. In the rotating right geocentric coordinate system xi,i=1÷3, the coordinates of the center of viewing are X1′=recos(φ)-hsin(φ),X2′=0,X3′=resin(φ)+hcos(φ), where re 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 angley1 and, in the local vertical plane, by an angle y2. These angles are defined by
(6.2)y1=arctan(ba),y2=arctan(c(a2+b2)1/2),
where
(6.3)a=cos(λ)cos(φ)X1+sin(λ)cos(φ)X2+sin(φ)X3-re,b=-sin(λ)X1+cos(λ)X2,c=-sin(φ)cos(λ)X1+sin(φ)sin(λ)X2+sin(φ)X3-h.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 Xi,i=1÷3.Letτ be the time interval between two consecutive observations of pairs of viewing angles of the space object. After the first k observations, we construct integral observation results yint,1,yint,2. Here, yint,1 is the sum of the first k fixed quantities y1; similarly, yint,2 is the sum of the first k fixed quantities y2. In like manner, we construct integral pairs yint,3,yint,4,…,yint,2s-1,yint,2s,…, and so forth. The sequence of these random variables is a sequence of inputs to the recurrent algorithm at the instants T,2T,…,rT,…, where T=kτ. The summation operator reduces the effect of random tracking errors on the accuracy of the estimates.Choosing an integers, we introduce the following definitions. Let us divide the sequence of pairs of integral observation results into a sequence of segments, each containing s 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 instantsT after sk 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 instant2sT after 2sk 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 powers2s of integral observation results; the sum of these powers does not exceed the given integer d. For d=1, the dimension of the basic vector 1 is equal to 6+2s; for d=2, the dimension of the basic vector 1 is equal to 6+2s+(2s+1)s, and so forth.A priori scattering of the vectorsX(0),V(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 X(0)=X(0)N+δX(0),V(0)=V(0)N+δV(0), where X(0)N is the nominal initial radius-vector of the observed space object and V(0)N is the nominal initial vector of its velocity. It is assumed that V(0)N is the velocity vector of an object moving on a circular orbit of radius |X(0)N|. Uniform scattering of the components of random vectors δX(0),δV(0) is determined by the inequalities
(6.4)-ρ|X(0)|≤δX1,δX2,δX3≤ρ|X(0)|,-ρ|V(0)|≤δV1,δV2,δV3≤ρ|V(0)|,
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ΩX,ΩV, which are calculated by the Monte-Carlo method.After a priori statistical characteristics of the basic vector 1 have been defined,2s random variables are sent to the input of the recurrent algorithm, namely, 2s integral observation results for the viewing angles on the time interval sT. After 2s steps of calculations, the recurrent algorithm determines the iteration vector 1, which is the estimate vector ZX(sT),ZV(sT) for the state vector of the space object at the instant sT.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. SetX(0)=ZX(sT)+δX(0),V(0)=ZV(sT)+δV(0) for
(6.5)-ρ3σx,i≤δXi≤ρ3σx,i,-ρ3σv,i≤δVi≤ρ3σv,i,
where i=1÷3.The components of a priori initial state vector correspond to a satellite on a circular orbit:(6.6)X1(0)=26×106m,X2(0)=X3(0)=0,V1(0)=0,V2(0)=3.685×103m/sec,V3(0)=0.The viewing system is defined by the relationst=0:λ(0)=0,ϕ=0,h=0(y1(0)=y2(0)=0).The components of the actual initial state vector are determined by the inequalities (6.4) for ρ=0.2, |X(0)|=26×106,|V(0)|=3.685×103.In the process of simulation, the recurrent algorithm was operating under the conditionsτ=10 seconds, k=20,s=6. These conditions mean that each 20 pairs of viewing angles y1,y2 observed during the sequence of time-intervals of length T=200seconds are replaced by pairs yint,1,yint,2 of integral observation results. Iteration vectors 1,2,3,… form the output of the algorithm at the instants 6T,12T,18T,….In the process of simulation, the method of polynomial approximation for the components of the estimate vector was used ford=1 and d=2. For d=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 d=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σx,i in meters and σv,i in meters per second attained in the process of several iterations.A priori MSDsσx,ia,σv,ia of the components of the state vector of the space object at the instant 6T have the form
(6.7)σx,1aσx,2aσx,3aσv,1aσv,2aσv,3a3.76×1063.73×1063.67×1065.45×1025.39×1025.23×102
Iteration 1, d=2:
(6.8)σx,1σx,2σx,3σv,1σv,2σv,35.00×1051.60×1057.89×1043.44×106.55×101.62×10
Iteration 2, d=2:
(6.9)σx,1σx,2σx,3σv,1σv,2σv,36.57×1043.18×1049.84×1016.26×10-18.348.94×10-3
Iteration 3, d=2:
(6.10)σx,1σx,2σx,3σv,1σv,2σv,32.03×1021.39×1021.65×10-38.71×10-33.25×10-21.53×10-7
Iteration 4, d=1:
(6.11)σx,1σx,2σx,3σv,1σv,2σv,35.804.228.78×10-94.43×10-49.40×10-43.88×10-121
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 (t=24T), 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 timeT, an amplifier with factor K1, an integrating link 1, a nonlinear link a gap 1  (magnitude of the gap is c1), an amplifier with factor K2, an integrating link 2, a nonlinear link a gap 2  (magnitude of the gap is c2), a transmitter of a feedback, which exit a variable x5. A variable U is an input of the servo-mechanism which should be repeated by an exit of the servo-mechanism—a variable x4.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(7.1)dx1dt=(1T)(-x1+Δ),dx2dt=K1x1,x3=F(x1,x2,c1),dx4dt=K2x3,x5=F(x3,x4,c2),Δ=U-x1-x2-αx5,
where F(a,b,c)=b-c, if a>0, F(a,b,c)=b+c, if a<0, and   α=0, if the system is open loop, α=1 if the system is closed loop.We suppose thatU(t)=sin(2πt)—to the servo-mechanism set harmonic motions with frequency 1 hertz and amplitude 1, a priori data about gaps: 0<c1,c2<0.1, remaining a priori data: T=T(n)(1+θ3),K1=K1(n)(1+θ4),K2=K2(n)(1+θ5), where T(n),K1(n),K2(n) are nominal (a priori) magnitudes, θi=ρνi,0<ρ<1,-1<νi<1,i=3,…,5. Further we suppose ρ=0.2 (we know parameters with errors of not surpassing 20 percents).Outcomes of observations are variablesx5 measure with frequency 1000 hertz: y1=x5(t1)+ε1,…,yk=x5(tk)+εk,…,yN=x5(tN)+εn They are used for identification of parameters. Random errors of observations εk are uniformly distributed on a segment -[ε0,ε0]. The sequence of observations decomposes on segments in length L. The sums of outcomes of observations, belonging each segment, are sequence of inputs in MPA algorithm, which estimates 5 parameters. Parameters c1,c2,T,K1,K2 have numbers 1,…,5.The exactitude of identification is presented by magnitudesσ(i)/σa(i),i=1,…,5, and the rationes are certain in Section 4.AtT(n)=0.01seconds,K1(n)=60,K2(n)=30,d=1,L=20,α=0,ε0=0.0001 and various times of observations tN the mentioned rationes are presented as follows:
(7.2)σ(i)/σa(i)12345tN=2sec0.110730.000300.142570.129810.00091tN=4sec0.056320.000190.078600.067380.00012


## 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 Sections5 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.



---

*Source: 801475-2009-07-19.xml*