CN111193528B - Gaussian filtering method based on non-linear network system under non-ideal condition - Google Patents

Abstract

The invention relates to a Gaussian filtering method of a nonlinear network system based on a Gaussian filtering method of the nonlinear network system under a non-ideal condition. The invention aims to solve the problems that the existing method does not consider the related noise, one-step random delay measurement and data packet loss which may occur in a nonlinear network system, and the problem that the estimation accuracy of a filter is reduced or even diverged due to model-based linear approximation or neglect of delay measurement. The Gaussian filtering method based on the nonlinear network system under the nonideal condition comprises the following steps: firstly, establishing a system model and a sensor measurement model; step two, providing hypothesis and lemma; thirdly, designing a Gaussian filter based on the second step; and step four, based on a third-order sphere diameter volume rule, approximating the Gaussian weighted integral in the step three to obtain a numerical form of the designed filter. The invention can be applied to the technical field of spacecraft and aircraft navigation.

Description

Gaussian filtering method based on non-linear network system under non-ideal condition

Technical Field

The invention relates to a Gaussian filtering method of a nonlinear network system, in particular to a state estimation method for a nonlinear system with correlated noise, one-step random delay measurement and data packet loss.

Background

In recent years, the estimation problem of network systems has attracted extensive attention^[1-3]([1]L.Schenato，“Optimal estimation in networked control systems subject to random delay and packet drop，”IEEE transactions on automatic control，vol.53，no.5，pp.1311，2008.[2]W.A.Zhang，L.Yu，G.Feng，“Optimal linear estimation for networked systems with communication constraints，”Automatica，vol.47，no.9，pp.1992-2000，2011.[3]R.Caballero-

Hermoso-Carazo, J.Linares-P rez, "Optimal state estimation for network Systems with random parameter matrices, corrected non-esses and delayed measurements," International Journal of General Systems, vol.44, No.2, pp.142-154, 2015.). Kalman Filter (KF)^[4](R.E.Kalman，“A new approach to linear fiteration and prediction schemes, "Journal of basic Engineering, vol.82, No.1, pp.35-45, 1960") are effective methods for solving the state estimation problem. However, it is only suitable for ideal linear systems. In fact, non-linearity, correlated noise, random delay and packet loss are ubiquitous. In this context, the estimation problem of a non-linear network system with synchronous correlated noise and one-step random delay measurement and multi-packet loss is considered.

For non-linear systems, there are many estimation methods. For systems with weak non-linearity, Extended KF (EKF) based on first order Taylor series expansion can achieve acceptable accuracy^[5](Y.Bar-Shalom，X.R.Li，T.Kirubarajan，“Estimation with applications to tracking and navigation:theory algorithms and software，”John Wiley&Sons, 2004.). Differential Filters (DDF) unlike EKF algorithms which require computation of the Jacobian matrix^[6](M.

N.k. poulsen, o.ravn, "New definitions in state estimation for nonlinear systems," Automatica, vol.36, No.11, pp.1627-1638, 2000.) the problem of the EKF algorithm possibly falling into local linearization is overcome by approximating the nonlinear function using stirling interpolation. Based on the fact that it is easier to approximate a probability distribution than to approximate any non-linear function or transformation, a series of sigma point filters, such as Particle Filters (PF) have been proposed^[7](A.Doucet, S.Godsill, C.Andrieu, "On sequential Monte Carlo sampling methods for Bayesian filtering," Statistics and computing, vol.10, No.3, pp.197-208, 2000.), Unscented KF (UKF)^[8](S.J.Julie, J.K.Uhlmann, "unknown filtering and nonlinear evaluation," Proceedings of the IEEE, vol.92, No.3, pp.401-422, 2004.) and the volume KF (CKF)^[9](I.Arasaratnam, S.Haykin, "Cubauture kalman filters," IEEE Transactions on automatic control, vol.54, No.6, pp.1254-1269, 2009.). Meanwhile, in order to adapt to different application scenarios, many improved versions of the above algorithm are proposed. E.g. based on second order Taylor expansionKF, unscented PF, DDF with second order stirling interpolation, adaptive UKF, high order CKF, square root CKF, etc. Because the UKF and the CKF have simpler forms and the numerical stability of the CKF algorithm is much higher than that of the UKF, the CKF algorithm is widely applied to the field of practical engineering, such as target tracking and mobile terminal positioning.

