An Algorithm for Parameters Identification of an Aircraft’s Dynamics

Development of efficient parameter identification methods for the model of a dynamic system based on real-time measurements of some components of its state vector should be taken as one of the most important problems of applied statistics and computational mathematics. Calculating themotion of the system given the initial conditions and its mathematical model is conventionally called the direct problem of dynamics. Then, the inverse problem of dynamics would be the problem of identifying the system model parameters based on measurements of certain components of the state vector provided that the general structural scheme of the model is known from physical considerations. Such an inverse problem corresponds to identification problem for the dynamic system representing an aircraft. In this case, the general structural scheme of the model (motion equations) follows from the fundamental laws of aerodynamics.


Introduction
Development of efficient parameter identification methods for the model of a dynamic system based on real-time measurements of some components of its state vector should be taken as one of the most important problems of applied statistics and computational mathematics. Calculating the motion of the system given the initial conditions and its mathematical model is conventionally called the direct problem of dynamics. Then, the inverse problem of dynamics would be the problem of identifying the system model parameters based on measurements of certain components of the state vector provided that the general structural scheme of the model is known from physical considerations. Such an inverse problem corresponds to identification problem for the dynamic system representing an aircraft. In this case, the general structural scheme of the model (motion equations) follows from the fundamental laws of aerodynamics.
In many cases, modern computational methods and wind tunnel experiments can provide sufficient data on nominal parameters of the mathematical model -nominal aerodynamic characteristics of the aircraft. Nevertheless, there exist problems [1] that require correcting nominal parameters based on measurements taken in real flights. These imply (1) verifying and interpreting theoretical predictions and results of wind tunnel experiments (flight data can also be used to improve ground prediction methods), (2) obtaining more exact and complete mathematical models of the aircraft dynamics to be applied in designing stability enhancement methods and flight control systems, (3) designing flight simulators that require more accurate dynamic aircraft profile in all flight modes (many motions of aircrafts and flight conditions can be neither reconstructed in the wind tunnel nor calculated analytically up to sufficient accuracy or efficiency), (4) extending the range of flight modes for new aircrafts, which can include quantitative determination of stability and impact of control when the configuration is changed or when special flight conditions are realized, (5) testing whether the aircraft specification is compliant.
Furthermore, dimensionless numbers at the nodes of one-or two-dimensional tables found in wind tunnel experiments serve as nominal values in the aerodynamic parameter identification problem of the aircraft. This causes the vector that corrects these parameters determined by the algorithm processing digital data flows received from the aircraft sensors to have a significant dimension of the order about several tens or hundreds.
It is worth noting that the USA (NASA) is doing extensive work on theoretical and practical aircraft identification by test flights. In 2006 alone, in addition to many journal publications, American Institute of Aeronatics and Astronatics (AIAA) published three fundamental monographs [1-3] on the subject. An implementation of multiple NASA recommended algorithms for identification problems, SIDPAS (Systems Identification Programs for Aircraft) software package written in MATLAB M-files language is available on the Internet as an appendix to [1]. Various existing identification methods published in monographs on statistics and computational mathematics are widely reviewed in [1].
For the most general identification method, one should take the known nonlinear least squares method [4] that forms the sum of errors squared -differences between the real measurements and their calculated analogues obtained by numerical integration of motion equations of the system for some realization of the vector of unknown parameters.
Successful identification yields the vector of parameters that delivers the global minimum to the above mentioned sum of errors squared. Still, this criterion is statistically valid only for linear identification problems, in which measurements are linear with respect to the unknown vector of parameters.
Implementing the nonlinear least squares method to correct nominal parameters of the aircraft based on its test flight data involves computational challenges. These arise when the dimension of the correction vector is big and the sum of errors squared as the function of the correction vector has multiple relative minimums or when variations of the Newton's method are applied, with the sequence of local linearizations performed to find stationary points of this function. In [1], the regression method supported by lesq.m, smoo.m, derive.m, and xstep.m files in SIDPAS is recommended for practical applications.
Suppose the motion equations of the system and the sequence of measurements have the form where x(t k ) is the n × 1-dimensional vector of the system states at the current instant t and at the given instants t k , k = 1, ..., N, ϑ is the r × 1-vector of nominal (known) parameters of the system, η is the vector of unknown parameters that serves as the correction vector for the nominal vector ϑ after the results of measurements are stochastically processed, u is the control vector of the system, f (...) is the given vector-function, y k is the sequence of vectors-results of measurements, H k (...) is the given vector function, and ξ k , k = 1, ..., N is the sequence of random vectors-errors of measurements with the given random generator for the mathematical simulation.
We can state the identification problem for the vector η as follows. Find the estimate as the function of the vector Y N formed of the results of all measurements y 1 , ...yN.

