CN105975747A

CN105975747A - CSTR (Continuous Stirred Tank Reactor) model parameter identification method based on unscented Kalman filtering algorithm

Info

Publication number: CN105975747A
Application number: CN201610272832.8A
Authority: CN
Inventors: 周琪; 于占东; 王焕清; 王巍
Original assignee: Bohai University
Current assignee: Bohai University
Priority date: 2016-04-27
Filing date: 2016-04-27
Publication date: 2016-09-28

Abstract

The invention discloses a CSTR model parameter identification method based on an unscented Kalman filter algorithm. This method first obtains the state space expression of the state component including the parameters to be identified based on the CSTR continuous system model; then, the obtained nonlinear continuous state space expression is discretized with the help of the Euler algorithm, and the corresponding discrete Iterate the model. Finally, the unscented Kalman filter algorithm is used for multiple iterations of identification, and accurate identification results are obtained. The algorithm has good convergence and is easy to combine with existing software, so it has a good prospect of engineering application.

Description

A CSTR Model Parameter Identification Method Based on Unscented Kalman Filter Algorithm

技术领域technical field

本发明涉及一种基于无迹卡尔曼滤波算法的CSTR模型参数辨识方法，属于系统建模与参数辨识技术领域。The invention relates to a CSTR model parameter identification method based on an unscented Kalman filter algorithm, and belongs to the technical field of system modeling and parameter identification.

背景技术Background technique

连续搅拌釜式反应器(CSTR)是化工生产过程中典型的、高度非线性的化学反应系统。在化工生产的核心设备中占有相当重要的地位，在染料、医药试剂、食品及合成材料工业中，CSTR系统得到了广泛的应用。The continuous stirred tank reactor (CSTR) is a typical and highly nonlinear chemical reaction system in the chemical production process. It occupies a very important position in the core equipment of chemical production. In the dyestuff, pharmaceutical reagent, food and synthetic material industries, CSTR system has been widely used.

由于CSTR在化工生产中的重要作用，所以有必要对该过程进行详细的研究。在对CSTR建立模型进行分析时，模型的参数有时是未知的，所以对CSTR的模型参数辨识方法进行研究，具有重要的工程意义。然而，已有的方法如扩展卡尔曼滤波等，这些方法辨识结果有时会出现发散，得不到正确的结果。为了提高辨识效率和精度，研究运用无迹卡尔曼滤波算法进行CSTR模型参数辨识，具有重要的意义。Due to the important role of CSTR in chemical production, it is necessary to study the process in detail. When analyzing the model of CSTR, the parameters of the model are sometimes unknown, so it is of great engineering significance to study the identification method of the model parameters of CSTR. However, the existing methods such as extended Kalman filter, etc., the identification results of these methods sometimes diverge, and the correct results cannot be obtained. In order to improve the identification efficiency and accuracy, it is of great significance to study the use of unscented Kalman filter algorithm for CSTR model parameter identification.

发明内容Contents of the invention

为了有效的了解化工生产中CSTR化学反应系统，本发明提出了一种基于无迹卡尔曼滤波算法的CSTR模型参数辨识方法，有效的实现了CSTR模型的参数辨识。In order to effectively understand the CSTR chemical reaction system in chemical production, the present invention proposes a CSTR model parameter identification method based on an unscented Kalman filter algorithm, which effectively realizes the parameter identification of the CSTR model.

本发明的技术解决方案是：一种基于无迹卡尔曼滤波算法的CSTR模型参数辨识方法，其步骤如下：The technical solution of the present invention is: a kind of CSTR model parameter identification method based on unscented Kalman filter algorithm, its steps are as follows:

(1)、获取扩维状态向量中包含CSTR模型待辨识参数的状态空间表达式；(1) Obtain the state space expression containing the parameters to be identified of the CSTR model in the expanded dimension state vector;

(2)、运用欧拉算法对连续的状态空间表达式进行离散化，获得离散化的状态空间表达式；(2) Using the Euler algorithm to discretize the continuous state space expression to obtain the discretized state space expression;

