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

Abstract

The invention discloses a CSTR (Continuous Stirred Tank Reactor) model parameter identification method based on an unscented Kalman filtering algorithm. The method comprises the following steps: according to a CSTR continuous system model, obtaining a state spatial expression of which the state component contains a parameter to be identified; then, in virtue of an Euler algorithm, carrying out discretization processing on the obtained non-linear continuous state spatial expression to obtain a corresponding discrete iteration model; and finally, applying the unscented Kalman filtering algorithm to carry out multi-time iteration identification, and obtaining an accurate identification result. The algorithm has good astringency, is easy in combining with traditional software and has a good engineering application prospect.

Description

A kind of CSTR model parameter identification method based on Unscented kalman filtering algorithm

Technical field

The present invention relates to a kind of CSTR model parameter identification method based on Unscented kalman filtering algorithm, belong to System modelling and parameter identification technique field.

Background technology

CSTR (CSTR) is typical, the chemistry of nonlinearity in chemical process Response system.In the nucleus equipment of Chemical Manufacture, occupy considerable status, dyestuff, pharmaceutical reagent, In food and synthetic material industry, CSTR system is widely used.

Due to CSTR important function in Chemical Manufacture, it is therefore necessary to this process is carried out detailed grinding Study carefully.When CSTR setting up model and being analyzed, the parameter of model is unknown sometimes, so to CSTR Model parameter identification method study, there is important engineering significance.But, existing method such as extension Kalman filterings etc., these method identification results, it sometimes appear that dissipate, can not get correct result.In order to carry High identification efficiency and precision, research application Unscented kalman filtering algorithm carries out CSTR identification of Model Parameters, tool There is important meaning.

Summary of the invention

In order to effectively understand CSTR chemical reaction system in Chemical Manufacture, the present invention proposes a kind of based on nothing The CSTR model parameter identification method of mark Kalman filtering algorithm, effectively achieves the parameter of CSTR model Identification.

The technical solution of the present invention is: a kind of CSTR model parameter based on Unscented kalman filtering algorithm Discrimination method, its step is as follows:

(1), acquisition is augmented in state vector the state-space expression comprising CSTR model parameter to be identified；

(2), use euler algorithm that continuous print state-space expression is carried out discretization, it is thus achieved that the shape of discretization State space expression formula；

(3), initialize, including: the initial value of setup parameter identificationInitial parameter Identification Errors covariance And covariance matrix Q and R that process noise and measurement noise are met, algorithm iteration number of times maximum L；

(4), choosing the sigma point in k-1 moment, computing formula is:

${\mathrm{\δ}}_{k-1}=\left[\begin{array}{ccc}{\hat{x}}_{k-1}& {\hat{x}}_{k-1}+\mathrm{\γ}\·\sqrt{{\hat{p}}_{x,k-1}}& {\hat{x}}_{k-1}-\mathrm{\γ}\·\sqrt{{\hat{p}}_{x,k-1}}\end{array}\right]$

$\mathrm{\γ}=\sqrt{n+\mathrm{\λ}},\mathrm{\λ}={\mathrm{\α}}^2\·(n+k_f)-n$

In formula,Represent the state estimation in k-1 moment,Represent the state estimation error association side in k-1 moment Difference, γ represents that scale parameter, n representDimension；Constant α determines that sigma point is around averageRipple Dynamic scope, usual α ∈ [10^-4,1]；Constant k_fIt is another scale parameter, when state estimation and parameter identification Generally take 0.

(5), on the basis of previous step, calculating the point of increase of the sigma point in k-1 moment, computing formula is:

${\left(f_{k-1}^{*}\right)}_i=f({\left({\mathrm{\δ}}_{k-1}\right)}_i,u_{k-1})$

${\left(h_{k-1}^{*}\right)}_i=h({\left({\mathrm{\δ}}_{k-1}\right)}_i,u_{k-1})$

In formula, f () is the nonlinear function of corresponding particular problem system equation, and h () is corresponding particular problem output side Nonlinear function in journey, u_k-1Being to input control matrix the k-1 moment, subscript i represents corresponding to i-th sigma The relevant value of point, i=0 ... 2n.

(6), calculating the state vector average in k-1 moment and covariance, computing formula is:

${\hat{x}}_{k-1}^{-}=\underset{i=0}{\overset{2n}{\Σ}}{w}_i^m{\left(f_{k-1}^{*}\right)}_i$