For systems with correlated noise, in^[10](X.X.Wang, Y.Liang, Q.Pan, et al, "A Gaussian adaptive iterative filter for nonlinear systems with corrected noise," Automatica, vol.48, No.9, pp.2290-2297, September, 2012.), a new pseudo process equation was reconstructed in which the process noise is uncorrelated with the observation noise. In that^[11](X.X.Wang, Y.Liang, Q.Pan, et al, "Design and implementation of Gaussian filter for nonlinear system with random delayed measurements and corrected reactions," Applied Mathematics and calculation, vol.232, pp.1011-1024, 2014.) A GASF of nonlinear system is proposed based on the projection theorem. Then, [11 ]]The method of (1) for designing a Gaussian filter for a nonlinear system with random delay measurements and associated noise^[12](G.B. Chang, "Comments on" A Gaussian adaptive feedback filter for nonlinear systems with corrected noises, "Automatica, vol.50, No.2, pp.655-656, February, 2014.). In that^[13](G.B.Chang，“Alternative formulation of the Kalman filter for correlated process and observation noise，”IET Science，Measurement&Technology, vol.8, No.5, pp.310-318, September, 2014.) [13]And [14]The two filtering algorithms in (a) have proven to be theoretically equivalent in a linear system. In that^[14](H.Yu，X.J.Zhang，S.Wang，et al.，“Alternative framework of the Gaussian filter for non-linear systems with synchronously correlated noises，”IET Science，Measurement&Technology, vol.10, No.4, pp.306-315, July, 2016.) two alternative frameworks were developed based on state enhancement and conditional Gaussian distributions and demonstrated [14 []And [10-11]The equivalence of the algorithm (c). Then, at [15 ]](K.Zhao，S.M.Song，“Alternative framework of the Gaussian filter for non-linear systems with randomly delayed measurements and correlated noises，”IET Science，Measurement&Technology, vol.12, No.2, pp.161-168, 2017.) and [16](S.L.Sun, L.Xie, W.Xiao, et al, "Optimal linear estimation for systems with multiple packet routes," Automatica, vol.44, No.5, pp.1333-1342, May,2008.) the alternative framework of conditional Gaussian distributions is extended to nonlinear systems with correlated noise and one-step delay measurements.

For systems with one-step random delay measurement, the delay is typically described using a bernoulli-distributed random variable. For systems with multiple packet losses, zero input compensation, hold input compensation, and prediction compensation are separately proposed. For the problem of observed packets sent from the sensor to the filter, one or more packets arriving at the filter within one sampling period are discussed in [17] (S.L. Sun, "Optimal line filters for discrete-time Systems with random delay and lost measures with/without time stages," IEEE Transactions on Automatic Control, vol.58, No.6, pp.1551-1556, June,2013.), and [18] (S.L. Sun, G.H. Wang, "model and estimation for networked Systems with multiple transmission delays and packet locations)," Systems & controls, vol.73, pp.6-16, 20146, 2014. In [19] (S.L.Sun, J.Ma, "Linear estimation for network control systems with random transmission delays and packet routes," Information Sciences, vol.269, pp.349-365, June,2014.), a more general case is considered in which data packets are transmitted from the sensor and the controller to the filter and the actuator, respectively. In the above case, when packet loss occurs, the measurement is invalid. To avoid this, a hold-in compensation mechanism was developed. In [20] (y.liang, t.w.chen, q.pan, "Optimal linear state estimator with multiple packet drop," IEEE Transactions on Automatic Control, vol.55, No.6, pp.1428-1433, June,2010.), in the case where there is a plurality of packets lost, input compensation is kept for the best linear estimation problem. In [21] (j.ma, s.l.sun, "Information Fusion estimators for systems with multiple sensors of differential packet drop rates," Information Fusion, vol.12, No.3, pp.213-222, July,2011.), the same approach is considered in the centralized and distributed Fusion estimation problem in the sense of Linear Minimum Variance (LMV). In [22] (S.L.Sun, L.Xie, W.Xiao, et al, "Optimal linearity estimation for systems with multiple packet routes," Automatica, vol.44, No.5, pp.1333-1342, May,2008.), an Optimal linear estimator was developed based on a novel packet loss model. However, keeping the input compensation does not take into account the latest information of the system, and therefore a new compensation mechanism is proposed in which the latest observed predicted value is used as compensation. In [23] (J.Ding, S.L.Sun, J.Ma, N.Li, "Fusion optimization for Multi-Sensor network Systems with Packet Loss Compensation," Information Fusion, vol.45, pp.138-149, January,2019.), a centralized and distributed Fusion estimator was designed based on a new Compensation mechanism. In [24] (e.i. silver, m.a. solis, "alternative look at the constant-gain Kalman filter for state estimation over channels," IEEE Transactions on Automatic Control, vol.58, No.12, pp.3259-3265, Decumber, 2013.), a type of handover estimator is considered in which missing data is replaced by the best estimate.

It is worth mentioning that in [25-27] ([25] j.ma, s.l.sun, "Linear estimators for networked systems with one-step transmission delay and multiple packet drop based prediction compensation," IET Signal Processing, vol.11, No.2, pp.197-204, April,2017.[26] c.zhu, y.xia, l.xie, et al, "" optical line estimation for systems with transmission delay and packet, "IET Signal Processing, vol.7, No.9, pp.814-823, Decumber, 2013.[27] s.l.m. transmission delay and packet," audio transmission, pp.354, 349, latency loss, random loss, multiple data is taken into account. In [27], possible packet loss, delay measurement and compensation values are described by introducing two uncorrelated bernoulli random variables, which ensure that the system can receive one of the three quantities as a measurement value at any time in order to maintain the accuracy of the filtering algorithm. However, in [27], it does not consider the case where the delay measurement and the real-time measurement arrive in synchronization. In [26], the measurement model is changed because adjacent measurements may arrive at the same time, where the measurement received at the most recent time in the data processing center is used as a compensation value for the current epoch, and the superiority of the algorithm in [26] compared to [27] is shown in [26 ]. Further, in [25], a new measurement model is proposed by changing the compensation mechanism in which the measurements of the current epoch are used as a compensator instead of the one-step prediction of the latest measurements received by the data processing center. It has been demonstrated that the algorithm proposed in [25] has higher estimation accuracy and less computational burden than [26], since the model in [25] always uses the latest measurement information.

Because the existing method does not consider that the data processing center can receive two adjacent measurements simultaneously in the nonlinear system, information loss may occur when the actual system processes such problems, thereby affecting the estimation accuracy of the system, and even possibly causing filter divergence.

Disclosure of Invention

The invention aims to solve the problems that the existing method does not consider the related noise, one-step random delay measurement and data packet loss which may occur in a nonlinear network system, and the problem that the estimation accuracy of a filter is reduced or even diverged due to model linear approximation or neglect of delay measurement, and provides a Gaussian filtering method based on the nonlinear network system under the nonideal condition.

The Gaussian filtering method based on the non-linear network system under the non-ideal condition comprises the following specific processes:

firstly, establishing a system model and a sensor measurement model;

step two, providing hypothesis and lemma;

thirdly, designing a Gaussian filter based on the second step;

and step four, based on a third-order sphere diameter volume rule, approximating the Gaussian weighted integral in the step three to obtain a numerical form of the designed filter.

The invention has the beneficial effects that:

the invention provides a state estimation method of a nonlinear network system with synchronous correlated noise, one-step random delay and a plurality of data packet losses, which considers that the system is a general nonlinear system, designs a Gaussian recursive filtering algorithm aiming at the synchronous correlated noise, one-step random delay measurement and data packet losses which may occur in the system, can ensure that the system obtains a high-precision estimation value, avoids the divergence of the system and ensures the stability of the system.

Drawings

FIG. 1 is a diagram of the estimation of the state in the UNGM model by the algorithm of the present invention and the algorithm of document [23 ];

FIG. 2 is a diagram of the estimated root mean square error of the state in the UNGM model for the algorithm of the present invention and the algorithm in document [23 ];

FIG. 3 is a graph of the estimated error of the state in the UNGM model for the algorithm of the present invention and the algorithm of document [23 ];

FIG. 4 is a diagram of the RMS error estimation of the algorithm of the present invention and the algorithm of document [23] with respect to states in a strongly non-linear model;

FIG. 5 is a diagram of the estimation error of the algorithm of the present invention and the algorithm of document [23] with respect to a state in a strongly non-linear model.

Detailed Description

The first embodiment is as follows: the Gaussian filtering method based on the non-linear network system under the non-ideal condition in the embodiment comprises the following specific processes:

firstly, establishing a system model and a sensor measurement model;

step two, providing hypothesis and lemma;

thirdly, designing a Gaussian filter based on the second step;

and step four, based on a third-order sphere diameter volume rule, approximating the Gaussian weighted integral in the step three to obtain a numerical form of the designed filter.

The second embodiment is as follows: the first embodiment is different from the first embodiment in that a system model and a sensor measurement model are established in the first step; the specific process is as follows:

establishing a nonlinear discrete time system model with correlated noise:

x_k+1＝f(x_k)+ω_k (7)

establishing a general nonlinear measurement model:

z_k＝h(x_k)+υ_k (8)

in the formula, x_k+1System state at time k +1, x_kSystem state at time k, x_k,x_k+1∈Rⁿ，RⁿIs an n-dimensional real number space; z is a radical of_kIs the sensor model at time k, z_k∈R^m，R^mIs m-dimensional real number space; f (-) and h (-) are known non-linear functions; omega_k∈RⁿAnd upsilon_k∈R^mIs correlated zero mean white Gaussian noise and has a covariance of

In the formula, delta_klIs the Kronecker delta function, Q_kAnd R_kRespectively process noise and measurement noise covariance, S_kIs the cross-covariance, l is the time l, ω_l∈RⁿAnd upsilon_l∈R^mIs correlated zero mean gaussian white noise;

in the present invention, the following is considered: one step of random delay and packet loss may occur during data transmission; the data packet is sent only once; the estimator may receive a maximum of two measurement data simultaneously. That is, the estimator may receive zero, one or two packets. In practice, the measurement equation can be described by the following model. Considering communication bandwidth, delay measurement and data packet loss, a general nonlinear measurement model is further established as follows:

in the formula, z_kIs the sensor model at time k; z is a radical of_k-1Is the sensor model at time k; z is a radical of_kk-1Is when z is_kThe compensation amount when the compensation is lost is predicted for the measurement value at the k moment in one step; gamma ray_kAnd η_kIs an uncorrelated Bernoulli distribution variable and satisfies

P is the probability;

is an intermediate variable; y is_kIs a k-time sensor model with measurement skew and packet loss.

Other steps and parameters are the same as those in the first embodiment.

The third concrete implementation mode: the second embodiment is different from the first or second embodiment in that assumptions and reasoning are given in the second step; the specific process is as follows:

corresponding hypothesis and lemma are given, wherein the hypothesis is the premise of filter design, and the lemma is used for facilitating filter derivation;

hypothesis 1. hypothesis ω_k,υ_k,γ_kAnd η_kAnd x₀Is not related, and x₀Satisfy the requirement of

In the formula, x₀Is the initial value of the number of the first,

an estimated value of the initial value, E [ ]]To meet expectations (·)^TIs composed of

T is the transpose of the first image,

the initial value corresponds to the covariance;

theory 1.A ═ a_ij]_n×nIs a real valued matrix, B ═ diag { B }₁,…,b_nC and C ═ diag { C }₁,…,c_nIs a diagonal random matrix, defining

In the formula, a_ijIs the ith row and jth column element of the A matrix, [ a ]_ij]_n×nIs the ith row and the jth column element is a_ijN × n matrix of (d), diag denotes arranging the following elements thereof as a diagonal matrix, b₁Is the 1 st element of the diagonal of the matrix, b_nIs the nth element of the diagonal of the matrix, c₁Is the 1 st element of the diagonal of the matrix, c_nFor the nth element of the diagonal of the matrix, E { BAC^TIs a first pair matrix BAC^TMaking a multiplication operation and then asking for an expectation, an Hadamard product defining [ M ] N for arbitrary matrices M and N of the same dimension]_ij＝M_ij·N_ijIt is clear that in the corresponding equation (12), even though A contains uncertainty, (12) still satisfies that A needs to be E [ A ] as long as A is uncorrelated with B, C]Replacing;