(3)、初始化，包括：设定参数辨识的初值初始参数辨识误差协方差以及过程噪声和量测噪声所满足的协方差矩阵Q和R，算法迭代次数最大值L；(3) Initialization, including: setting the initial value of parameter identification initial parameter identification error covariance And the covariance matrices Q and R satisfied by the process noise and measurement noise, the maximum number of algorithm iterations L;

(4)、选取k-1时刻的sigma点，计算公式为：(4), select the sigma point at time k-1, the calculation formula is:

${δ δ}_{k k - - 11} = = [\begin{matrix} {\overset{^^}{x x}}_{k k - - 11} & {\overset{^^}{x x}}_{k k - - 11} + + γ γ \cdot \cdot \sqrt{{\overset{^^}{p p}}_{x x,, k k - - 11}} & {\overset{^^}{x x}}_{k k - - 11} - - γ γ \cdot \cdot \sqrt{{\overset{^^}{p p}}_{x x,, k k - - 11}} \end{matrix}]$

$γ γ = = \sqrt{n no + + λ λ},, λ λ = = {α α}^{22} \cdot \cdot ((n no + + {k k}_{f f})) - - n no$

式中，表示k-1时刻的状态估计值，表示k-1时刻的状态估计误差协方差，γ表示尺度参数，n表示的维数；常数α决定sigma点围绕均值的波动范围，通常α∈[10^-4,1]；常数k_f是另一尺度参数，用于状态估计和参数辨识时通常取0。In the formula, Indicates the estimated value of the state at time k-1, Indicates the state estimation error covariance at time k-1, γ indicates the scale parameter, and n indicates The dimension; the constant α determines the sigma point around the mean The fluctuation range of , usually α∈[10 ^-4 ,1]; the constant k _f is another scale parameter, which is usually taken as 0 when used for state estimation and parameter identification.

(5)、在上一步基础上，计算k-1时刻的sigma点的增值点，计算公式为：(5), on the basis of the previous step, calculate the value-added point of the sigma point at k-1 time, the calculation formula is:

${(({f f}_{k k - - 11}^{* *}))}_{i i} = = f f (({(({δ δ}_{k k - - 11}))}_{i i},, {u u}_{k k - - 11}))$

${(({h h}_{k k - - 11}^{* *}))}_{i i} = = h h (({(({δ δ}_{k k - - 11}))}_{i i},, {u u}_{k k - - 11}))$

式中，f(·)是对应具体问题系统方程的非线性函数，h(·)是对应具体问题输出方程中的非线性函数，u_k-1是k-1时刻输入控制矩阵，下标i表示对应于第i个sigma点的相关取值，i＝0…2n。In the formula, f(·) is the nonlinear function corresponding to the system equation of the specific problem, h(·) is the nonlinear function in the output equation corresponding to the specific problem, u _k-1 is the input control matrix at time k-1, and the subscript i Indicates the correlation value corresponding to the i-th sigma point, i=0...2n.

(6)、计算k-1时刻的状态向量均值和协方差，计算公式为：(6) Calculate the mean value and covariance of the state vector at time k-1, the calculation formula is:

${\overset{^^}{x x}}_{k k - - 11}^{- -} = = {Σ Σ}_{i i = = 00}^{22 n no} {w w}_{i i}^{m m} {(({f f}_{k k - - 11}^{* *}))}_{i i}$

${\overset{^^}{p p}}_{x x,, k k - - 11}^{- -} = = {Σ Σ}_{i i = = 00}^{22 n no} {w w}_{i i}^{c c} (({(({f f}_{k k - - 11}^{* *}))}_{i i} - - {\overset{^^}{x x}}_{k k - - 11}^{- -})) {(({(({f f}_{k k - - 11}^{* *}))}_{i i} - - {\overset{^^}{x x}}_{k k - - 11}^{- -}))}^{T T} + + {Q Q}_{k k - - 11}$