The regression method given in [1] solves this problem under the following limitations
(1) all components of the state vector can be measured : y k = x(t k )+ξ k , (2) at the measurement instants t k , the algorithm constructs the estimate of the vector of derivatives dx/dt, (3) the vector function f (x, ϑ + η, u) linearly depends on the vector η.
These fundamental limitations of the regression method duplicate features of the identification algorithm from [5]. The substantial drawback of the algorithm [5] and the algorithm of the regression method is that they do not allow using the mathematical model to analyze theoretically (without applying the Monte-Carlo method) observability conditions of components of the vector of parameters to be identified for the preliminary given control law for the test flight of the aircraft and information on random errors of its sensors. Note that this is the drawback of all known numerical methods that solve nonlinear identification problems.
Relations (0.1) and (0.2) show that when conditions (1)-(3) are met and N is sufficiently big, the estimation vector satisfies the overdetermind system of linear algebraic equations, with methods to solve it being well known. The given conditions seem to be rather rigid and may be hard-to-implement. For instance, it is arguable whether one can construct the vector of derivatives dx/dt sufficiently accurately given the real turbulent atmosphere conditions, which imply that the outputs of the angle of attack and sideslip sensors inevitably include random and unpredictable frequency components.
All this justifies the development of new identification algorithms that can be applied to dynamic systems of a rather general class and do not possess drawbacks of NASA algorithms. The proposed multipolynomial approximation algorithm (MPA algorithm) serves as such a new identification algorithm.