for the convenience of derivation hereinafter, the augmentation system is

Then, there are

In the formula, X_k+1In order to augment the system with a wide range of systems,

phi and

are all intermediate variables;

wherein the content of the first and second substances,

and is

Wherein I is an identity matrix, 0 is a zero matrix,

is the expectation of phi, E phi]In order to make the expectation for phi,

is composed of

In the expectation that the position of the target is not changed,

is a pair of

Calculating expectation;

a ≠ b means that a is not related to b. (a) (.)^TWhere represents the same amount as a. T represents the rank when it is the upper-label. Y is_k＝L{y_k,y_k-1,…,y₁Where L { · } represents a linear space formed by a · leaf. 0 and I represent the zero matrix and identity matrix of the appropriate dimensions. E [. C]Where E represents the mathematical expectation. E [ a | b ]]Representing the condition expectation of a under the condition of b. P represents a covariance matrix. N (-) represents the distribution function. Integral operation is represented by ^ integral.

Other steps and parameters are the same as those in the first or second embodiment.

The fourth concrete implementation mode: the difference between this embodiment and one of the first to third embodiments is that a gaussian filter is designed in the third step; the specific process is as follows:

the one-step prediction mean and covariance matrix are given as follows

X_k+1|k＝X_k+1|k-1+K_kε_k (1)