${\hat{p}}_{x,k-1}^{-}=\underset{i=0}{\overset{2n}{\Σ}}{w}_i^c({\left(f_{k-1}^{*}\right)}_i-{\hat{x}}_{k-1}^{-}){({\left(f_{k-1}^{*}\right)}_i-{\hat{x}}_{k-1}^{-})}^T+Q_{k-1}$

In formula,Represent the state vector average in k-1 moment,Represent the state covariance in k-1 moment, Q_k-1 Represent the covariance matrix that k-1 moment system noise is met, weight coefficientWithThe computing formula of value As follows:

$w_0^m=\frac{\mathrm{\λ}}{n+\mathrm{\λ}},w_0^c=\frac{\mathrm{\λ}}{n+\mathrm{\λ}}+(1-{\mathrm{\α}}^2+\mathrm{\β})$

$w_i^m=w_i^c=\frac{1}{2(n+\mathrm{\λ})},i=1...2n$

In formula, β is typically to comprise the priori of x distribution, and for Gauss distribution, its optimal value typically takes 2.

(7), calculating the k-1 moment and measures vectorial average and covariance, computing formula is:

$y_{k-1}^{-}=\underset{i=0}{\overset{2n}{\Σ}}{w}_i^m{\left(h_{k-1}^{*}\right)}_i$

${\hat{p}}_{y,k-1}^{-}=\underset{i=0}{\overset{2n}{\Σ}}{w}_i^c({\left(h_{k-1}^{*}\right)}_i-y_{k-1}^{-}){({\left(h_{k-1}^{*}\right)}_i-y_{k-1}^{-})}^T+R_{k-1}$

In formulaRepresent that the k-1 moment measures vector average,Represent that k-1 measures the covariance of vector, R_k-1Table Show the covariance matrix that the measurement noise in k-1 moment is met.

(8), calculating cross-covariance, computing formula is as follows:

${\hat{p}}_{xy,k-1}^{-}=\underset{i=0}{\overset{2n}{\Σ}}{w}_i^c({\left(f_{k-1}^{*}\right)}_i-{\hat{x}}_{k-1}^{-}){({\left(h_{k-1}^{*}\right)}_i-y_{k-1}^{-})}^T$

In formula,Represent the cross-covariance in k-1 moment.

(9), on the basis of previous step, the Kalman filtering gain in k-1 moment, its calculating followed are calculated Formula is:

$K_{k-1}={\hat{p}}_{xy,k-1}^{-}\·{\left({\hat{p}}_{y,k-1}^{-}\right)}^{-1}$

In formula, K_k-1Represent the Kalman filtering gain in k-1 moment.

(10), Unscented kalman filtering is used to update step, it is thus achieved that the state estimation in k moment and covariance, Computing formula is:

${\hat{x}}_k={\hat{x}}_{k-1}^{-}+K_{k-1}\·(y_{k-1}-y_{k-1}^{-})$

${\hat{p}}_{x,k}={\hat{p}}_{x,k-1}^{-}-K_{k-1}\·{\hat{p}}_{y,k-1}^{-}\·K_{k-1}^T$

In formula,Represent the state estimation in k moment, y_k-1Represent that the k-1 moment measures output actual value,Represent The estimate covariance in k moment.

(11), successive ignition identification is carried out according to above-mentioned steps, until during k >=L, terminating iterative process, defeated Go out identification result.

The invention has the beneficial effects as follows: utilization Unscented kalman filtering is to CSTR identification of Model Parameters, due to nothing Mark Kalman filtering have employed the linearization technique being different from EKF, therefore, overcomes expansion card Kalman Filtering identification process occurs situation about dissipating, improves identification efficiency and precision.

Accompanying drawing explanation

Fig. 1 is the method flow diagram of the embodiment of the present invention；

Fig. 2 is the identification result that embodiment uses institute of the present invention extracting method；

Fig. 3 is the relative error of embodiment identification result；

Detailed description of the invention

Below in conjunction with specific embodiment, it is further elucidated with the present invention, it should be understood that these embodiments are merely to illustrate this Invention rather than restriction the scope of the present invention, after having read the present invention, those skilled in the art are to this The amendment of the bright various equivalent form of values all falls within the application claims limited range.

As it is shown in figure 1, CSTR model parameter identification method, it comprises the steps of:

(1) state-space expression comprising CSTR model parameter to be identified in state component, is obtained.

(2), use euler algorithm that continuous print state-space expression is carried out discretization, it is thus achieved that the shape of discretization State space expression formula.

(3), initialize, including: the initial value of setup parameter identificationInitial parameter Identification Errors covariance And covariance matrix Q and R that process noise and measurement noise are met, algorithm iteration number of times maximum L.