Statement of the problem and basic scheme of the proposed identification algorithm
The general scheme for identifying aerodynamic characteristics of the aircraft by the test flight data is as follows [1]. Motion equations of the aircraft (0.1) and system (0.2) of measurements of motion characteristics of the aircraft are given. The vector ϑ is the vector of nominal aerodynamic parameters determined in the wind tunnel experiment. Calculated by the results of real (test) flight, the vector η is used to correct the vector ϑ.
When the aircraft flies, its computer fixes the digital array of initial conditions and time functions, viz. current control surface angles and measurements of some motion parameters of the aircraft (some components of the vector x(t) of the state of the aircraft) received from its sensors. Note that selecting the criterion for optimal or, at least, rational mode to control the test flight is a separate problem and lies beyond our further consideration. The current motion characteristics measured as the time function such as angles of attack and sideslip and components of the vector of angular velocity and g-load obtained by the inertial system of the aircraft are registered for real (not known for sure) aerodynamic parameters of the aircraft (parameters ϑ + η) and can be called measured characteristics of the perturbed motion.
Once the flight under the mentioned (given) initial conditions and time functions (control surface angles) is completed, nominal motion equations (equations of form (1) for η = 0) are integrated numerically for the nominal aerodynamic parameters of the aircraft. For the calculated characteristics of the nominal motion of the aircraft one should take the obtained data -components of the state vector of the aircraft as the function of discrete time. Differences between measurable characteristics of the perturbed motion and calculated characteristics of the nominal motion serve as carriers of data on the unknown vector η that shows the difference between real and nominal aerodynamic parameters.
The input of the MPA identification algorithm receives the vector of initial conditions and control surface angles as functions of time and arrays of characteristics of nominal and perturbed motions.
The output of the algorithm isη(Y N ) , which is the correction vector for nominal aerodynamic parameters.
The identification algorithm is efficient if the motion equations integrated numerically with the corrected aerodynamic parameters yield such motion characteristics ϑ +η(Y N ) (corrected characteristics, in what follows) that are close to real (measurable) characteristics.
In this work, we consider the technology of applying the Bayes MPA algorithm [6,7] to solve identification problems on the example of the aircraft, for which nominal aerodynamic parameters of the pitching motion are the nominal parameters of one of an "pseudo" F-16 aircraft.
We replace real flights by mathematical simulation, with characteristics of the perturbed motion obtained by integrating the motion equations of the aircraft numerically. In these equations, nominal aerodynamic parameters at the nodes of the corresponding tables are changed to random values that do not exceed in modulus the given 25 ÷ 50 percents of nominal values at these nodes.
Fundamentally, the MPA algorithm assumes that the vector of unknown parameters η is random on the set of possible flights. We assume that the a priori statistical-generator for computer generated random vectors η and ξ k is given. This generator makes the algorithm estimating components of the vector η (the identification algorithm) Bayesian. Further, for particular calculations, we assume that random components of the mentioned vectors are distributed uniformly and can be called by the standard Random program in Turbo Pascal.
The MPA algorithm provides the approximation method we implement with the multidimensional power series of the vector E(η|Y N ) of the conditional mathematical expectation of the vector η if the vector of measurements Y N is fixed and a priori statistical data on random vectors η and ξ k are given.
The vector E(η|Y N ) is known to be optimal, in root-mean-square sense, estimate of the random vector η.
We describe the steps of operation of the MPA algorithm when it identifies the vector η [6,7].
Step 1. Suppose d is a given positive integer number and the set of integer numbers a 1 , ..., a N consists of all nonnegative solutions of the integer inequality a 1 + ... + a N ≤ d, the number of which we denote by m(d, N). The value m(d, N) is given by the recurrent formula proved by We obtain the vector W N (d) of dimension m(d, N) × 1, the components w 1 , ..., w m (d, N) of which are all possible values y a 1 1 ...y a n N of the form that represent the powers of measurable values.
Then, we construct the base vector Step 2. We use a known statistical generator of random vectors η and ξ k to solve repeatedly the Cauchy problem for Eq.(1) for given initial conditions x(0), a control law u(t) and various realizations of random vectors η and xi k .
We apply the Monte-Carlo method to find the prior first and second statistical moments of the vector V(d, N), i.e., the mathematical expectationV(d, N), and the covariance matrix Implementation of step 2 is a learning process for the algorithm, adjusting it to solve the particular problem described by Eqs. (1) and (2).
Step 3. For given d and N and a fixed vector Y N , we assign the vector η(W N (d)) to be the solution to the estimation problem. This vector gives an approximate estimate of the vector E(η|Y N ) that is optimal in the root-mean-square sense on the set of vector linear combinations of components of the vector W N 1 (d) The vectorV(d, N) and the matrix C V (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 coefficients λ(a 1 , ..., a N ) for (1.1). Moreover, we determine the matrix C(d, N), which is the covariance matrix of the estimation errors for the vector E(η N |Y N ) of conditional mathematical expectation estimated by the vector η(W N (d)).
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 W N (d) is close to linear combination of its previous components. To implement the recursion, we process the components of the vector W N (d) one after another. However, the adjustment of the algorithm performed by applying the Monte-Carlo method to find the vectorV(d, N) and the matrix C V (d, N) takes into account a priori ideas on stochastic structure of components of the whole set of possible vectors W N (d) that can appear in any realizations of the random vectors η and ξ k allowed by the 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 differ 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.
In [6], a multidimensional analogue of the K. Weierstrass theorem (the corollary of the M. Stone theorem [9]) is used to prove that when the integer d is increases then the error estimates of the vector E(η|Y N ) the vector | η(W N (d)) − E(η|Y N )| tend to zero uniformly on some region. Formulas of the recurrent algorithm are given and justified in [6,7] and in the Appendix.
This scheme for the MPA algorithm operation shows that it can be applied to identify parameters of almost any dynamic system provided that the structures of the motion equations and measurements of form (0.1) and (0.2) and prior statistical generators of random unknown parameters and errors of measurements are given. The MPA algorithm is devoid of the above listed limitations and drawbacks, which gives it substantial advantages over NASA identification algorithms. Apart from errors of computations, the algorithm does not add any other errors (such as errors due to linearization of nonlinear functions) into the identified parameters. Therefore, one should expect that the priori spread of identifiable parameters to be always greater than the posterior spread. This is why we can use iterations.
Let us compare the sequential steps of the standard discrete Kalman filter and the MPA algorithm.
(1) The Kalman filter identifies the vector η, which can be represented by part of components of the state vector of the linear dynamic system for the observations that linearly depend on state vectors. The a priori data are the first and second moments of components of random initial state vectors, uncorrelated random vectors of perturbations and observation errors.
We need these data for sequential (recurrent) construction of the estimation vector that is root-mean-square optimal. Usually assigned, a priori data can be also determined by the Monte-Carlo method if the complex mechanism of their appearance is given.
(2) To find an asymptotic solution to the nonlinear identification problem, the MPA algorithm, unlike the Kalman filter, requires a priori statistical data on both the initial and all hypothesized future state vectors of the dynamic system and observations. These a priori data are represented by the first and second statistical moments for the random vector V(d, N): the vector V(d, N)) and the matrix C V (d, N). These moments are calculated using the Monte-Carlo method. However, there are cases when they can be obtained by numerical multidimensional region integration.
(1.1) Once conditions from (1) are met, the Kalman filter constructs the recurrent process, at every step of which the current estimation vector optimal in the root-mean-square sense and the covariance matrix of errors of the estimate are calculated.
(2.1) Based on (2), the MPA algorithm implements the recurrent computational process that do not require matrix inversion. At each step of the process, we construct i. the current estimation vectorη(W N (d)) linear with respect to components of the vector W N (d) and optimal in the root-mean-square sense on the set of linear combinations of components of this vector; moreover, the uniform convergenceη( is attained on some region, ii. the current covariance matrix of estimation errors (we emphasize that known numerical methods of constructing approximations of the vector of nonlinear estimates cannot calculate current covariance matrices of estimation errors).
Implementation of items 2 and 2.1 makes the MPA algorithm more efficient than any known linear identification algorithm since it i. does not involve linearization, ii. does not apply variants of the Newton method to solve systems of nonlinear algebraic equations, iii. forms the estimation vector that tends uniformly to the vector of conditional mathematical expectation for the growing integer d, iiii. obtains the covariance matrix of estimation errors.
It is worth emphasizing that in this work we just develop the fundamental ground of computational technique for solving the complex problem of aircraft parameter identification.

Testing the MPA algorithm: Problem reconstruction (identification of parameters) for the attractor from units of an electrical chain
We consider the boundary inverse problem for the attractor whose equations are presented in [8 ]. The three parameters are the initial conditions: The six parameters η 3+i , i = 1, ..., 6 correspond to combinations of the inductance, the resistances and the two capacitances of a circuit.
The equations of the mathematical model of the circuit take the following form [8]: where X 1 [k] corresponds to a voltage, X 2 [k] to a current and X 3 [k] to another voltage.

Identification of aerodynamic coefficients of the pitching motion for an pseudo f-16 aircraft
We illustrate efficiency of offered MPA algorithm on an example of identification of 48 dimensionless aerodynamic coefficients for the aircraft of near F-16. The aircraft we shall conditionally name " pseudo F-16 ". The term "near" is justified by that, what is the coefficients are taken from SIDPAS [1], but are perturbed by addition of some random numbers.
The tables resulted below, show, that errors of identification are small also modules of their relative values do not surpass several hundredth. The considered problem corresponds to minimization of object function of 48 variables, which it is made of the sum of squares of differences of actual and computational angles of attack, g-load, pitch angles, observable with frequency 10 hertz during 25 sec. flight of the aircraft maneuvering in a vertical plane.

Pitching motion equations
We use the rectangular coordinate system XYZ adopted in NASA. Then for the unperturbed atmosphere and conditions V = const, pitching motion equations have the form [1]: where V=300 ft/sec,H=20000 ft, α is the angle of attack, N Z is the g-load, which is the vector of aerodynamic forces projected onto the axis Z and divided by the weight of the aircraft, M Y is the vector of the moment of aerodynamic forces projected onto the axis Y, ω is the vector of the angular velocity of the aircraft projected onto the axis Y,θ is the angle between the the axis X and the horizontal plane, q is the value of the dynamic pressure, G is the weight, J Y is the moment of inertia with respect to the axis Y, S is the area of the surface generating aerodynamic forces, b is the mean aerodynamic of the wing, C Z (α, δ) and C m (α, δ) are dimensionless coefficients of the aerodynamic force and moment,δ s is the angle of the stabilator devlection measured in degrees.

Parametric model aerodynamic forces and moments
Nominal values of 4 functions of the angle of attack C Z 0 (α), C m 0 (α), C Z q (α), C m q (α) are given with the argument step (55 − 1)/12 degree at 12 nodes (Table 1) To determine values of functions between the nodes, we use linear interpolation. Having analyzed Table 1, we can see that functions nonlinear. Table 2 confirms this visual impression. In it increments are presented 4 functions on each step of Table 1. Apparently, increments noticeably vary.
We study the identification problem for the perturbed analogues of the functions C Z 0 (α), C m 0 (α), C Z q (α), C m q (α). The number of nominal coefficients that determine these functions is 12+12+12+12 = 48. Let us single out the problem which is the most complex for the MPA algorithm, when the actual coefficients differs from the nominal coefficients by the unknown bounded by the prior limits value η i at each point of the table. Then, for accumulated results of measurements of parameters of the perturbed motion, the MPA algorithm is to estimate 48 components of the vector of random estimates, -the vector of differences between actual and nominal coefficients.
Suppose ϑ i and B i are the i-th components of the nominal and actual (perturbed) vectors of aerodynamic coefficients , i = 1, ..., 48, i.e. the number of actual coefficients to be identified is 48 in this case. We assume that the parametric model holds. The vector η serves as the vector of perturbations of nominal data errors of aerodynamic parameters, and identification yields the estimates of its components. We give the structure of these components by the formula The positive number ρ i gives the maximum value that, by identification conditions, can be attained by the ratio of the absolute values of the random value of perturbations η i and nominal coefficients ϑ i .

Transient processes of characteristics of nominal motions
We wish to identify-estimate -during one test flight the 48 unknown aerodynamic coefficients for 12 nodes-12 the set angles of attack α i , i = 1..., 12. For a testing maneuver the characteristics α(t), N Z (t), θ(t) of Transient Processes are carrier of information of the the identified coefficients. Therefore during flight the aircraft should "visit" vicinities of angles of attack −10 • ≤ α ≤ 45 •  Table 3. The characteristics α(t), N Z (t), θ(t) of the nominal motions for the chosen control law δ s (t).

Estimating identification accuracy of 48 errors of aerodynamic parameters of the aircraft
Primary task of MPA algorithm consists in identification -estimation-48 increments of 4 functions. If entry conditions and increments are determined, values of the unknown coefficients follow from obvious recurrent formulas.
To estimate the accuracy, we assume that the current values of α, N Y , θ are measured every 0.1 sec. during 25 seconds .We assume that random errors of measurement represent the discrete white noise bounded by the true measurable value multiplied by the given value ǫ.
An amount of the primary observations equal 3*250=750.  Table 4. The Relative errors of the identifications of C Z 0 (α i ) by ρ = 0.25  Table 5. The Relative errors of the identifications of C m 0 (α i ) by ρ = 0.25 We compress primary observations for a smoothing the high-frequency errors and reduction of a dimension of matrixes covariance . The file of the primary observations is divided into 12 groups and as an input of the algorithm of the identification the vector of the dimension 12 × 1 serves. Components of this vector are the sums of elements of each of 12 groups.
To characterize the accuracy of identification of the random parameter η i the degree of perturbation of the aerodynamic coefficients ϑ , we determine the relative errors of estimation (η i −η i )/η i for every component the identifiable functions . The relative errors designate δ(C Z 0 (α i )), δ(C m 0 (α i )), δ(C Z q (α i )), δ(C m q (α i )), i = 1, ..., 12.
Apparently, relative errors of identification are small and do not surpass several hundredth at ρ = 0.25  Table 6. The Relative errors of the identifications of C Z q (α i ) by ρ = 0.25  Table 7. The Relative errors of the identifications of C m q (α i ) by ρ = 0.25  Table 8. The Relative errors of the identifications of C Z 0 (α i ) by ρ = 0.50  Table 9. The Relative errors of the identifications of C m 0 (α i ) by ρ = 0.50  Table 10. The Relative errors of the identifications of C Z q (α i ) by ρ = 0.25  Table 11. The Relative errors of the identifications of C m q (α i ) by ρ = 0.50

Conclusions
The presented data show that the multipolynomial approximation algorithm can provide a computational ground for developing an efficient parameter identification technique for the nonlinear dynamic system, including identification of aerodynamic parameters of an aircraft. We emphasize that tables characterizing a sufficiently high accuracy of aerodynamic parameter identification are obtained when there are no iterations and d = 1, which corresponds to the case when the estimation vector(ϑ + η)(W N (d)) is represented by the vector linear combination of measured data that is optimal on the family of linear operators over the vector of measurements. This is due to good (in terms of the identification problem) properties of the parametric system of equations of the pitching motion of the "pseudo F-16 " aircraft. It can become much more complicated when it comes to the identification problem of the parametric system of equations of complete (spatial) motion of the aircraft. In such case, we may need to use polynomials of the power d > 1 and increase requirements on the computer performance and RAM. This was the case for identification attempts made for some parameters of F-16 complete motion equations. We emphasize that the inputs of the MPA algorithm we considered were not real (were not the results of operation of real sensors of the aircraft during its test flight); they were determined by mathematical simulation -by means the numerically integrations motion equations for perturbed parameters of aerodynamic forces and moments.

A.1. An algorithm fundamental (AF)
We consider the algorithm fundamental (AF) for solving the problem of finding the estimate of the vector E(η|Y N ) that is optimal in the root-mean-square sense. This vector is known to be the estimate optimal in the root-mean-square sense of the vector η once the vector Y N is fixed. Therefore, it is justified that it is the vector of conditional expectation that AF tends to estimate.
We construct AF that ensures polynomial approximation of the vector E(η|Y N ). To do this, we find the approximate estimate of the vector E(η|Y N ), which is linear with respect to components of the vector W N (d) and optimal in the root-mean-square sense. We denote the vector of this estimate by η(W N (d)) . To obtain the explicit expression for the estimation vector, we calculate elements of the vector V(d, N) and the covariance matrix C V (d, N) that are the first and second (centered) statistical moments for the vector V(d, N). These vector and matrix can be divided into blocks of the following structure The right-hand sides of these blocks are the first and second (centered) statistical moments calculated by the Monte-Carlo method. However, their left-hand sides also serve as the first and second (centered) statistical moments of components of the vector of conditional mathematical expectations. Hence, we can use mathematical models of form (0.1) and (0.2) to find these statistical moments experimentally for vectors of conditional expectations as well. This obvious proposition gives us the basis for practical implementation of the computational procedure of estimating the vector of the conditional expectation.
We introduce We also introduce The lemma follows from the identity Corollary of the lemma. For the vector E(η|Y N , the vector η(W N (d)) is the estimate optimal in the root-mean-square sense among the set of estimates linear with respect to components of the vector W N (d).IfQ N (d) > 0, the estimation vector is unique and The covariance matrix C(d, N) of estimation errors of the vector E(η|Y N ) is given by the formula If Q N (d) ≥ 0, the vectors that provide linear and optimal in the root-mean-square sense estimate are not unique; however, the variances of components of the difference between these vectors are zeros.
Formula (A.1) gives explicit expressions for the vector coefficients of the form λ(a 1 , ..., a N ) in (1.1). To find these relations, we open the explicit expressions for components of the vector W N (d) and the right-hand side of (A.1) and equate them to the right-hand side of formula (1.1).
We consider asymptotic estimation errors when we use (A.1). Suppose the vector Y N is fixed. We assume that the vector E(η|Y N is given by the function of Y N on some a priori region that is a compact; the function is continuous on this region. Then, the following theorem holds. Theorem. Proof. The multidimensional analogue of the K. Weierstrass theorem, which is the corollary of the M. Stone theorem [9], states that for any number ε > 0 there exists a multidimensional polynomial P(W N (d ε )) such that We can rewrite this relation as We assume that C is the covariance matrix of the random vector P(W N (d) − E(θ|Y N ) : It follows from (A.7) that By construction, the vector η(W N (d)) provides the estimate of the vector η that is linear with respect to components W N (d) and optimal in the root-mean-square sense. However, it follows from the lemma that for any other non-optimal linear estimate, including estimates of the form P(W N (d), the relation C ≥ C(d, N) holds. Hence, taking into account (A.8), we obtain Proposition (A.9) is equivalent to (A.6) if we recall that where p(η, Y N ) is the joint probability density of the random vectors η and Y N . The theorem is proved.
Thus, by (A.1), AF determines the vector series that, with the increasing number m(d, N) of its terms, approximates the vector of conditional mathematical expectation of the vector θ of the estimated parameters with an arbitrary uniformly small root-mean-square error.

A.2. Recurrent (Realizable) MPA algorithm
To use formula (A.1), we need to find the matrix inverse to the matrix Q N (d). When the dimension m(d, N) × m(d, N) of the matrix Q N (d) is high and Q N (d) is close to the singular matrix, it is difficult to calculate elements of the inverse matrix. Below, we give the recurrent computational process based on the principle of decomposing observations, described in [6,7]. Above, we specified the vector W N (d) of dimension m(d, N) × 1 with the components w 1 , ..., w m(d,N) . The computational process consists of m(d, N) successive steps. At each step, we use new updated prior data to find the new estimate of the vector θ and perform the prediction, which provides estimates for the rest part of the observation vector. At the same time, we determine the covariance matrix of the estimation errors attained at this step. There is no prediction at the last m(d, N)− th step, and the vector θ is refined for the last time. At step 1, we use the particular case of formulas of form (A.1) and (A.5) to calculate the vector V 1 that estimates the vector V 1 and is optimal in the root-mean-square sense and linear with respect to w 1 , and the covariance matrix of the estimation errors C(V 1 ).
The estimation vector is formed of the estimate of the vector of conditional mathematical expectation E(η|Y N ) and the vector of dimension m(d, N) − 1) × 1. Once we fix the value w 1 , the latter becomes the vector of statistical prediction of the mean values of "future" values w 2 , ..., w m(d,N) . We emphasize that calculations performed at step 1 are based on the preliminary found prior ,V(d, N), C V (d, N).
Suppose steps 1, ..., i of the computational process yielded the vectorV(d, N) and the matrix C(V i ) after the values w 1 , ...w i , wi were fixed. At step i + 1, we have from the particular case of formulas (A.1) and (A.5) the vectorV i+1 that estimates the vector V i+1 and is optimal in the root-mean-square sense and linear with respect to w 1 , ..., w i+1 , and the covariance matrix C(V i+1 ) of estimation errors. The vectorV i+1 is still formed of the estimate of the vector of conditional mathematical expectation E(η|w 1 , ..., w i+1 (first r components of the vectorV i+1 ) and the vector of statistical prediction of mean values of "future" -h values w i+2 , ..., w m(d,N) after w 1 , ..., w i+1 (the rest m(d, N) − (i + 1) components of the vectorV i+1 ) are fixed. We emphasize that calculations at step i + 1 are based on the preliminary foundV i and C(V i ), which can be naturally called the first and second statistical moments for "future" random values w i+1 , ..., w m(d,N) . These vectors and matrices represent a priori data on statistical moments of components of the vector V i+1 before the algorithm receives the value w i+1 at its input.
Recurrent formulas that corresponds exactly to the above given qualitative description of the computational process have the form where the scalar z w i+1 is the (r + 1)-th component of the vectorV i , the scalar is the linear and optimal in the root-mean-square sense estimate of the component after the algorithm has processed the components w 1 , ..., w i ,V −1 i is the vector obtained from the vectorV i by eliminating its component z w i+1 , the scalar q i is the (r + 1)-th diagonal element of the matrix C(V i ), which is the variance of the estimation error of the component w i+1 after components w 1 , ..., w i were processed, C(V i ) 1 is the matrix formed of C(V i ) after the (r + 1)-th row vector and (r + 1)-th column vector were excluded, and b i is the (r + 1)-th column vector of the matrix C(V i ) with its (r + 1)-th component deleted.
If the scalar q i turned out to be close to zero, the component w i+1 corresponds to a linear combination of components w 1 , ..., w i . Then, w i+1 do not give any new information on θ and should be excluded from the computational process. Note that the sequence of random variables like (w i+1 − z w i+1 ) forms an updating sequence. The upper left block of the (r × r)-matrix C(V i ) includes the covariance matrix C(d, i) of estimation errors of the vector E(η|w i , ..., w i ) after the algorithm processed the vector W i (d). To test this MPA algorithm, we solved numerically several problems of estimating the components of the state vector for essentially nonlinear dynamic systems. The estimated components are unknown random constant parameters η 1 , ..., η r of the dynamic system.
As for particular applied problems, we considered smoothing problems and the filtration problem.
In the above examples, we applied the Monte-Carlo method for the number of random realizations lying within 5000 -10000. This number does not affect the estimation errors provided by the MPA algorithm significantly. The estimated random parameters are assumed to be statistically independent and are a priori uniformly distributed. The value of the root-mean-square deviation σ(i, theo) is determined theoretically by calculating variances :the diagonal elements of the covariance matrix C(d, N). The value of the root-mean-square deviation σ(i, exp) is obtained experimentally by applying the Monte-Carlo method for 5000 realizations. Experimental and theoretical root-mean-square deviations almost coincide, which proves that the above given formulas of the MPA algorithm are correct.

Acknowledgments
This study was financially supported by the Russian Foundation for Basic Research (grant no. 10-08-00415a).