In the formula, X_k+1|k-1In order to perform the two-step prediction,

predicting the covariance matrix, ε, for two steps_kIn order to be a new message,

is an innovation covariance matrix, T is transposed, K_kFor the gain matrix, the following is defined

In the formula (I), the compound is shown in the specification,

is X under the measurement condition of time k-1_k+1And epsilon_kCross covariance matrix of (a);

the mean and covariance matrices after measurement modification are as follows:

X_k+1|k+1＝X_k+1|k+M_k+1ε_k+1 (4)

in the formula (I), the compound is shown in the specification,

as an innovation covariance matrix, ε_k+1To be new, M_k+1Is a gain matrix;

in the formula (I), the compound is shown in the specification,

is X under the measurement condition of time k_k+1And epsilon_k+1Cross covariance matrix of (a);

A. prediction

Theorem 1 for systems (13) - (14), the one-step prediction mean and covariance are as follows

X_k+1|k＝X_k+1|k-1+K_kε_k (17)

X_k+1|k-1And

the compounds are given in (20) and (24);

new message epsilon_kSum covariance

Has been given at the last moment;

and (3) proving that:

the certification process is divided into two parts:

1) the two-step prediction mean and covariance are calculated.

According to GASF, there are

X_k+1|k＝X_k+1|k-1+K_kε_k (19)

Because of omega_k⊥L{y_k-1,…,y₁Thus, E [ omega ]_k|Y_k-1]Is equal to 0, so

Bringing (20) into (19) to obtain (17);

then calculate

From (19) to

A covariance-based definition sum (21) of

Rewriting X_k+1|k-1And with X_k+1Are subtracted to obtain

According to the definition of covariance, there are

Here, the

Substituting (24) into (22) to obtain (18).

2) A gain matrix is calculated.

The gain matrix is defined as

Here, the first and second liquid crystal display panels are,

then calculate

And

first, calculate

Using ε in (46)_k+1For ε_kIs provided with

Here, the

In (31), N₁As follows

In N₁In, remove

All but a few of the amounts are known,

the calculation is as follows,

here, ω_k-1|k-1And upsilon_k-1|k-1The calculation is as follows.

y_k-1And ω_k-1The joint distribution of (A) can be as follows

According to the conditional Gaussian distribution theorem, obtaining

For upsilon_k-1Available as above

Of note is ε_k-1＝y_k-1-y_k-1|k-2And

has been calculated at time k-1;

then calculate

Here, the first and second liquid crystal display panels are,

bringing (27) and (38) into (26) and (26) into (25), K can be obtained_k；

The certification is over.

B. Correction

Theorem 2 for systems (13) - (14), the state corrections for the mean and covariance of the Gaussian filter are as follows:

X_k+1|k+1＝X_k+1|k+M_k+1ε_k+1 (43)

wherein the content of the first and second substances,

in the formula, y_k+1For the k +1 moment sensor model with measurement lag and packet loss, N (-) is a distribution function, x_k+1|kState x at time k +1_k+1The one-step prediction value of (1),

is a state x_k+1One step of predicting covariance, x_kkIs state x at time k_kIs determined by the estimated value of (c),

is state x at time k_kCovariance matrix of (v)_k|kMeasuring noise upsilon for time k_kAn estimated value of (d);

the certification process is divided into two parts:

1) calculating innovation epsilon_k+1Sum covariance

According to new message epsilon_k+1Is defined as

ε_k+1＝y_k+1-y_k+1|k (47)

Here, the

In (48), x_k+1|k,

x_kkAnd

bring (48) into (47), get (46).

Overwrite y_k+1|kIs provided with

Here upsilon_k|kThe calculation process is as follows (36)_k-1|k-1. Then, the user can use the device to perform the operation,

is obviously provided with

Thus, it is possible to provide

To calculate

According to introduction 1, define

E[φAφ^T]＝Φ⊙E[A] (52)

Order to

Then the

Through a series of algebraic operations, obtain

Hereinafter, the term on the right side is calculated (51).

Item 1.

Here, the first and second liquid crystal display panels are,

Item 2.

here, the first and second liquid crystal display panels are,

and is

Can obtain the product

and

Where N is₂And N₁In the same form, but N₁K in (1) needs to be replaced by k +1,

