CC BY
16
MSC 93E13
DOI: 10.14529/ mmp200210
PARAMETRIC IDENTIFICATION BASED ON THE ADAPTIVE UNSCENTED KALMAN FILTER
V.M. Chubich1, O.S. Chernikova1
Novosibirsk State Technical University, Novosibirsk, Russian Federation E-mails: chubich@ami.nstu.ru, chernikova@corp.nstu.ru
The detailed adaptive unscented Kalman filter algorithm is provided. Step-by-step schemes of filtering algorithms used for the software development are given. Nonlinear filtering algorithm efficiency is investigated with considering an example of a nonlinear continuous-discrete model. The statistic estimator based on the continuous-discrete adaptive unscented Kalman filter with noise is proposed for the nonlinear system parameters estimation. The solution to the problem of solar radiation parameters estimation based on the maximum likelihood method and the adaptive unscented Kalman filter is shown. The obtained results lead to significant improvement of satellite trajectory prediction quality.
Keywords: nonlinear stochastic continuous-discrete system; adaptive unscented Kalman filter; parametric identification; ML method; spacecraft motion model; solar radiation model.
Introduction
In practice, it is often necessary to work with different classes of nonlinear dynamical systems. The description of such system in terms of nonlinear mathematical models provides researchers with additional opportunities in qualitative analysis of objects and allows taking into account the accompanying factors caused by the nonlinear essence of the nature laws.
Consider the following model of a stochastic nonlinear continuous-discrete system in a state space:
where x(t) G Rn is the state process; Rn is the n-dimensional Euclidean space; u(t) G Rr is deterministic control (input) vector; w(t) G Rp is the process noise vector; y(tk+1) G Rm is the measurement (output) vector; v(tk+1) G Rm is the measurement error vector; f (•), h(-) are nonlinear functions. Suppose that
• the random vectors w(t) and v(tk+1) form white Gaussian noises with unknown covariance matrices of system and measurements noises
t G [t0,tN],
y(tfc+i) = h{x{tk+i),e) + v(tfc+1), k = 0,1,..., N - 1,
(1) (2)
E[w(t)] = 0, E[w(t)wT(t)] = Q(t)S(t - t), E [v (tk+i )] = 0, E [v (tk+i )vT (ti+i)] = R(tk+i)Ski, E [v (tk+i )wT (t )]=0, k, i = 0,1,..., N - 1, t G [t0)tN ];
• the initial state x(t0) is normally distributed with parameters
E[x(t0)] = x(t0), E{[x(t0) - x(t0)][x(t0) - x(t0)]T} = P(t0)
and has no correlation with w(t), v(tk+1) for all values of k; • 6 = (91,92, ■■■,0s) G H© is the vector of unknown parameters.
In our work we consider the problem of parametric identification of model (1), (2) taking into account the a priori assumptions.
The same problem for the case when the functions f (•) and h(-) are linear, was solved by R. Mehra, N. Gupta [1] and K. Astrom [2].
Despite that a lot of solutions to this problem are provided in different books and papers, the problem of nonlinear models identification is much more complicated and important in a practical way. One class of the methods to solve the mentioned problem is basically the continuous-discrete extended Kalman filter (CD-EKF) applied to a linearized system [3]. Although the CD-EKF is widely used, this filter has some drawbacks. The filter applies the standard linear Kalman filter technique to linearize a nonlinear model. Hence, it requires the sufficient differentiability of the dynamic state and the susceptibility to biasing and to divergence of the state estimates. This approach is sub-optimal and can easily lead to the divergence. The CD-EKF achieves only the first-order accuracy and produces a good result only if the initial estimation error and disturbing noises are small.
These difficulties can be successfully overcome with such nonlinear filters as the continuous-discrete cubature Kalman filter [4, 5] and the continuous-discrete unscented Kalman filter (CD-UKF) that is used in this research.
S.J. Julier et al. [6] proposed the unscented Kalman filter (UKF) as a derivativefree alternative to the extended Kalman filter in the framework of state estimation. The UKF has been developed for the case of highly nonlinear state estimation problems. The UKF performs a Gaussian approximation with a limited number of points (sigma points), using the unscented transform. This technique is used to linearize a nonlinear function of a random variable via the linear regression based on the points drawn from the prior distribution of the random variable. The UKF has the same computational complexity as the EKF has. The UKF does not require the Jacobians computing and can achieve the second-order accuracy of the Taylor expansion. Modification of the UKF for continuous-discrete models is given in [7].
When solving practical problems, statistical parameters of noise are set inaccurately or they are completely unknown. The presence of outliers in the measurement data makes the further determination of such characteristics complicated. When using the incorrect a priori information about the noise properties of the system and/or the measurements, the obtained estimates may be biased. Usually the covariance matrices of the system and the measurements noises are selected accordingly to the results of some empirical data analysis or the various situations modeling. Often the correct specification of the statistical noise parameters determines the accuracy of the state vector estimation.
One of the possible solutions to this problem is using adaptive methods for the measurement data processing, which, along with the state vector estimation, can restore the statistical characteristics of noises. In this research, the sub-optimal Sage-Husa estimator [8] is combined with the CD-UKF algorithm in order to estimate and improve the statistical properties of the process noise. Such improvement reduces the model error, suppresses the filtering divergence and improves the filtering accuracy.
This paper is organized as follows. In Section 2 the statistic estimator based on the continuous-discrete adaptive unscented Kalman filter (CD-AUKF) with noise is proposed for the nonlinear systems parameters estimation. The application of the solar-radiation model parametric identification algorithm is presented in Section 3. The conclusion is provided in Section 4.
1. Parameter Identification
In this work we investigate the problem of estimating unknown parameters of the model structure that has been obtained from the physical modeling of the process. To obtain a model with good predictive properties we need informative measurement data and a suitable model structure able to describe the dynamics of the process.
Estimation of unknown parameters of mathematical model is carried out according to observation data S and identification criterion x. The a priori assumptions allow using the maximum likelihood (ML) method for the parameters estimation. Under mild conditions, ML estimates have such practically important properties, as asymptotic unbiasedness, consistency, asymptotic efficiency and asymptotic normality.
Let YT = yT(t1),yT(t2), ■ ■ ■ ,yT(tN) be the output signal corresponding to the input
signal U = u(t),t £ [t0,tN] . As the result of an identification experiment a set S =
{U, Y} is generated. According to ML method it is necessary to find such values of <9 parameters, for which
0 = arg min [x(6; S)] = arg min [- In L(0; S)], (3)
where
Nm 1 N-i 1 N-i
X(6; S) = In 2tt + - J] eT(tk+l)Py \tk+l)e(tk+l) + - J] In det Py{tk+l), (4)
k=0 k=0
s(tfc+1) and PY(tk+1) are defined based on the corresponding equations of the CD-AUKF.
Incorrect mathematical models and inaccurate noise statistic properties often lead to the filter divergence. In order to solve this problem, the adaptive filtering technology has been studied extensively [9-12]. The estimation and correction of the unknown time-varying noise statistics are carried out with using the noise statistics estimator. Such estimator can cover the drawbacks of the traditional UKF for the noise statistics of unknown time-varying filter. The Sage-Husa algorithm is a noise estimator, which is easy to understand and simple in terms of computing. This algorithm can estimate the first and the second moments of noise at the same time. Following [10,11], we combine the CD-UKF with the Sage-Husa noise statistics estimator.
Algorithm CD-AUKF
Initialization:
• Set the values [13]
f = 0, 001, n = 2, p = k = 0, b = 0, 998.
• Define the initial values
x{t0 | t0) = x{t0),P{t0 | t0) = P(to), Q(to), R{h).
• Calculate l
I = £2(n + f) ~ ao
l + n
1
(n + 1) + (1 - f2 + n)
, a
a = [ao ,a1, ■ ■ ■, a2n]
2(n +1)
T
= A,i = 1,
2n,
A = (i - [a | ... | a] ) • diag(/30, A, • • • > fan) • (j ~ [a \ ■ ■ ■ \ a] )
T
2n+1
2n+1
For k = 0,N - 1 Prediction:
• Define X(tk+1 | tk) and P(tk+1 | tk) as the result of differential equations (5),(6) integration
d
—x(t | tk) = Xf(t | tk)a,tk < t < tk+1,
d
-Pit | tfc) = Xs{t | ifc)AX; (t | tk) + | tfc) {t | ifc) +
+G(t)QQ(t)GT(t),tfc < t < tfc+1, where the transformed set of vectors is identified as
Xf (t | tfc )= fx0(t | tfc ),u(t)) | f(x$(t | tk ),u(t)) | ■■■
(5)
(6)
I ^x2n(t | tk),u(t)
nx (2n+1)
sigma points xf(t \ tk), i = 1 ,n are computed in accordance with the following formula xf (t | tk) =
x(t | tk), i = 0,
x{t | tk) + + lDi(t \ tk), i = 1, n,
x(
(i | tk) — yjn + lDi_n(t \ tk), i = n + l,2n
(7)
Xf(t I tk) = xf (t I tk) I xf(t I tk) I ■ ■ ■ I xfn(t I tk)]nx(2n+1), Dj is the i-th row of the lower triangular matrix obtained by the Cholesky decomposition
P(t I tk).
Updating:
• Find the set Xf (tk+1 I tk) using (7) with the substitution t = tk+1.
• Calculate
Yh(tk+1 I tk) = [h(xf (tk+1 I tk)) I h(xf (tk+1 I tk)) I ■ ■ ■
■■■ I Mxfn(tk+1 I tk)
mx(2n+1)
e(tk+1) = y(tk+1) - Yh(tk+1 I tk)a,
Tk
1-6 1 - 6fc+!
^(tk+1) = (1 - Tk)-R(tk) + Tk e(tk+1)eT(tk+1)-
1
2n t
J] h(xf (ifc+i | tk)) - Yh(tfc+i | tk)a) (h(xf (tfc+i | tk)) - Yh(tk+i | ifc)a) ],
i=0
Py(tk+i) = Yh(tk+i | tk)AYhT(tk+i | tk) + R(tk+i),
PXY (tk+i) = Xf (tk+i 1 tk )AYhT (tk+i 1 tk)-
• Given these predicted values, the state x(tk+i | tk+i) and covariance estimates P(tk+i | tk+i) are computed according to the equations
K (tk+i) = PXY (tk+i) PY i (tk+i) j
x(tk+i | tk+i) = x(tk+i | tk) + K(tk+i)e(tk+i), P(tk+i | tk+i) = P(tk+i | tk) - K(tk+i)PY(tk+i)KT(tk+i),
r(tk+i) = {GT (tk+i)G(tk+i^ GT (tk+i)j
Q(tk+i) = (i - Tk)QQ(tk)+
+r(tk+i^Tk [K(tk+i)e(tk+i)eT(tk+i)KT(tk+i) + P(tk+i | tk+i)-
2n
f(x?(tk+i | tk),u(tk+i^ - x(tk+i | tk)) f(xS(tk+i | tk),u(tk+i^-
i=0
T T
-X(ifc+1 | tk)) }rT(tk+i).
The cost function (3) is known to have many local optima. There are many algorithms available for this kind of problems, for instance, Newton's method and various quasiNewton methods, which are the local ones. In case of using a gradient based local optimization method there is the large risk that the minimum found is not the global one, unless the initial values are chosen close enough to the global minimum. In general, if the obtained parameters estimates give a bad fit, there is no way to understand if the reason is either the convergence to a local minimum or the insufficient model structure. One approach to solving these problems with many local minima is to use the global optimization methods. To find the optima of the problem (3), in this work the global optimization approach based on the SQP method is used.
2. Simulation Results
Consider the following model of a stochastic nonlinear continuous-discrete dynamic system:
Jt = -jrjT^^t) + 9i(r(t)) + 92(r(t)) +
+0a(r(t),r(t),0)+ w(t), t G [to, tN], (8)
s(tk+i) = r(tk+i) + v(tk+i), k = 0,1,..., N - 1.
Here r(t) = (x(t),y(t),z(t))T is the coordinate vector of the navigation satellite in an inertial coordinate system; r(t) = (VX(t),Vy(t),VZ(t))T is the velocity vector of the navigation satellite in an inertial coordinate system; ^ is the gravitational constant; ME is
the mass of the Earth; ||r(t)|| = \Jx2(t) + y2(t) + z2(t) is the radius of the orbit; Qiijit)) is the perturbations, caused by the non-sphericity of the Earth's geopotential (see, for example, O. Montenbruck [14]); g2(r(t)) is the perturbations, caused by the gravitational influence of the Moon, the Sun and/or the other planets (also see O. Montenbruck [14]); g3(r(t), r(t), 0) is perturbations from the solar radiation (SR); s(tk+1) is the measurement vector.
To compute g3(r(t), r(t), 0) in an inertial coordinate system, the following SR model in the object-centered coordinate system has been used (this model has been prepared in the processing centers of the International GNSS Service):
g(xoc(T),y0C(T),zoc(t), 0) = c • A(xoc(t),y0C(T),zoc(t)) • p-2(xoc(t),yoc(t),zoc(t))• •[xoc(t) • (01 + 02 cos a(xoc(T),y0C(T),zoc(t)) + 03 sina(xoc(t),yoc(t),zoc(t))) + +yRP(t) • (04 + 05 cos a(xoc(T),y0C(T),zoc(t)) + 06 sina(xoc(t),y0C(T),zoc(t))) +
+ZRP(t) • (07 + 08 cos a(xoc(T),y0C(T),zoc(t)) + 09 sin a(xoc(t),y0C(T),zoc(t)))]■ (9)
Here xoc(t),yoc(T),zoc(t) are the coordinates of the satellite in the object-centered coordinate system; c is the factor depending on the form of the satellite, its mass, reflectivity and absorption of the materials of its surface; A(roc(t)) is the eclipse factor; p(roc(t)) the distance between the satellite and the Sun; a(roc(T)) is the argument of the latitude for the navigation satellite.
Usually, the following parameter values are considered (see, for example, J. Kouba [15])
0 1 = (01,01,01,04,01,01,01, <M) = (1, 0, 0, 0, 0, 0, 0, 0, 0)-
As the measurement data we have taken the rapid ephemeris of the GPS from July 14, 2016, obtained by the international GNSS service. In this case, the satellite makes more than one revolution around the Earth (passes through the various light zones). At the initial time, we compute the velocity of the satellite on the basis of rapid ephemeris using Everett interpolation. Estimation of the SR parameters of the model (9) can be carried out using the maximum likelihood method according to the trajectory observations in areas of total illumination and penumbra zones. As the result, we have obtained the following:
02 = (02,022,022,022,02,022,02,022,022) = (1,06906362,0,06858372,0,04729392,
0,11131100, 0,06708731, 0,09272575, 0,09483054, 0,13523427, 0,10924206)^
To compare the predicted trajectories of the PG01 satellite orbital movement with the final ephemeris from July 15, 2016 we have computed
\/j2i=o \\s{tk+i) - S^tk+i)
iN-1
„ v —....................• -, ...
51 = --=-, i= 1,2,
y/N
where ||-1| is the Euclidean vector norm; |s(ik+1), k = 0,1,..., N —1} is the final ephemeris; (s 1(tk+1),k = 0,1,...,N — 1} is the predicted trajectory for the filter equation at 9 1; {s 2(tk+1), k = 0,1,..., N — 1} is the predicted trajectory for the filter equation at 92.
Finally, we have obtained ^ = 8, 6228e—06 km., SS, = 3, 9537e—08 km. Thus, the result of the sunlight SR parameters specification is that it is possible to significantly increase (by two orders of magnitude) the accuracy of the satellite trajectory forecasting. Figure represents the obtained results. The i-th curve in the figure corresponds to ||s(tk+1) — ^+1)11.
Conclusion
In this research, the statistic estimator based on the CD-AUKF with noise is used for the nonlinear continuous-discrete systems parameters estimation. This approach improves the robustness of the conventional UKF with the respect to the variable noise distribution. The results show that in case of uncertain or time-varying noise statistic the adaptive UKF is more efficient than the conventional UKF in terms of the fast convergence and the state estimation accuracy. Such efficiency is achieved due to applying the noise statistic estimator for the noise statistic calibration.
Using the developed parametric identification procedures allows finding the estimates of the solar radiation parameters of the spacecraft motion model in the inertial coordinate system. The accuracy of the navigation satellite motion prediction has been significantly improved.
Acknowledgements. The work was supported by the Ministry of Education and Science of the Russian Federation (project No 2.7996.2017/8.9).
References
1. Gupta N.K., Mehra R.K. Computational Aspects of Maximum Likelihood Estimation and Reduction in Sensitivity Function Calculations. IEEE Transactions on Automatic Control, 1974, vol. 19, no. 6, pp. 774-783. DOI: 10.1109/tac.1974.1100714
2. Astrom K.J. Maximum Likelihood and Prediction Errors Methods. Automatica, 1980, vol. 16, no. 5, pp. 551-574.
3. Jazwinsky A. Stochastic Processes and Filtering Theory. Academic Press, N.Y., 1970.
4. Sarkka S., Solin A. On Continuous-Discrete Cubature Kalman Filtering. IFAC Proceedings Volumes, 2012, vol. 45, no. 16, pp. 1221-1226. DOI: 10.3182/20120711-3-be-2027.00188
5. Arasaratnam I., Haykin S., Hurd T.R. Cubature Kalman Filtering for Continuous-Discrete Systems: Theory and Simulation. IEEE Transactions on Signal Processing, 2010, vol. 58, no. 10, pp. 4977-4993.
6. Julier S.J., Uhlmann J.K., Durrant-Whyte H.F. A New Approach for the Nonlinear Systems. American Control Conference, Seattle, 1995, pp. 1628-1632.
7. Sarkka S. On Unscented Kalman Filtering for State Estimation of Continuous-Time Nonlinear Systems. IEEE Transactions on Automatic Control, 2007, vol. 52, no. 9, pp. 1631-1641. DOI: 10.1109/TAC.2007.904453
8. Julier S.J., Uhlmann J.K. A New Extension of the Kalman Filter to Nonlinear Systems. The 11th International Symposium on Aerospace/Defence, Sensing, Simulation and Controls, 1997, pp. 12. DOI: 10.1117/12.280797.
9. Mohamed A.H., Schwarz K.P. Adaptive Kalman Filtering for INS/GPS. Journal of Geodesy,
1999, vol. 73, pp. 193-203.
10. Qijun Xia, Ming Rao, Yiqun Ying, Xuemin Shen. Adaptive Fading Kalman Filter with an Application. Automatica, 1994, vol. 30, no. 8, pp. 1333-1338. DOI: 10.1016/0005-1098(94)90112-0
11. Sarkka S., Nummenmaa A. Recursive Noise Adaptive Kalman Filtering by Variational Bayesian Approximations. IEEE Transactions on Automatic Control, 2009, vol. 54, no. 3, pp. 596-600. DOI: 10.1109/tac.2008.2008348
12. Izanloo R., Fakoorian S.A., Yazdi H.S., Simon D. Kalman Filtering Based on the Maximum Correntropy Criterion in the Presence of Non-Gaussian Noise. Annual Conference on Information Science and Systems (CISS), Princeton, 2016, pp. 500-505. DOI: 10.1109/ciss.2016.7460553
13. Wei Gao, Jingchun Li, Guangtao Zhou, Qian Li. Adaptive Kalman Filtering with Recursive Noise Estimator for Integrated Sins/Dvl Systems. The Journal of Navigation, 2015, vol. 68, no. 1, pp. 142-161. DOI: 10.1017/s0373463314000484
14. Montenbruck O., Gill E. Satellite Orbits: Models, Methods and Applications. Berlin, Springer,
2000.
15. Kouba J. A Guide to Using International Gnss Service (IGS) Products Geodetic Survey Division Natural Resources, Ottawa, Natural Resources Canada, 2009.
Received October 8, 2019
УДК 51-74 БЭТ: 10.14529/mmp200210
ПАРАМЕТРИЧЕСКАЯ ИДЕНТИФИКАЦИЯ НА ОСНОВЕ
АДАПТИВНОГО СИГМА-ТОЧЕЧНОГО ФИЛЬТРА КАЛМАНА
В.М. Чубин1, О.С. Черникова1
1 Новосибирский государственный технический университет, г. Новосибирск,
Российская Федерация
Представлен подробный алгоритм адаптивного сигма-точечного фильтра Калма-на. Приведена пошаговая схема алгоритма фильтрации, используемая при решении задачи параметрической идентификации стохастических непрерывно-дискретных систем. На примере математической модели движения навигационного спутника показана эффективность процедуры параметрической идентификации с использованием адаптивного сигма-точечного фильтра Калмана. Полученные результаты позволяют значительно улучшить качество прогнозирования траектории движения спутника.
Ключевые слова: нелинейная стохастическая непрерывно-дискретная система; адаптивный сигма-точечный фильтр Калмана; параметрическая идентификация; метод ИЬ; модель движения космического аппарата; модель солнечного излучения.
Литература
1. Gupta, N.K. Computational Aspects of Maximum Likelihood Estimation and Reduction in Sensitivity Function Calculations / N.K. Gupta, R.K. Mehra // IEEE Transactions on Automatic Control. - 1974. - V. 19, № 6. - P. 774-783.
2. Astrom, K.J. Maximum Likelihood and Prediction Errors Methods / K.J. Astrom // Automatica. - 1980. - V. 16, № 5. - P. 551-574.
3. Jazwinsky, A. Stochastic Processes and Filtering Theory / A. Jazwinsky. - New York: Academic Press, 1970.
4. Sarkka, S. On Continuous-Discrete Cubature Kalman Filtering / S. Sarkka, A. Solin // IFAC Proceedings Volumes. - 2012. - V. 45, № 16. - P. 1221-1226.
5. Arasaratnam, I. Cubature Kalman Filtering for Continuous-Discrete Systems: Theory and Simulation / I. Arasaratnam, S. Haykin, T.R. Hurd // IEEE Transactions on Signal Processing. - 2010. - V. 58, № 10. - P. 4977-4993.
6. Julier, S.J. A New Approach for the Nonlinear Systems / S.J. Julier, J.K. Uhlmann, H.F. Durrant-Whyte // American Control Conference. - Seattle, 1995. - P. 1628-1632.
7. Sarkka, S. On Unscented Kalman Filtering for State Estimation of Continuous-Time Nonlinear Systems / S. Sarkka // IEEE Transactions on Automatic Control. - 2007. - V. 52, № 9. - P. 1631-1641.
8. Julier, S.J. A New Extension of the Kalman Filter to Nonlinear Systems / S.J. Julier, J.K. Uhlmann // The 11th International Symposium on Aerospace/Defence, Sensing, Simulation and Controls. - 1997. - P. 12.
9. Mohamed, A.H. Adaptive Kalman Filtering for INS/GPS / A.H. Mohamed, K.P. Schwarz // Journal of Geodesy. - 1999. - V. 73. - P. 193-203.
10. Qijun Xia. Adaptive Fading Kalman Filter with an Application / Qijun Xia, Ming Rao, Yiqun Ying, Xuemin Shen // Automatica. - 1994. - V. 30, № 8. - P. 1333-1338.
11. Sarkka, S. Recursive Noise Adaptive Kalman Filtering by Variational Bayesian Approximations / S. Sarkka, A. Nummenmaa // IEEE Transactions on Automatic Control. -2009. - V. 54, № 3. - P. 596-600.
12. Izanloo, R. Kalman Filtering Based on the Maximum Correntropy Criterion in the Presence of Non-Gaussian Noise / R. Izanloo, S.A. Fakoorian, H.S. Yazdi, D. Simon // Annual Conference on Information Science and Systems. - Princeton, 2016. - P. 500-505.
13. Gao, Wei. Adaptive Kalman Filtering with Recursive Noise Estimator for Integrated Sins/Dvl Systems / Wei Gao, Jingchun Li, Guangtao Zhou, Qian Li // The Journal of Navigation. -2015. - V. 68, № 1. - P. 142-161.
14. Montenbruck, O. Satellite Orbits: Models, Methods and Applications / O. Montenbruck, E. Gill. - Berlin: Springer, 2000.
15. Kouba, J. A Guide to Using International GNSS Service (IGS) Products Geodetic Survey Division Natural Resources / J. Kouba. - Ottawa: Natural Resources Canada, 2009.
Владимир Михайлович Чубич, доктор технических наук, профессор, кафедра «Теоретическая и прикладная информатика:», Новосибирский государственный технический университет (г. Новосибирск, Российская Федерация), chubich@ami.nstu.ru.
Оксана Сергеевна Черникова, кандидат технических наук, доцент, кафедра «Теоретическая и прикладная информатика: , Новосибирский государственный технический университет (г. Новосибирск, Российская Федерация), chernikova@corp.nstu.ru.
Поступила в редакцию 8 октября 2019 г.