式中，表示k-1时刻的状态向量均值，表示k-1时刻的状态协方差，Q_k-1表示k-1时刻系统噪声所满足的协方差矩阵，权重系数和取值的计算公式如下：In the formula, Indicates the mean value of the state vector at time k-1, Represents the state covariance at time k-1, Q _k-1 represents the covariance matrix satisfied by the system noise at time k-1, and the weight coefficient and The formula for calculating the value is as follows:

${w w}_{00}^{m m} = = \frac{λ λ}{n no + + λ λ},, {w w}_{00}^{c c} = = \frac{λ λ}{n no + + λ λ} + + ((11 - - {α α}^{22} + + β β))$

${w w}_{i i}^{m m} = = {w w}_{i i}^{c c} = = \frac{11}{22 ((n no + + λ λ))},, i i = = 11 ... ... 22 n no$

式中β通常是包含x分布的先验知识，对高斯分布来说，其最优值一般取2。In the formula, β usually contains prior knowledge of x distribution, and for Gaussian distribution, its optimal value is generally 2.

(7)、计算k-1时刻量测向量均值和协方差，计算公式为：(7) Calculate the mean value and covariance of the measurement vector at time k-1, the calculation formula is:

${y the y}_{k k - - 11}^{- -} = = {Σ Σ}_{i i = = 00}^{22 n no} {w w}_{i i}^{m m} {(({h h}_{k k - - 11}^{* *}))}_{i i}$

${\overset{^^}{p p}}_{y the y,, k k - - 11}^{- -} = = {Σ Σ}_{i i = = 00}^{22 n no} {w w}_{i i}^{c c} (({(({h h}_{k k - - 11}^{* *}))}_{i i} - - {y the y}_{k k - - 11}^{- -})) {(({(({h h}_{k k - - 11}^{* *}))}_{i i} - - {y the y}_{k k - - 11}^{- -}))}^{T T} + + {R R}_{k k - - 11}$

式中表示k-1时刻量测向量均值，表示k-1量测向量的协方差，R_k-1表示k-1时刻的量测噪声所满足的协方差矩阵。In the formula Indicates the mean value of the measured vector at time k-1, Represents the covariance of k-1 measurement vectors, and R _k-1 represents the covariance matrix satisfied by the measurement noise at time k-1.

(8)、计算交互协方差，计算公式如下：(8) Calculate the interaction covariance, the calculation formula is as follows:

${\overset{^^}{p p}}_{x x y the y,, k k - - 11}^{- -} = = {Σ Σ}_{i i = = 00}^{22 n no} {w w}_{i i}^{c c} (({(({f f}_{k k - - 11}^{* *}))}_{i i} - - {\overset{^^}{x x}}_{k k - - 11}^{- -})) {(({(({h h}_{k k - - 11}^{* *}))}_{i i} - - {y the y}_{k k - - 11}^{- -}))}^{T T}$

式中，表示k-1时刻的交互协方差。In the formula, Indicates the interaction covariance at time k-1.

(9)、在上一步的基础上，计算k-1时刻的卡尔曼滤波增益，其遵循的计算公式为：(9), on the basis of the previous step, calculate the Kalman filter gain at k-1 moment, the calculation formula that it follows is:

${K K}_{k k - - 11} = = {\overset{^^}{p p}}_{x x y the y,, k k - - 11}^{- -} \cdot &Center Dot; {(({\overset{^^}{p p}}_{y the y,, k k - - 11}^{- -}))}^{- - 11}$

式中，K_k-1表示k-1时刻的卡尔曼滤波增益。In the formula, K _k-1 represents the Kalman filter gain at time k-1.

(10)、运用无迹卡尔曼滤波更新步，获得k时刻的状态估计值和协方差，计算公式为：(10), use the unscented Kalman filter update step to obtain the state estimate and covariance at time k, the calculation formula is:

${\overset{^^}{x x}}_{k k} = = {\overset{^^}{x x}}_{k k - - 11}^{- -} + + {K K}_{k k - - 11} \cdot &Center Dot; (({y the y}_{k k - - 11} - - {y the y}_{k k - - 11}^{- -}))$