Item 3.

Item 4.

Item 5.

Item 6.

Item 7.

Item 8.

Item 9.

Item 10.

because of upsilon_k+1And upsilon_kIs not related, available

Item 11.

Item 12.

Item 13.

Bringing (61), (66), (68) - (78) into (51) makes it possible to obtain

2) Computing a gain matrix M_k+1.

From (45), only calculation is needed

Here, the first and second liquid crystal display panels are,

can be similar to

And (6) obtaining.

Here, the first and second liquid crystal display panels are,

bringing (82) - (84) into (81) can obtain

Will be provided with

And

bring in (79), can obtain

The certification is complete.

Other steps and parameters are the same as those in one of the first to third embodiments.

The fifth concrete implementation mode: the difference between this embodiment and the first to the fourth embodiment is that the gaussian weighted integral in step four is approximated based on the third-order sphere diameter volume rule to obtain the numerical form of the designed filter; the specific process is as follows:

based on the sphere diameter volume rule algorithm, the numerical implementation of the algorithm provided by the invention is as follows.

And (3) prediction:

1. decomposition of

In the formula (I), the compound is shown in the specification,

one-step predictive covariance matrix, M, for the state at time k_k|k-1Is the intermediate variable(s) of the variable,

augmenting system for time k

The one-step prediction of the covariance matrix,

is an intermediate variable;

in N₁In the middle, let

In the formula, N₁Is a Gaussian distribution, Γ_k,k-1|k-1Being intermediate variables, is the state x of the construct_kAnd noise v_k-1Estimation based on measurements at time k-1, x_k|k-1One-step prediction of state at time k, v_k-1|k-1For measuring the estimated value of the noise at the time k-1, Π_k,k-1|k-1Is the intermediate variable(s) of the variable,

the covariance matrix is predicted for the state at time k in one step,

is a state x_kAnd noise v_k-1Based on the cross covariance matrix measured at time k-1,

covariance matrix, M, for the measured noise at time k-1_k,k-1|k-1Is an intermediate variable;

2. calculating volume points

x_i,k|k-1＝M_k|k-1ξ_i+x_k|k-1,i＝1,…2n (88)

In the formula, x_i,k|k-1Volumetric point, ξ, for time k with respect to a one-step predictor_i、ζ_i、

Sigma points and dimensions of 2n, 4n and 2(n + m), respectively;

for component x in joint estimation_k|k-1The volume point of (a) is,

for component x in joint estimation_k-1|k-1Volume point of (1), X_i,k|k-1Augmenting system federation states for time k

One step of predicting the volume point, δ_i,k,k-1|k-1To correspond to gamma_k,k-1|k-1In x_k|k-1The ith volume point of the part, gamma_i,k,k-1|k-1To correspond to gamma_k,k-1|k-1Is on v_k-1|k-1Ith volume point of the part, M_k,k-1|k-1Being an intermediate variable, Γ_k,k-1|k-1Is a middleThe variable, n is the system state dimension, and m is the measurement noise dimension;

3. calculating propagated volume points

In the formula (I), the compound is shown in the specification,

is x_i,k|k-1The volume points after propagation through the system model,

is x_i,kk-1The volume point, h (-) after propagation through the metrology model is a known non-linear function,

the volume points, f (-) after propagation through the system model are known non-linear functions,

is composed of

The volume points after propagation through the metrology model,

is delta_i,k,k-1|k-1Volume points after propagation through the system model;

definition of

In the formula I₁、l₂、l₃、l₄、l₅、l₆、l₇、l₈、l₉、l₁₀、l₁₁、l₁₂、l₁₃、l₁₄、l₁₅Is an intermediate variable;

4. calculating the corresponding quantities in equations 1 and 2

Wherein the content of the first and second substances,

where ω is_k-1|k-1And upsilon_k-1|k-1Given in (34) and (36);

bring in (103) - (104)

To obtain

For epsilon obtained at the last moment_kAnd

carry in (1) - (2) to obtain X directly_k+1|kAnd

in the formula, x_k|k-1The mean is predicted for one step of the state at time k,

the covariance is predicted for one step of the state at time k,

is a state x under the measurement condition at the time k-1_k+1The cross-covariance matrix with the innovation at time k,

is a state x under the measurement condition at the time k-1_kCross covariance matrix with innovation at time k, Q_kProcess noise at time k;

and (3) correction:

1. decomposition of

In the formula (I), the compound is shown in the specification,

one-step prediction of covariance matrix, M, for state at time k +1_k+1|kIs the intermediate variable(s) of the variable,

estimating a covariance matrix, M, for a state at time k_k|kIs the intermediate variable(s) of the variable,

augmenting system joint states for k +1 moments

The one-step prediction of the covariance matrix,

is an intermediate variable;

in N₂In the middle, let

In the formula, N₂Is a Gaussian distribution, Γ_k+1,k|kIs an intermediate variable, Π_k+1,k|kIs gamma_k+1,k|kThe corresponding covariance matrix is then used as a basis,

the cross covariance matrix is predicted for one step at time k +1,

is a state x under the measurement condition of time k_k+1And measure noise v_kThe cross-covariance matrix of (a) is,

measuring a posterior covariance matrix of the noise at the time k;

2. calculating volume points

x_i,k+1|k＝M_k+1|kξ_i+x_k+1|k,i＝1,…2n (109)

x_i,k|k＝M_k|kξ_i+x_k|k,i＝1,…2n (110)