(4) the sigma point in k-1 moment, is chosen.

(5), on the basis of previous step, the point of increase of the sigma point in k-1 moment is calculated.

(6) state vector average and the covariance in k-1 moment, are calculated.

(7), the calculating k-1 moment measures vector average and covariance.

(8) cross-covariance in k-1 moment, is calculated.

(9), on the basis of previous step, the Kalman filtering gain in k-1 moment is calculated.

(10), Unscented kalman filtering is used to update step, it is thus achieved that the state estimation in k moment and covariance.

(11), successive ignition identification is carried out according to above-mentioned steps, until during k >=L, terminating iterative process, defeated Go out identification result.

Assume to occur heat release, irreversible reaction in CSTR (CSTR).Reactant be A, Product is B.According to material balance and energy balance relations, obtain following reaction procedure model (CSTR shape State space describes):

$\left\{\begin{array}{c}{\stackrel{\·}{C}}_A=\frac{q}{V}(C_{Af}-C_A)-k_0{e}^{-\frac{E}{RT}}{C}_A\\ \stackrel{\·}{T}=\frac{q}{V}(T_f-T)+\frac{(-\mathrm{\Δ}H)}{{\mathrm{\ρC}}_p}{k}_0{e}^{-\frac{E}{RT}}{C}_A+\frac{UA}{{\mathrm{V\ρC}}_p}(T_c-T)\end{array}\right.$

This nonlinear system has two state variables i.e.: the concentration C of component A in reactor_A(state variable x₁), instead Answer temperature T (state variable x₂)；One controls input variable T_c(input u).At following embodiment simulation analysis Time assume ρ andBeing parameter to be identified, their true value is 1000 and 8750 respectively.Other ginsengs in model The physical significance of number and value such as table 1.

The physical significance of relevant parameter and value in table 1 CSTR model

In order to use Unscented kalman filtering that the unknown parameter in model is carried out identification, it is necessary first to obtain state Component comprises the state-space expression of parameter to be identified, assumes that for this ρ is state component x₃、For state Component x₄, then the state-space expression that is augmented after arranging is (bringing relevant model parameter into):

State equation is:

Output equation is:

W in formula_i(i=1,2,3,4) is system noise, v_i(i=1,2) is measurement noise, and they are the white Gaussian of zero-mean Noise, and meet covariance matrix Q respectively_kAnd R_k, it may be assumed that

$Q_k=E\[w_k{w}_k^T\],R_k=E\[v_k{v}_k^T\]$

Use euler algorithm that above-mentioned state equation is carried out discretization, obtain the discrete state equations of correspondence, at this base On plinth, Unscented kalman filtering the most just can be used to carry out successive ignition identification.In identification process, without mark card The relevant parameter value of Kalman Filtering is:

$Q=\left[\begin{array}{cccc}{10}^{-5}& 0& 0& 0\\ 0& {10}^{-3}& 0& 0\\ 0& 0& {10}^{-5}& 0\\ 0& 0& 0& {10}^{-5}\end{array}\right],R=\left[\begin{array}{cc}{10}^{-4}& 0\\ 0& {10}^{-2}\end{array}\right],u=290$

Based on above-mentioned analysis, by successive ignition identification, it is thus achieved that identification result accurately.Fig. 1 is embodiment Algorithm flow chart used, Fig. 2 is the parameter identification result of embodiment, and Fig. 3 is that utilization is proposed by the invention The method relative error to embodiment identification result.

Claims (1)

1. a CSTR model parameter identification method based on Unscented kalman filtering algorithm, it is characterised in that Comprise the steps of:

(1) state-space expression comprising CSTR model parameter to be identified in state component, is obtained.

(2), use euler algorithm that continuous print state-space expression is carried out discretization, it is thus achieved that discrete state Spatial expression.

(3), initialize, including: the initial value of setup parameter identificationInitial parameter Identification Errors covariance And covariance matrix Q and R that process noise and measurement noise are met, algorithm iteration number of times maximum L.

(4), choosing the sigma point in k-1 moment, computing formula is:

${\mathrm{\δ}}_{k-1}=\left[\begin{array}{ccc}{\hat{x}}_{k-1}& {\hat{x}}_{k-1}+\mathrm{\γ}\·\sqrt{{\hat{p}}_{x,k-1}}& {\hat{x}}_{k-1}-\mathrm{\γ}\·\sqrt{{\hat{p}}_{x,k-1}}\end{array}\right]$

$\mathrm{\γ}=\sqrt{n+\mathrm{\λ}},\mathrm{\λ}={\mathrm{\α}}^2\·(n+k_f)-n$

In formula,Represent the state estimation in k-1 moment,Represent the state estimation error association side in k-1 moment Difference, γ represents that scale parameter, n representDimension；Constant α determines that sigma point is around averageRipple Dynamic scope, usual α ∈ [10^-4,1]；Constant k_fIt is another scale parameter, when state estimation and parameter identification Generally take 0.

(5), on the basis of previous step, calculating the point of increase of the sigma point in k-1 moment, computing formula is:

${\left(f_{k-1}^{*}\right)}_i=f({\left({\mathrm{\δ}}_{k-1}\right)}_i,u_{k-1})$

${\left(h_{k-1}^{*}\right)}_i=h({\left({\mathrm{\δ}}_{k-1}\right)}_i,u_{k-1})$

In formula, f () is the nonlinear function of corresponding particular problem system equation, and h () is corresponding particular problem output side Nonlinear function in journey, u_k-1Being to input control matrix the k-1 moment, subscript i represents corresponding to i-th sigma The relevant value of point, i=0 ... 2n.

(6), calculating the state vector average in k-1 moment and covariance, computing formula is:

${\hat{x}}_{k-1}^{-}=\underset{i=0}{\overset{2n}{\Σ}}{w}_i^m{\left(f_{k-1}^{*}\right)}_i$

${\hat{p}}_{x,k-1}^{-}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{w}_i^c({\left(f_{k-1}^{*}\right)}_i-{\hat{x}}_{k-1}^{-}){({\left(f_{k-1}^{*}\right)}_i-{\hat{x}}_{k-1}^{-})}^T+Q_{k-1}$

In formula,Represent the state vector average in k-1 moment,Represent the state covariance in k-1 moment, Q_k-1 Represent the covariance matrix that k-1 moment system noise is met, weight coefficientWithThe computing formula of value As follows:

$w_0^m=\frac{\mathrm{\λ}}{n+\mathrm{\λ}},w_0^c=\frac{\mathrm{\λ}}{n+\mathrm{\λ}}+(1-{\mathrm{\α}}^2+\mathrm{\β})$

$w_i^m=w_i^c=\frac{1}{2(n+\mathrm{\λ})},i=1...2n$

In formula, β is typically to comprise the priori of x distribution, and for Gauss distribution, its optimal value typically takes 2.

(7), calculating the k-1 moment and measures vectorial average and covariance, computing formula is:

$y_{k-1}^{-}=\underset{i=0}{\overset{2n}{\Σ}}{w}_i^m{\left(h_{k-1}^{*}\right)}_i$

${\hat{p}}_{y,k-1}^{-}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{w}_i^c({\left(h_{k-1}^{*}\right)}_i-y_{k-1}^{-}){({\left(h_{k-1}^{*}\right)}_i-y_{k-1}^{-})}^T+R_{k-1}$

In formulaRepresent that the k-1 moment measures vector average,Represent that k-1 measures the covariance of vector, R_k-1Table Show the covariance matrix that the measurement noise in k-1 moment is met.

(8), calculating cross-covariance, computing formula is as follows:

${\hat{p}}_{xy,k-1}^{-}=\underset{i=0}{\overset{2n}{\Σ}}{w}_i^c({\left(f_{k-1}^{*}\right)}_i-{\hat{x}}_{k-1}^{-}){({\left(h_{k-1}^{*}\right)}_i-y_{k-1}^{-})}^T$

In formula,Represent the cross-covariance in k-1 moment.

(9), on the basis of previous step, the Kalman filtering gain in k-1 moment, its calculating followed are calculated Formula is:

$K_{k-1}={\hat{p}}_{xy,k-1}^{-}\·{\left({\hat{p}}_{y,k-1}^{-}\right)}^{-1}$

In formula, K_k-1Represent the Kalman filtering gain in k-1 moment.

(10), Unscented kalman filtering is used to update step, it is thus achieved that the state estimation in k moment and covariance, Computing formula is:

${\hat{x}}_k={\hat{x}}_{k-1}^{-}+K_{k-1}\·(y_{k-1}-y_{k-1}^{-})$

${\hat{p}}_{x,k}={\hat{p}}_{x,k-1}^{-}-K_{k-1}\·{\hat{p}}_{y,k-1}^{-}\·K_{k-1}^T$

In formula,Represent the state estimation in k moment, y_k-1Represent that the k-1 moment measures output actual value,Represent The estimate covariance in k moment.

(11), successive ignition identification is carried out according to above-mentioned steps, until during k >=L, terminating iterative process, defeated Go out identification result.