${\overset{^^}{p p}}_{x x,, k k} = = {\overset{^^}{p p}}_{x x,, k k - - 11}^{- -} - - {K K}_{k k - - 11} \cdot &Center Dot; {\overset{^^}{p p}}_{y the y,, k k - - 11}^{- -} \cdot &Center Dot; {K K}_{k k - - 11}^{T T}$

式中，表示k时刻的状态估计值，y_k-1表示k-1时刻量测输出真实值，表示k时刻的估计协方差。In the formula, Indicates the estimated value of the state at time k, y _k-1 indicates the true value of the measured output at time k-1, Denotes the estimated covariance at time k.

(11)、依据上述步骤进行多次迭代辨识，直至k≥L时，结束迭代过程，输出辨识结果。(11) Carry out multiple iterative identification according to the above steps until k≥L, end the iterative process, and output the identification result.

本发明的有益效果是：运用无迹卡尔曼滤波对CSTR模型参数辨识，由于无迹卡尔曼滤波采用了不同于扩展卡尔曼滤波的线性化方法，因此，克服了扩展卡尔曼滤波辨识过程中出现发散的情况，提高了辨识效率和精度。The beneficial effect of the present invention is: use the unscented Kalman filter to identify the CSTR model parameters, because the unscented Kalman filter adopts a linearization method different from the extended Kalman filter, therefore, overcomes the occurrence of the extended Kalman filter identification process. In the case of divergence, the identification efficiency and accuracy are improved.

附图说明Description of drawings

图1为本发明实施例的方法流程图；Fig. 1 is the method flowchart of the embodiment of the present invention;

图2为实施例采用本发明所提方法的辨识结果；Fig. 2 is the identification result of the embodiment adopting the method proposed by the present invention;

图3为实施例辨识结果的相对误差；Fig. 3 is the relative error of embodiment identification result;

具体实施方式detailed description

下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。Below in conjunction with specific embodiment, further illustrate the present invention, should be understood that these embodiments are only used to illustrate the present invention and are not intended to limit the scope of the present invention, after having read the present invention, those skilled in the art will understand various equivalent forms of the present invention All modifications fall within the scope defined by the appended claims of the present application.

如图1所示，CSTR模型参数辨识方法，其包含如下步骤：As shown in Figure 1, the CSTR model parameter identification method includes the following steps:

(1)、获取状态分量中包含CSTR模型待辨识参数的状态空间表达式。(1) Obtain the state space expression including the parameters to be identified of the CSTR model in the state component.

(2)、运用欧拉算法对连续的状态空间表达式进行离散化，获得离散化的状态空间表达式。(2) Using the Euler algorithm to discretize the continuous state space expression to obtain the discretized state space expression.

(3)、初始化，包括：设定参数辨识的初值初始参数辨识误差协方差以及过程噪声和量测噪声所满足的协方差矩阵Q和R，算法迭代次数最大值L。(3) Initialization, including: setting the initial value of parameter identification initial parameter identification error covariance And the covariance matrices Q and R satisfied by the process noise and measurement noise, and the maximum number of algorithm iterations L.

(4)、选取k-1时刻的sigma点。(4) Select the sigma point at time k-1.

(5)、在上一步基础上，计算k-1时刻的sigma点的增值点。(5) On the basis of the previous step, calculate the value-added point of the sigma point at time k-1.

(6)、计算k-1时刻的状态向量均值和协方差。(6) Calculate the mean value and covariance of the state vector at time k-1.

(7)、计算k-1时刻量测向量均值和协方差。(7) Calculate the mean value and covariance of the measurement vector at time k-1.

(8)、计算k-1时刻的交互协方差。(8) Calculate the interaction covariance at time k-1.

(9)、在上一步的基础上，计算k-1时刻的卡尔曼滤波增益。(9) On the basis of the previous step, calculate the Kalman filter gain at time k-1.