In the formula, x_i,k+1|kOne-step prediction of the volume point, x, for the state quantity at the time k +1_k+1|kFor one-step prediction of state quantities at the time k +1, x_i,k|kCorresponding product point, x, for the estimated value of state quantity at time k_k|kIs an estimate of the state quantity at time k, X_i,k+1|kAugmenting system joint state quantities for k +1 moments

The volume point is predicted in one step,

augmenting system joint state quantities for k +1 moments

Corresponding to x_k+1|kThe volume point of the part of the volume,

augmenting system joint state quantities for k +1 moments

Corresponding to x_k|kThe volume point of the part of the volume,

is an intermediate variable, X_k+1|kIs a joint state quantity at the time of k +1

One step prediction, δ_i,k+1,k|kIs a combined quantity

Corresponding to x_k+1|kPartial volume point, gamma_i,k+1,k|kIs a combined quantity

Corresponding upsilon_k|kPartial volume points, M_k+1,k|kBeing an intermediate variable, Γ_k+1,k|kIs an intermediate variable;

3. calculating propagated volume points

In the formula (I), the compound is shown in the specification,

is x_i,k+1|kThe volume points propagated through the system model,

is x_i，k|kThe volume points propagated through the system model,

is x_i,k+1|kThe volume points propagated through the metrology model,

is x_i,k|kThe volume points propagated through the metrology model,

is composed of

The volume points propagated through the metrology model,

is composed of

The volume points propagated through the metrology model,

is delta_i,k+1,k|kVolume points propagated through the metrology model;

definition of

Of formula (II) to'₁、l′₂、l′₃、l′₄、l′₅、l′₆、l′₇、l′₈、l′₉、l′₁₀、l′₁₁、l′₁₂、l′₁₃、l′₁₄、l′₁₅、l′₁₆、l′₁₇Is an intermediate variable;

4. calculating the corresponding quantities in equations 4, 5 and 6

In the formula, y_k+1A k +1 moment sensor model with measurement time lag and data packet loss exists;

order to

Then there is

In the formula (I), the compound is shown in the specification,

an innovation covariance matrix at the time k +1, phi is an intermediate variable, Hadamard product, psi is an intermediate variable, and xi is an intermediate variable,

is a state x measured at time k_kThe cross-covariance matrix with the innovation at time k +1,

a cross-covariance matrix of the state at time k +1 and the innovation at time k +1 measured at time k, S_kFor the cross-correlation noise at time k, omega_k|kMeasuring process noise omega for time k_kAn estimated value of (d);

bringing (127) - (128) into

To obtain

From (125) - (126), epsilon is obtained_k+1And

will epsilon_k+1And

carry in (43) - (44) to obtain X_k+1|k+1And

X_k+1|k+1＝X_k+1|k+M_k+1ε_k+1 (43)

other steps and parameters are the same as in one of the first to fourth embodiments.

The first embodiment is as follows:

numerical simulation analysis

To verify the validity of the control algorithm designed by the present invention, two non-linear models were used to test the validity of the proposed algorithm. Wherein, 'original', 'this paper' and '23' represent reference signals, the estimation result of the algorithm proposed herein and the estimation result of the algorithm in [23] are extended to a nonlinear system based on EKF.

1) Single variable non-static growth model (UNGM)

UNGM is represented as follows:

where ω is_kAnd upsilon_kIs zero mean white Gaussian noise and the covariance is satisfied

Initial value of x₀Filter initial value is x ═ 0.3_0|00 and P_0|0Probability satisfied by random variables of Bernoulli distribution

And

a sub-independent monte carlo simulation is performed. And the Error at time k (Error) and the Root Mean Square Error (RMSE) are defined as follows

Here, the

And

showing the true and estimated values of the nth Monte Carlo method at time k, N being the total number of simulations, FIGS. 1-3 are corresponding simulation results, where FIG. 1 illustrates that the estimation results herein are closer to the true values and have better results than the literature [23]]Fig. 2 and 3 show the algorithm and document [23]]The advantages of the algorithm are further illustrated by the estimation error and the estimation root mean square error of the algorithm, the algorithm fluctuation is small, and the algorithm robustness is good.

2) Strongly non-linear model

The strong nonlinear model is given as follows:

z_k＝cos(x_1,k)+x_2,kx_3,k+υ_k (135)

here, x_kAnd z_kIs the system state and quantity measurement, ω_kAnd upsilon_kIs zero mean white Gaussian noise and the covariance is satisfied

Initial state is x₀＝[-0.7 1 1]^T. The initial value of the filter is x_0|0＝[0 0 0]^TAnd P_0|0＝I₃。

And

as described above.

Fig. 4-5 are corresponding simulation results, where fig. 4 and 5 are respectively an estimation error and an estimated root mean square error of the algorithm in this document and the algorithm in document [23], and the simulation results show that the proposed algorithm still has a strong processing capability even for a strong nonlinear system, while the algorithm in document [23] hardly obtains a satisfactory true value, which illustrates that the algorithm has a stronger capability of processing problems when applied to an actual system.

The present invention is capable of other embodiments and its several details are capable of modifications in various obvious respects, all without departing from the spirit and scope of the present invention.

Claims (1)

1. The Gaussian filtering method based on the nonlinear network system under the nonideal condition is characterized by comprising the following steps: the method comprises the following specific processes:

firstly, establishing a system model and a sensor measurement model;

step two, providing hypothesis and lemma;

thirdly, designing a Gaussian filter based on the second step;

step four, based on a third-order sphere diameter volume rule, approximating the Gaussian weighted integral in the step three to obtain a numerical form of the designed filter;

establishing a system model and a sensor measurement model in the first step; the specific process is as follows:

establishing a nonlinear discrete time system model with correlated noise:

x_k+1＝f(x_k)+ω_k (7)

establishing a nonlinear measurement model:

z_k＝h(x_k)+υ_k (8)

in the formula, x_k+1System state at time k +1, x_kSystem state at time k, x_k,x_k+1∈Rⁿ，RⁿIs an n-dimensional real number space; z is a radical of_kIs the sensor model at time k, z_k∈R^m，R^mIs m-dimensional real number space; f (-) and h (-) are known non-linear functions; omega_k∈RⁿAnd upsilon_k∈R^mIs correlated zero mean white Gaussian noise and has a covariance of

In the formula, delta_klIs the Kronecker delta function, Q_kAnd R_kRespectively process noise and measurement noise covariance, S_kIs the cross-covariance, l is the time l, ω_l∈RⁿAnd upsilon_l∈R^mIs correlated zero mean gaussian white noise;

considering communication bandwidth, delay measurement and data packet loss, the nonlinear measurement model is further established as follows:

in the formula, z_kIs the sensor model at time k; z is a radical of_k-1Is a k-1 moment sensor model; z is a radical of_k|k-1Is when z is_kThe compensation amount when the compensation is lost is predicted for the measurement value at the k moment in one step; gamma ray_kAnd η_kIs an uncorrelated Bernoulli distribution variable and satisfies

P is the probability;

is an intermediate variable; y is_kThe k-time sensor model has measurement time lag and data packet loss;

giving assumptions and lemmas in the second step; the specific process is as follows:

hypothesis 1. hypothesis ω_k,υ_k,γ_kAnd η_kAnd x₀Is not related, and x₀Satisfy the requirement of

In the formula, x₀Is the initial value of the number of the first,

an estimated value of the initial value, E [ ]]To meet expectations (·)^TIs composed of

T is the transpose of the first image,

the initial value corresponds to the covariance;

theory 1.A ═ a_ij]_n×nIs a real valued matrix, B ═ diag { B }₁,…,b_nC and C ═ diag { C }₁,…,c_nIs a diagonal random matrix, defining

In the formula, a_ijIs the ith row and jth column element of the A matrix, [ a ]_ij]_n×nIs the ith row and the jth column element is a_ijN × n matrix of (d), diag denotes arranging the following elements thereof as a diagonal matrix, b₁Is the 1 st element of the diagonal of the matrix, b_nIs the nth element of the diagonal of the matrix, c₁Is the 1 st element of the diagonal of the matrix, c_nFor the nth element of the diagonal of the matrix, E { BAC^TIs a first pair matrix BAC^TA multiplication operation is performed and then an expectation is asserted, an Hadamard product;

the augmentation system is

Then, there are

In the formula, X_k+1In order to augment the system with a wide range of systems,

phi and

are all intermediate variables;

wherein the content of the first and second substances,

and is

Wherein I is an identity matrix, 0 is a zero matrix,

is the expectation of phi, E phi]In order to make the expectation for phi,

is composed of

In the expectation that the position of the target is not changed,

is a pair of

Calculating expectation;

designing a Gaussian filter in the third step; the specific process is as follows:

the one-step prediction mean and covariance matrix are given as follows

X_k+1|k＝X_k+1|k-1+K_kε_k (1)

In the formula, X_k+1|k-1In order to perform the two-step prediction,

predicting the covariance matrix, ε, for two steps_kIn order to be a new message,

is an innovation covariance matrix, T is transposed, K_kFor the gain matrix, the following is defined

In the formula (I), the compound is shown in the specification,

is X under the measurement condition of time k-1_k+1And epsilon_kCross covariance matrix of (a);

the mean and covariance matrices after measurement modification are as follows:

X_k+1|k+1＝X_k+1|k+M_k+1ε_k+1 (4)

in the formula (I), the compound is shown in the specification,

as an innovation covariance matrix, ε_k+1To be new, M_k+1Is a gain matrix;

in the formula (I), the compound is shown in the specification,

is X under the measurement condition of time k_k+1And epsilon_k+1Cross covariance matrix of (a);

based on the third-order sphere diameter volume rule in the fourth step, the Gaussian weighted integral in the third step is approximated to obtain the numerical form of the designed filter; the specific process is as follows:

and (3) prediction:

1. decomposition of

In the formula (I), the compound is shown in the specification,

one-step predictive covariance matrix, M, for the state at time k_k|k-1Is the intermediate variable(s) of the variable,

augmenting system for time k

The one-step prediction of the covariance matrix,

is an intermediate variable;

in N₁In the middle, let

In the formula, N₁Is a Gaussian distribution, Γ_k,k-1|k-1Being intermediate variables, is the state x of the construct_kAnd noise v_k-1Estimation based on measurements at time k-1, x_k|k-1One-step prediction of state at time k, v_k-1|k-1For measuring the estimated value of the noise at the time k-1, Π_k,k-1|k-1Is the intermediate variable(s) of the variable,

the covariance matrix is predicted for the state at time k in one step,

is a state x_kAnd noise v_k-1Based on the cross covariance matrix measured at time k-1,

covariance matrix, M, for the measured noise at time k-1_k,k-1|k-1Is an intermediate variable;

2. calculating volume points

x_i,k|k-1＝M_k|k-1ξ_i+x_k|k-1,i＝1,…2n (88)

In the formula, x_i,k|k-1Volumetric point, ξ, for time k with respect to a one-step predictor_i、ζ_i、

Sigma points and dimensions of 2n, 4n and 2(n + m), respectively;

for component x in joint estimation_k|k-1The volume point of (a) is,

for component x in joint estimation_k-1|k-1Volume point of (1), X_i,k|k-1Augmenting system federation states for time k