(10)、运用无迹卡尔曼滤波更新步，获得k时刻的状态估计值和协方差。(10), using the unscented Kalman filter update step to obtain the state estimation value and covariance at time k.

假设在连续搅拌釜式反应器(CSTR)内发生放热的、不可逆反应。反应物为A、生成物为B。根据物料平衡和能量平衡关系，得到如下的反应过程模型(CSTR状态空间描述)：An exothermic, irreversible reaction is assumed to occur in a continuous stirred tank reactor (CSTR). The reactant is A and the product is B. According to the material balance and energy balance relationship, the following reaction process model (CSTR state space description) is obtained:

$\{\begin{matrix} {\overset{\cdot &Center Dot;}{C C}}_{A A} = = \frac{q q}{V V} (({C C}_{A A f f} - - {C C}_{A A})) - - {k k}_{00} {e e}^{- - \frac{E E.}{R R T T}} {C C}_{A A} \\ \overset{\cdot \cdot}{T T} = = \frac{q q}{V V} (({T T}_{f f} - - T T)) + + \frac{((- - Δ Δ H h))}{{ρC ρC}_{p p}} {k k}_{00} {e e}^{- - \frac{E E.}{R R T T}} {C C}_{A A} + + \frac{U u A A}{{VρC VρC}_{p p}} (({T T}_{c c} - - T T)) \end{matrix}$

该非线性系统具有两个状态变量即：反应器中组分A的浓度C_A(状态变量x₁)，反应温度T(状态变量x₂)；一个控制输入变量T_c(输入u)。在下面的实施例仿真分析时假设ρ和是待辨识的参数，它们的真值分别是1000和8750。模型中其他参数的物理意义及取值如表1。The nonlinear system has two state variables: the concentration CA of component _A in the reactor (state variable x ₁ ), the reaction temperature T (state variable x ₂ ); and a control input variable T _c (input u). Assume that ρ and are the parameters to be identified, and their true values are 1000 and 8750, respectively. The physical meanings and values of other parameters in the model are shown in Table 1.

表1 CSTR模型中相关参数的物理意义及取值Table 1 The physical meaning and value of relevant parameters in the CSTR model

为了运用无迹卡尔曼滤波对模型中的未知参数进行辨识，首先需要获取状态分量中包含待辨识参数的状态空间表达式，为此假设ρ为状态分量x₃、为状态分量x₄，则整理后的扩维状态空间表达式为(带入相关的模型参数)：In order to use the unscented Kalman filter to identify the unknown parameters in the model, it is first necessary to obtain the state space expression containing the parameters to be identified in the state components. For this reason, it is assumed that ρ is the state components x ₃ , is the state component x ₄ , then the expanded dimension state space expression after sorting out is (introducing relevant model parameters):

状态方程为： The state equation is:

输出方程为： The output equation is:

式中w_i(i＝1,2,3,4)是系统噪声，v_i(i＝1,2)是量测噪声，它们均为零均值的高斯白噪声，且分别满足协方差矩阵Q_k和R_k，即：In the formula, w _i (i=1,2,3,4) is the system noise, v _i (i=1,2) is the measurement noise, they are Gaussian white noise with zero mean, and respectively satisfy the covariance matrix Q _k and R _k , namely:

${Q Q}_{k k} = = E E. [[{w w}_{k k} {w w}_{k k}^{T T}]],, {R R}_{k k} = = E E. [[{v v}_{k k} {v v}_{k k}^{T T}]]$

采用欧拉算法对上述状态方程进行离散化，得到了对应的离散状态方程，在此基础上，则就可以运用无迹卡尔曼滤波进行多次迭代辨识。在辨识过程中，无迹卡尔曼滤波的相关参数取值为：The Euler algorithm is used to discretize the above state equation, and the corresponding discrete state equation is obtained. On this basis, the unscented Kalman filter can be used for multiple iterative identification. In the identification process, the relevant parameters of the unscented Kalman filter are:

$Q Q = = [\begin{matrix} 1010^{- - 55} & 00 & 00 & 00 \\ 00 & 1010^{- - 33} & 00 & 00 \\ 00 & 00 & 1010^{- - 55} & 00 \\ 00 & 00 & 00 & 1010^{- - 55} \end{matrix}],, R R = = [\begin{matrix} 1010^{- - 44} & 00 \\ 00 & 1010^{- - 22} \end{matrix}],, u u = = 290290$

基于上述分析，通过多次迭代辨识，获得了准确的辨识结果。图1为实施例所用的算法流程图，图2为实施例的参数辨识结果，图3为运用本发明所提出的方法对实施例辨识结果的相对误差。Based on the above analysis, through multiple iterations of identification, accurate identification results are obtained. Fig. 1 is a flow chart of the algorithm used in the embodiment, Fig. 2 is the parameter identification result of the embodiment, and Fig. 3 is the relative error of the identification result of the embodiment using the method proposed by the present invention.

Claims

1. A CSTR model parameter identification method based on unscented Kalman filter algorithm, is characterized in that, comprises the steps:

(1) Obtain the state space expression including the parameters to be identified of the CSTR model in the state component.

(2) Using the Euler algorithm to discretize the continuous state space expression to obtain a discrete state space expression.

(3) Initialization, including: setting the initial value of parameter identification initial parameter identification error covariance And the covariance matrices Q and R satisfied by the process noise and measurement noise, and the maximum number of algorithm iterations L.

(4), select the sigma point at time k-1, the calculation formula is:

{δ δ}_{k k - - 11} = = [\begin{matrix} {\overset{^^}{x x}}_{k k - - 11} & {\overset{^^}{x x}}_{k k - - 11} + + γ γ \cdot \cdot \sqrt{{\overset{^^}{p p}}_{x x,, k k - - 11}} & {\overset{^^}{x x}}_{k k - - 11} - - γ γ \cdot &Center Dot; \sqrt{{\overset{^^}{p p}}_{x x,, k k - - 11}} \end{matrix}]

γ γ = = \sqrt{n no + + λ λ},, λ λ = = {α α}^{22} \cdot &Center Dot; ((n no + + {k k}_{f f})) - - n no

In the formula, Indicates the estimated value of the state at time k-1, Indicates the state estimation error covariance at time k-1, γ indicates the scale parameter, and n indicates The dimension; the constant α determines the sigma point around the mean The fluctuation range of , usually α∈[10 ^-4 ,1]; the constant k _f is another scale parameter, which is usually taken as 0 when used for state estimation and parameter identification.

(5), on the basis of the previous step, calculate the value-added point of the sigma point at k-1 time, the calculation formula is:

{(({f f}_{k k - - 11}^{* *}))}_{i i} = = f f (({(({δ δ}_{k k - - 11}))}_{i i},, {u u}_{k k - - 11}))

{(({h h}_{k k - - 11}^{* *}))}_{i i} = = h h (({(({δ δ}_{k k - - 11}))}_{i i},, {u u}_{k k - - 11}))

In the formula, f(·) is the nonlinear function corresponding to the system equation of the specific problem, h(·) is the nonlinear function in the output equation corresponding to the specific problem, u _k-1 is the input control matrix at time k-1, and the subscript i Indicates the correlation value corresponding to the i-th sigma point, i=0...2n.

(6), calculate the mean value and covariance of the state vector at k-1 time, the calculation formula is:

{\overset{^^}{x x}}_{k k - - 11}^{- -} = = {Σ Σ}_{i i = = 00}^{22 n no} {w w}_{i i}^{m m} {(({f f}_{k k - - 11}^{* *}))}_{i i}

{\overset{^^}{p p}}_{x x,, k k - - 11}^{- -} = = {Σ Σ}_{i i = = 00}^{22 n no} {w w}_{i i}^{c c} (({(({f f}_{k k - - 11}^{* *}))}_{i i} - - {\overset{^^}{x x}}_{k k - - 11}^{- -})) {(({(({f f}_{k k - - 11}^{* *}))}_{i i} - - {\overset{^^}{x x}}_{k k - - 11}^{- -}))}^{T T} + + {Q Q}_{k k - - 11}

In the formula, Indicates the mean value of the state vector at time k-1, Represents the state covariance at time k-1, Q _k-1 represents the covariance matrix satisfied by the system noise at time k-1, and the weight coefficient and The formula for calculating the value is as follows:

{w w}_{00}^{m m} = = \frac{λ λ}{n no + + λ λ},, {w w}_{00}^{c c} = = \frac{λ λ}{n no + + λ λ} + + ((11 - - {α α}^{22} + + β β))

{w w}_{i i}^{m m} = = {w w}_{i i}^{c c} = = \frac{11}{22 ((n no + + λ λ))},, i i = = 11 ... ... 22 n no

In the formula, β usually contains prior knowledge of x distribution, and for Gaussian distribution, its optimal value is generally 2.

(7) Calculate the mean value and covariance of the measurement vector at time k-1, the calculation formula is:

{y the y}_{k k - - 11}^{- -} = = {Σ Σ}_{i i = = 00}^{22 n no} {w w}_{i i}^{m m} {(({h h}_{k k - - 11}^{* *}))}_{i i}

{\overset{^^}{p p}}_{y the y,, k k - - 11}^{- -} = = {Σ Σ}_{i i = = 00}^{22 n no} {w w}_{i i}^{c c} (({(({h h}_{k k - - 11}^{* *}))}_{i i} - - {y the y}_{k k - - 11}^{- -})) {(({(({h h}_{k k - - 11}^{* *}))}_{i i} - - {y the y}_{k k - - 11}^{- -}))}^{T T} + + {R R}_{k k - - 11}

In the formula Indicates the mean value of the measured vector at time k-1, Represents the covariance of k-1 measurement vectors, and R _k-1 represents the covariance matrix satisfied by the measurement noise at time k-1.

(8) Calculate the interaction covariance, the calculation formula is as follows:

{\overset{^^}{p p}}_{x x y the y,, k k - - 11}^{- -} = = {Σ Σ}_{i i = = 00}^{22 n no} {w w}_{i i}^{c c} (({(({f f}_{k k - - 11}^{* *}))}_{i i} - - {\overset{^^}{x x}}_{k k - - 11}^{- -})) {(({(({h h}_{k k - - 11}^{* *}))}_{i i} - - {y the y}_{k k - - 11}^{- -}))}^{T T}

In the formula, Indicates the interaction covariance at time k-1.

(9), on the basis of the previous step, calculate the Kalman filter gain at k-1 moment, the calculation formula that it follows is:

{K K}_{k k - - 11} = = {\overset{^^}{p p}}_{x x y the y,, k k - - 11}^{- -} \cdot &Center Dot; {(({\overset{^^}{p p}}_{y the y,, k k - - 11}^{- -}))}^{- - 11}

In the formula, K _k-1 represents the Kalman filter gain at time k-1.

(10), use the unscented Kalman filter update step to obtain the state estimate and covariance at time k, the calculation formula is:

{\overset{^^}{x x}}_{k k} = = {\overset{^^}{x x}}_{k k - - 11}^{- -} + + {K K}_{k k - - 11} \cdot &Center Dot; (({y the y}_{k k - - 11} - - {y the y}_{k k - - 11}^{- -}))

{\overset{^^}{p p}}_{x x,, k k} = = {\overset{^^}{p p}}_{x x,, k k - - 11}^{- -} - - {K K}_{k k - - 11} \cdot &Center Dot; {\overset{^^}{p p}}_{y the y,, k k - - 11}^{- -} \cdot &Center Dot; {K K}_{k k - - 11}^{T T}

In the formula, Indicates the estimated value of the state at time k, y _k-1 indicates the true value of the measured output at time k-1, Denotes the estimated covariance at time k.

(11) Carry out multiple iterative identification according to the above steps until k≥L, end the iterative process, and output the identification result.