One step of predicting the volume point, δ_i,k,k-1|k-1To correspond to gamma_k,k-1|k-1In x_k|k-1The ith volume point of the part, gamma_i,k,k-1|k-1To correspond to gamma_k,k-1|k-1Is on v_k-1|k-1Ith volume point of the part, M_k,k-1|k-1Being an intermediate variable, Γ_k,k-1|k-1Is an intermediate variable, n is a system state dimension, and m is a measurement noise dimension;

3. calculating propagated volume points

In the formula (I), the compound is shown in the specification,

is x_i,k|k-1The volume points after propagation through the system model,

is x_i,k|k-1The volume point, h (-) after propagation through the metrology model is a known non-linear function,

is composed of

The volume points, f (-) after propagation through the system model are known non-linear functions,

is composed of

The volume points after propagation through the metrology model,

is composed of

Volume points after propagation through the system model;

definition of

In the formula I₁、l₂、l₃、l₄、l₅、l₆、l₇、l₈、l₉、l₁₀、l₁₁、l₁₂、l₁₃、l₁₄、l₁₅Is an intermediate variable;

4. calculating the corresponding quantities in equations 1 and 2

Wherein the content of the first and second substances,

bring in (103) - (104)

To obtain

For epsilon obtained at the last moment_kAnd

carry in (1) - (2) to obtain X directly_k+1|kAnd

in the formula, x_k|k-1The mean is predicted for one step of the state at time k,

the covariance is predicted for one step of the state at time k,

is a state x under the measurement condition at the time k-1_k+1The cross-covariance matrix with the innovation at time k,

is a state x under the measurement condition at the time k-1_kCross covariance matrix with innovation at time k, Q_kProcess noise at time k;

and (3) correction:

1. decomposition of

In the formula (I), the compound is shown in the specification,

one-step prediction of covariance matrix, M, for state at time k +1_k+1|kIs the intermediate variable(s) of the variable,

estimating a covariance matrix, M, for a state at time k_k|kIs the intermediate variable(s) of the variable,

augmenting system joint states for k +1 moments

The one-step prediction of the covariance matrix,

is an intermediate variable;

in N₂In the middle, let

In the formula, N₂Is a Gaussian distribution, Γ_k+1,k|kIs an intermediate variable, Π_k+1,k|kIs gamma_k+1,k|kThe corresponding covariance matrix is then used as a basis,

the cross covariance matrix is predicted for one step at time k +1,

is a state x under the measurement condition of time k_k+1And measure noise v_kThe cross-covariance matrix of (a) is,

measuring a posterior covariance matrix of the noise at the time k;

2. calculating volume points

x_i,k+1|k＝M_k+1|kξ_i+x_k+1|k,i＝1,…2n (109)

x_i,k|k＝M_k|kξ_i+x_k|k,i＝1,…2n (110)

In the formula, x_i,k+1|kIs k +1One-step prediction of volume point, x, of state quantity_k+1|kFor one-step prediction of state quantities at the time k +1, x_i,k|kCorresponding product point, x, for the estimated value of state quantity at time k_k|kIs an estimate of the state quantity at time k, X_i,k+1|kAugmenting system joint state quantities for k +1 moments

The volume point is predicted in one step,

augmenting system joint state quantities for k +1 moments

Corresponding to x_k+1|kThe volume point of the part of the volume,

augmenting system joint state quantities for k +1 moments

Corresponding to x_k|kThe volume point of the part of the volume,

is an intermediate variable, X_k+1|kIs a joint state quantity at the time of k +1

One step prediction, δ_i,k+1,k|kIs a combined quantity

Corresponding to x_k+1|kPartial volume point, gamma_i,k+1,k|kIs a combined quantity

Corresponding upsilon_k|kPartial volume points, M_k+1,k|kBeing an intermediate variable, Γ_k+1,k|kIs an intermediate variable;

3. calculating propagated volume points

In the formula (I), the compound is shown in the specification,

is x_i,k+1|kThe volume points propagated through the system model,

is x_i,k|kThe volume points propagated through the system model,

is x_i,k+1|kThe volume points propagated through the metrology model,

is x_i,k|kThe volume points propagated through the metrology model,

is composed of

The volume points propagated through the metrology model,

is composed of

The volume points propagated through the metrology model,

is delta_i,k+1,k|kVolume points propagated through the metrology model;

definition of

Of formula (II) to'₁、l′₂、l′₃、l′₄、l′₅、l′₆、l′₇、l′₈、l′₉、l′₁₀、l′₁₁、l′₁₂、l′₁₃、l′₁₄、l′₁₅、l′₁₆、l′₁₇Is an intermediate variable;

4. calculating the corresponding quantities in equations 4, 5 and 6

In the formula, y_k+1A k +1 moment sensor model with measurement time lag and data packet loss exists;

in the formula (I), the compound is shown in the specification,

an innovation covariance matrix at the time k +1, phi is an intermediate variable, Hadamard product, psi is an intermediate variable, and xi is an intermediate variable,

is a state x measured at time k_kThe cross-covariance matrix with the innovation at time k +1,

a cross-covariance matrix of the state at time k +1 and the innovation at time k +1 measured at time k, S_kFor the cross-correlation noise at time k, omega_k|kMeasuring process noise omega for time k_kAn estimated value of (d);

bringing (127) - (128) into

To obtain

From (125) - (126), epsilon is obtained_k+1And

will epsilon_k+1And

carry in (43) - (44) to obtain X_k+1|k+1And

X_k+1|k+1＝X_k+1|k+M_k+1ε_k+1 (43)

Images