CN105654053A - Improved constrained EKF algorithm-based dynamic oscillation signal parameter identification method - Google Patents

Abstract

The invention discloses an improved constrained EKF algorithm-based dynamic oscillation signal parameter identification method. According to the method, a state space expression of which the state components include parameters to be identified is constructed according to the mathematic model of dynamic oscillation signals; an improved constrained method and an EKF algorithm are combined, so that an improved constrained EKF algorithm can be designed, and problems in constrained optimization in an constrained EKF framework can be successfully solved by means of an improved PSO algorithm and a penalty function method; and the improved constrained EKF algorithm is utilized to perform a plurality of times of iterative identification, and therefore, dynamic oscillation signal parameter identification with parameter actual constraint conditions considered can be realized. The algorithm considers actual engineering background, and is simple and convenient, and has a certain engineering application value.

Description

Based on the dynamic oscillation signal parameter identification method improving constraint EKF algorithm

Technical field

The present invention relates to a kind of dynamic oscillation signal parameter identification method based on improving constraint EKF algorithm, belong to signal analysis and parameter identification technique field.

Background technology

In recent years, along with modern power network scale constantly expands, the raising day by day of the interconnected degree of electrical network, the dynamic oscillation that system produces after being subject to large and small disturbance has become one of topmost factor of restriction power network safety operation. Owing to these dynamic oscillation signals can provide the important information about Operation of Electric Systems pattern, thus find and accurately grasp these oscillation signal features for power system safety and stability operation significant.

In view of the importance of dynamic oscillation signal recognition, researchist proposes many discrimination methods, such as Matrix Pencil, maximum likelihood method, Pu Longnifa etc. But, these methods mostly have window attribute, cannot realize the real-time online identification of dynamic oscillation signal. Recently, researchist proposes a kind of parameter identification method based on EKF algorithm, and the shortcoming of many algorithms before overcoming effectively achieves the real-time online identification of oscillation signal. But, the method does not consider the physical constraint condition suffered by parameter to be identified, convergency in this way poor, and do not consider the practical background of engineering. Therefore, it is necessary to dynamic oscillation signal parameter identification problem when physical constraint is studied.

Summary of the invention

Goal of the invention: for problems of the prior art, in order to effectively solve real-time identification when dynamic oscillation signal parameter physical constraint, overcome the shortcoming based on traditional E KF algorithm identification, the present invention provides a kind of dynamic oscillation signal parameter identification method based on improving constraint EKF algorithm, effectively achieves the identification under dynamic oscillation signal parameter constraint condition.

Technical scheme: a kind of dynamic oscillation signal parameter identification method based on improving constraint EKF algorithm, the method realizes in a computer successively in accordance with the following steps:

(1) state-space expression comprising dynamic oscillation signal parameter to be identified in state component, is obtained;

(2), initialize, comprising: setting parameter identification initial valueInitial parameter identification error covarianceAnd covariance matrix Q and R that process noise and measurement noise are met, total algorithm iteration number of times maximum value S and population optimizing maximum iteration time M;

(3), by the state estimation in known k-1 moment and state estimation error covariance, utilizing the prediction step of tradition spreading kalman filtering, obtain status predication value and the status predication error covariance in k moment, calculation formula is:

${\tilde{x}}_k=f({\hat{x}}_{k-1},u_{k-1})$

${\tilde{P}}_k=F_{k-1}{\hat{P}}_{k-1}{F}_{k-1}^T+Q_{k-1}$

In formula,Represent the status predication value in k moment, the nonlinear function in the corresponding particular problem state equation of f (),Represent the state estimation vector in k-1 moment, u_k-1Represent the control inputs in k-1 moment.Represent the status predication error covariance in k moment,Represent nonlinear functionPlace refined gram than matrix,Representing the state estimation error covariance in k-1 moment, subscript T represents transposition, Q_k-1It it is the covariance matrix met in the system noise k-1 moment.

(4), on previous step basis, utilizing the filtering of spreading kalman filtering to walk, obtain the state estimation in k moment, calculation procedure is:

$K_k={\tilde{P}}_k{H}_k^T{(H_k{\tilde{P}}_k{H}_k^T+R_k)}^{-1}$

${\hat{P}}_k=(I-K_k{H}_k){\tilde{P}}_k$

${\hat{x}}_k={\tilde{x}}_k+K_k\[y_k-h\left({\tilde{x}}_k\right)\]$

In formula, K_kRepresent the Kalman filtering gain in k moment,Representing the status predication error covariance in k moment, subscript T represents transposition,Represent that nonlinear function h () existsRefined gram of place than matrix, the wherein nonlinear function in the corresponding particular problem output equation of h (). R_kBeing the covariance matrix met in the measurement noise k moment, I is the unit matrix identical with state vector dimension degree,Represent the state estimation vector in k moment, y_kIt it is the work output of k moment output equation.

(5), judge whether the parameter identification result in k moment meets corresponding physical constraint condition. If meeting, then directly use EKF iteration identification again.

(6) if not meeting, then using the bounding algorithm of improvement that this moment parameter identification result carries out constraint adjustment, obtaining constrained optimization objective function is:

${\widehat{\stackrel{\‾}{x}}}_k=\underset{x_k}{\arg \min }{(y_k-{\hat{y}}_k)}^T{W}_1(y_k-{\hat{y}}_k)+{(x_k-{\hat{x}}_k)}^T{W}_2(x_k-{\hat{x}}_k)$

In formula,Represent k moment optimizing state estimation to be asked, y_kIt is k moment work output, W₁And W₂It is the weight matrix of positive definite, for convenience, the embodiment of the present invention is chosen for unit matrix, x in formula_kMeet constraint condition as shown in the formula:

Dx��d

The capable full rank matrix of s �� n constant that D is known, s is the number by constraint condition parameter, and n is the dimension of state vector, it is clear that s��n, d are known constrainted constants conditions.

(7), in step (6) constrained optimization objective function, to the estimation work output in k momentWhen solving, in order to make full use of known measurement data information, the present invention adopts Kalman smoothing method, and calculation formula is as follows:

${\hat{x}}_{s,k-1}={\hat{x}}_{k-1}+L({\hat{x}}_k-{\tilde{x}}_k)$

$L={\hat{P}}_{k-1}{F}_{k-1}{\tilde{P}}_k$

${\hat{y}}_k=h\left(f\right({\hat{x}}_{s,k-1}\left)\right)$

In formula,Representing the level and smooth state estimation vector in k-1 moment, L is Kalman smoothing gain,It it is the estimation work output in k moment.

(8), on the basis of previous step, by penalty function method, by former objective function adds a penalty term, the optimization problem of constraint being converted into a unconfined optimization problem, the calculation formula of its conversion is as follows:

F (x)=f (x)+h (k) H (x), x �� Rⁿ

In formula, f (x) is former objective function, and h (k) is the penalty value dynamically updated, generallyOrK is current iteration number of times. H (x) is the punishment factor, and calculation formula is as follows:

$H\left(x\right)=\underset{i=1}{\overset{m}{\Σ}}\mathrm{\θ}\left(q_i\right(x\left)\right)q_i{\left(x\right)}^{\mathrm{\γ}\left(q_i\right(x\left)\right)}$

In formula, �� (q_i(x)) it is multistage partition function, q_iX () is the function relevant with violating constraint condition, q_i(x)=max{0, g_i(x) }, i=1 ..., s, wherein g_i(x)=Dx-d. �� (q_i(x)) represent penalty function effect. The rule that function value is followed is:

$\mathrm{\γ}\left(q_i\right(x\left)\right)=\left\{\begin{array}{cc}1& q_i\left(x\right)\≤1\\ 2& q_i\left(x\right)>1\end{array}\right.$

$\mathrm{\&theta;}\left(q_i\right(x\left)\right)=\left\{\begin{array}{cc}10& q_i\left(x\right)\&le;0.001\\ 20& 0.001<q_i\left(x\right)\&le;0.1\\ 100& 0.1<q_i\left(x\right)\&le;1\\ 300& q_i\left(x\right)>1\end{array}\right.$

(9), on the basis of previous step, utilize modified particle swarm optiziation to carry out successive ignition optimizing. Wherein, the rule that the population speed of improvement and location updating are followed is as follows:

$V_i^{iter+1}=w\left(L\right)\&times;V_i^{iter}+c_1{r}_{1,i}({\mathrm{Pbest}}_i^{iter}-X_i^{iter})+c_2{r}_{2,i}({\mathrm{Gbest}}^{iter}-X_i^{iter}))$

$X_i^{iter+1}=X_i^{iter}+\mathrm{\&alpha;}\&CenterDot;V_i^{iter+1}$

If search space is D dimension, population comprises N number of particle.In formulaWithRepresent i-th particle respectively at the speed in moment iter and position vector. �� is the contraction factor for controlling and limit speed, c₁And c₂It is respectively cognitive and society's coefficient (generally equal can be set to 2), r_1,iAnd r_2,iIt is two independent randomized numbers of value in [0,1] scope.Represent i-th particle by the iter moment to historical position optimum value. Gbest^iterRepresent in all particles by the historical position optimum value to the iter moment. W (L) is dynamic inertia weight, and for regulating search spatial dimension, optimum to obtaining the overall situation, its calculation formula is:

W (L)=w_start��(w_start-w_end)��(T_max-L)/T_max

In formula, w_startIt is initial inertia weight, w_endBe iteration to inertia weight during maximum times, L is the iteration number of times of the current optimizing of population, T_maxFor the iteration number of times upper limit of population optimizing, in general inertia weights get w_start=0.9, w_endWhen=0.4, population optimizing effect is better.

(10) if iteration number of times L > M, then population optimizing iteration terminates, now the optimizing result Gbest to k moment state estimation^LAs the state estimation in k moment, then subsequent time state is estimated;

(11) if k=k+1��S, then iteration continues, if k=k+1 is > S, then iteration terminates, and exports identification result.

Accompanying drawing explanation

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

The dynamic oscillation signal of Fig. 2 embodiment;

Fig. 3 is the signal frequency w identification result that embodiment adopts institute of the present invention extracting method;

Fig. 4 is the signal damping factor �� identification result that embodiment adopts institute of the present invention extracting method.

Embodiment

Below in conjunction with specific embodiment, illustrate the present invention further, these embodiments should be understood only be not used in for illustration of the present invention and limit the scope of the invention, after having read the present invention, the amendment of the various equivalent form of values of the present invention is all fallen within the application's claims limited range by those skilled in the art.

As shown in Figure 1, dynamic oscillation signal parameter identification method, it comprises following steps:

(1) state-space expression comprising dynamic oscillation signal parameter to be identified in state component, is obtained.

(2), initialize, comprising: setting parameter identification initial valueInitial parameter identification error covarianceAnd covariance matrix Q and R that process noise and measurement noise are met, total algorithm iteration number of times maximum value S and population optimizing maximum iteration time M.

(3), by the state estimation in the k-1 moment obtained and state estimation error covariance, utilizing the prediction of spreading kalman filtering to walk, obtain status predication value and the status predication error covariance in k moment, calculation formula is:

${\tilde{x}}_k=f({\hat{x}}_{k-1},u_{k-1})$

${\tilde{P}}_k=F_{k-1}{\hat{P}}_{k-1}{F}_{k-1}^T+Q_{k-1}$

In formula,Represent the status predication value in k moment, the nonlinear function in the corresponding particular problem state equation of f (),Represent the state estimation vector in k-1 moment, u_k-1Represent the control inputs in k-1 moment.Represent the status predication error covariance in k moment,Represent that nonlinear function f () existsPlace refined gram than matrix,Representing the state estimation error covariance in k-1 moment, subscript T represents transposition, Q_k-1It it is the covariance matrix met in the system noise k-1 moment.

(4), on previous step basis, utilize the filtering of spreading kalman filtering to walk, obtain the state estimation in k moment.

Calculation procedure is:

$K_k={\tilde{P}}_k{H}_k^T{(H_k{\tilde{P}}_k{H}_k^T+R_k)}^{-1}$

${\hat{P}}_k=(I-K_k{H}_k){\tilde{P}}_k$

${\hat{x}}_k={\tilde{x}}_k+K_k\&lsqb;y_k-h\left({\tilde{x}}_k\right)\&rsqb;$

In formula, K_kRepresent the Kalman filtering gain in k moment,Representing the status predication error covariance in k moment, subscript T represents transposition,Represent that nonlinear function h () existsRefined gram of place than matrix, the wherein nonlinear function in the corresponding particular problem output equation of h ().R_kIt is the covariance matrix met in the measurement noise k moment,Representing the state estimation error covariance in k moment, I is the unit matrix identical with state vector dimension degree,Represent the state estimation vector in k moment, y_kIt it is the work output of k moment output equation.

(5), judge whether the parameter identification result in k moment meets corresponding physical constraint condition. If meeting, then directly use EKF iteration identification again.

(6) if not meeting, then using the bounding algorithm of improvement this moment parameter identification result to be retrained, the state estimation problem of constraint equivalence being converted into the optimization problem of constraint.

Constrained optimization objective function is:

${\widehat{\stackrel{\&OverBar;}{x}}}_k=\underset{x_k}{\arg \min }{(y_k-{\hat{y}}_k)}^T{W}_1(y_k-{\hat{y}}_k)+{(x_k-{\hat{x}}_k)}^T{W}_2(x_k-{\hat{x}}_k)$

In formula,Represent k moment optimizing state estimation to be asked, y_kIt is k moment work output, W₁And W₂It is the weight matrix of positive definite, for convenience, the embodiment of the present invention is chosen for unit matrix, x in formula_kMeet constraint condition as shown in the formula:

Dx��d

The capable full rank matrix of s �� n constant that D is known, s is the number by constraint condition parameter, and n is the dimension of state vector, it is clear that s��n, d are known constrainted constants conditions.

In constrained optimization objective function, to the estimation work output in k momentWhen solving, in order to make full use of known measurement data information, adopting Kalman smoothing method, calculation formula is as follows:

${\hat{x}}_{s,k-1}={\hat{x}}_{k-1}+L({\hat{x}}_k-{\tilde{x}}_k)$

$L={\hat{P}}_{k-1}{F}_{k-1}{\tilde{P}}_k$

${\hat{y}}_k=h\left(f\right({\hat{x}}_{s,k-1}\left)\right)$

In formula,Representing the level and smooth state estimation vector in k-1 moment, L is Kalman smoothing gain,It it is the estimation work output in k moment.

(7), on the basis of previous step, by penalty function method, a penalty term constrained optimization problem is converted into a unconstrained optimization problem by being added by former objective function.

The calculation formula of conversion is as follows:

F (x)=f (x)+h (k) H (x), x �� Rⁿ

In formula, f (x) is former objective function, and h (k) is the penalty value dynamically updated, generallyOrK is current iteration number of times. H (x) is the punishment factor, and calculation formula is as follows:

$H\left(x\right)=\underset{i=1}{\overset{m}{\&Sigma;}}\mathrm{\&theta;}\left(q_i\right(x\left)\right)q_i{\left(x\right)}^{\mathrm{\&gamma;}\left(q_i\right(x\left)\right)}$

In formula, �� (q_i(x)) it is multistage partition function, q_iX () is the function relevant with violating constraint condition, q_i(x)=max{0, g_i(x) }, i=1 ..., s, wherein g_i(x)=Dx-d. �� (q_i(x)) represent penalty function effect. The rule that function value is followed is:

$\mathrm{\&gamma;}\left(q_i\right(x\left)\right)=\left\{\begin{array}{cc}1& q_i\left(x\right)\&le;1\\ 2& q_i\left(x\right)>1\end{array}\right.$

$\mathrm{\&theta;}\left(q_i\right(x\left)\right)=\left\{\begin{array}{cc}10& q_i\left(x\right)\&le;0.001\\ 20& 0.001<q_i\left(x\right)\&le;0.1\\ 100& 0.1<q_i\left(x\right)\&le;1\\ 300& q_i\left(x\right)>1\end{array}\right.$

(8), modified particle swarm optiziation then can be utilized on the basis of previous step to carry out successive ignition optimizing.

Wherein, the rule that the population speed of improvement and location updating are followed is as follows:

$V_i^{iter+1}=w\left(L\right)\&times;V_i^{iter}+c_1{r}_{1,i}({\mathrm{Pbest}}_i^{iter}-X_i^{iter})+c_2{r}_{2,i}({\mathrm{Gbest}}^{iter}-X_i^{iter}))$

$X_i^{iter+1}=X_i^{iter}+\mathrm{\&alpha;}\&CenterDot;V_i^{iter+1}$

If search space is D dimension, population comprises N number of particle. In formulaWithRepresent i-th particle respectively at the speed in moment iter and position vector. �� is the contraction factor for controlling and limit speed, c₁And c₂It is respectively cognitive and society's coefficient (generally equal can be set to 2), r_1,iAnd r_2,iIt is two independent randomized numbers of value in [0,1] scope.Represent i-th particle by the iter moment to historical position optimum value. Gbest^iterRepresent in all particles by the historical position optimum value to the iter moment. W (L) is dynamic inertia weight, and for regulating search spatial dimension, optimum to obtaining the overall situation, its calculation formula is:

W (L)=w_start��(w_start-w_end)��(T_max-L)/T_max

In formula, w_startIt is initial inertia weight, w_endBe iteration to inertia weight during maximum times, L is the iteration number of times of the current optimizing of population, T_maxFor the iteration number of times upper limit of population optimizing, in general inertia weights get w_start=0.9, w_endWhen=0.4, population optimizing effect is better.

(9) if iteration number of times L > M, then population optimizing iteration terminates, now the optimizing result Gbest to k moment state estimation^LAs the state estimation in k moment, then subsequent time state is estimated.

(10) if k=k+1��S, then iteration continues, if k=k+1 is > S, then iteration terminates, and exports identification result.

Generally dynamic oscillation signal can represent the sum of the positive string signal for multiple exponential attenuation, it is possible to is described as following form:

$y\left(t\right)=\underset{i=1}{\overset{N}{\&Sigma;}}{A}_i{e}^{-{\mathrm{\&delta;}}_{i}t}cos(w_{i}t+{\mathrm{\&phi;}}_i)+n\left(t\right)$

In formula, A_i,��_i,w_i,��_iBeing the unknown parameter of real number, n (t) is a zero-mean white Gaussian noise. Wherein, ��_iThe damping factor being called Dynamic Signal, w_iIt is the frequency of Dynamic Signal, wherein w_i, ��_iFor solve for parameter. Can obtain the state variables component of Dynamic Signal comprises the separate manufacturing firms model of solve for parameter through reasoning. Considering the Dynamic Signal being made up of the positive string signal summation of N number of exponential attenuation, its 4N state variables form can be expressed as follows:

$x_{4i-3,k}=A_i{e}^{-{\mathrm{\&delta;}}_i\frac{k}{f_s}}cos\left(w_i\frac{k}{f_s}\right)$

$x_{4i-2,k}=A_i{e}^{-{\mathrm{\&delta;}}_i\frac{k}{f_s}}\sin \left(w_i\frac{k}{f_s}\right)$

x_4i-1,k=w_i

x_4i,k=��_i

In formula, i represents these variablees and parameter is i-th attenuated sinusoidal signal belonging to Dynamic Signal. K represents the moment, f_sRepresent sample frequency. The state component in k+1 moment can be obtained according to reasoning:

$x_{4i-3,k+1}=e^{-\frac{x_{4i,k}}{f_s}}\&lsqb;x_{4i-3,k}cos\left(\frac{x_{4i-1,k}}{f_s}\right)-x_{4i-2,k}sin\left(\frac{x_{4i-1,k}}{f_s}\right)\&rsqb;+w_{4i-3,k}$

$x_{4i-2,k+1}=e^{-\frac{x_{4i,k}}{f_s}}\&lsqb;x_{4i-3,k}sin\left(\frac{x_{4i-1,k}}{f_s}\right)+x_{4i-2,k}cos\left(\frac{x_{4i-1,k}}{f_s}\right)\&rsqb;+w_{4i-2,k}$

x_4i-1,k+1=x_4i-1,k+w_4i-1,k

x_4i,k+1=x_4i,k+w_4i,k

Output equation is:

$y_k=\underset{i=1}{\overset{N}{\&Sigma;}}{k}_{2i-1}{x}_{4i-3,k}+k_{2i}{x}_{4i-2,k}+n_k$

In formula, k_2i-1=cos (��_i), k_2i=-sin (��_i), n_kFor average is the white Gaussian noise of zero. So, the state-space model of dynamic oscillation signal generally can represent and is:

$\left\{\begin{array}{c}{x}_{k+1}=f\left(x_k\right)+w_k\\ y_k=h\left(x_k\right)+v_k\end{array}\right.$

In formula, f () and h () represents the nonlinear function that can carry out linearizing according to Taylor series expansion, w_kAnd v_kBe average it is the Gaussian sequence of zero, meets covariance matrix Q respectively_kAnd R_k. Specifically, in dynamic oscillation signal:

$\begin{array}{cccc}f\left(x_k\right)=\left[\begin{array}{c}{M}_1\\ M_2\\ .\\ .\\ .\\ M_i\\ .\\ .\\ .\\ M_N\end{array}\right]& M_i=\left[\begin{array}{c}{x}_{4i-3,k}\\ x_{4i-2,k}\\ x_{4i-1,k}\\ x_{4i,k}\end{array}\right],& Q_k=E\&lsqb;w_k{w}_k^T\&rsqb;,& R_k=E\&lsqb;v_k{v}_k^T\&rsqb;\end{array}$

And function h (x_k) form can be expressed as:

H=(k₁k₂00��,k_2i-1k_2i00��,k_2N-1k_2N00)

h(x_k)=Hx_k

So far, the state-space expression comprising dynamic oscillation signal solve for parameter in state variables component is set up. On this basis, then can use the method that the present invention introduces, namely traditional spreading kalman filtering and improvement constrained procedure, penalty function method and modified particle swarm optiziation are combined, consider the physical constraint condition suffered by solve for parameter simultaneously, carry out dynamic oscillation Signal parameter estimation, obtain the identification result that tool is of practical significance.

Introduce one embodiment of the present of invention below:

Consider that dynamic oscillation signal is:

$\left\{\begin{array}{cc}y\left(k\right)=e^{-0.02k}sin\left(0.5k\right)+n_k,& 0\&le;k\&le;200\\ y\left(k\right)=e^{-0.01k}sin\left(0.6k\right)+n_k,& 200\&le;k\&le;400\end{array}\right.$

In formula, k is the signal sampling moment, n_kIt it is white Gaussian noise. As shown in Figure 2, this Dynamic Signal is made up of the positive string signal of an exponential attenuation, and, this Dynamic Signal was when sampling instant 200 seconds, and signal frequency and damping factor there occurs Spline smoothing. As shown in Figure 3-4, wherein within the moment of 0��k��200, the frequency of dynamic oscillation signal is w=0.5, and damping factor is ��=0.02. Within the moment of 200��k��400, the frequency of dynamic oscillation signal is w=0.6, and damping factor is ��=0.01. When using method proposed by the invention to carry out dynamic signal parameter identification, the relevant initial parameter value that spreading kalman filtering adopts is:

$\begin{array}{ccc}{\hat{P}}_0=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right],& Q_k=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& {10}^{-5}& 0\\ 0& 0& 0& {10}^{-7}\end{array}\right],& R_k={10}^{-5}\end{array}$

${\hat{x}}_0={\left[\begin{array}{cccc}0& 0& 0& 0\end{array}\right]}^T$

When this example is carried out identification, the improve PSO algorithm correlation parameter got is: ��=1, c₁=2, c₂=2, maximum optimizing iteration number of times M=200, it is N=100 that population comprises particle number, initial inertia weight w_start=1.2, termination inertia weight is w_end=0.4. Here, the physical constraint condition suffered by dynamic oscillation signal damping factor �� and frequency w is w >=0, and �� >=0.

Fig. 1 is embodiment algorithm flow figure used, Fig. 2 is the dynamic oscillation signal of embodiment, Fig. 3 uses method proposed by the invention to Dynamic Signal frequency w identification result, and Fig. 4 uses method proposed by the invention to the identification result of dynamic oscillation signal damping factor ��.

Claims (6)

1. the dynamic oscillation signal parameter identification method based on improvement constraint EKF algorithm, it is characterised in that, comprise following steps:

(1) state-space expression comprising dynamic oscillation signal parameter to be identified in state component, is obtained;

(2), initialize, comprising: setting parameter identification initial valueInitial parameter identification error covarianceAnd covariance matrix Q and R that process noise and measurement noise are met, total algorithm iteration number of times maximum value S and population optimizing maximum iteration time M;

(3), by the state estimation in the k-1 moment obtained and state estimation error covariance, utilize the prediction of spreading kalman filtering to walk, obtain status predication value and the status predication error covariance in k moment;

(4), on previous step basis, utilize the filtering of spreading kalman filtering to walk, obtain the state estimation in k moment;

(5), judge whether the parameter identification result in k moment meets corresponding physical constraint condition. If meeting, then directly use EKF iteration identification again;

(6) if not meeting, then using the bounding algorithm of improvement this moment parameter identification result to be retrained, the state estimation problem of constraint equivalence being converted into the optimization problem of constraint;

(7), on the basis of previous step, by penalty function method, a penalty term constrained optimization problem is converted into a unconstrained optimization problem by being added by former objective function;

(8), modified particle swarm optiziation then can be utilized on the basis of previous step to carry out successive ignition optimizing;

(9) if iteration number of times L > M, then population optimizing iteration terminates, now the optimizing result Gbest to k moment state estimation^LAs the state estimation in k moment, then subsequent time state is estimated;

(10) if k=k+1��S, then iteration continues, if k=k+1 is > S, then iteration terminates, and exports identification result.

2. as claimed in claim 1 based on the dynamic oscillation signal parameter identification method improving constraint EKF algorithm, it is characterised in that, the status predication value in k moment and status predication error covariance, calculation formula is:

${\tilde{x}}_k=f({\hat{x}}_{k-1},u_{k-1})$

${\tilde{P}}_k=F_{k-1}{\hat{P}}_{k-1}{F}_{k-1}^T+Q_{k-1}$

In formula,Represent the status predication value in k moment, the nonlinear function in the corresponding particular problem state equation of f (),Represent the state estimation vector in k-1 moment, u_k-1Represent the control inputs in k-1 moment.Represent the status predication error covariance in k moment,Represent that nonlinear function f () existsPlace refined gram than matrix,Representing the state estimation error covariance in k-1 moment, subscript T represents transposition, Q_k-1It it is the covariance matrix met in the system noise k-1 moment.

3. as claimed in claim 1 based on the dynamic oscillation signal parameter identification method improving constraint EKF algorithm, it is characterised in that, the calculation procedure of the state estimation in k moment is:

$K_k={\tilde{P}}_k{H}_k^T{(H_k{\tilde{P}}_k{H}_k^T+R_k)}^{-1}$

${\hat{P}}_k=(I-K_k{H}_k){\tilde{P}}_k$

${\hat{x}}_k={\tilde{x}}_k+K_k\&lsqb;y_k-h\left({\tilde{x}}_k\right)\&rsqb;$

In formula, K_kRepresent the Kalman filtering gain in k moment,Representing the status predication error covariance in k moment, subscript T represents transposition,Represent that nonlinear function h () existsRefined gram of place than matrix, the wherein nonlinear function in the corresponding particular problem output equation of h (). R_kBeing the covariance matrix met in the measurement noise k moment, I is the unit matrix identical with state vector dimension degree,Represent the state estimation vector in k moment, y_kIt it is the work output of k moment output equation.

4. as claimed in claim 1 based on the dynamic oscillation signal parameter identification method improving constraint EKF algorithm, it is characterised in that, constrained optimization objective function is:

${\widehat{\stackrel{\&OverBar;}{x}}}_k=\underset{x_k}{\arg min}{(y_k-{\hat{y}}_k)}^T{W}_1(y_k-{\hat{y}}_k)+{(x_k-{\hat{x}}_k)}^T{W}_2(x_k-{\hat{x}}_k)$

In formula,Represent k moment optimizing state estimation to be asked, y_kIt is k moment work output, W₁And W₂It is the weight matrix of positive definite, it is chosen for unit matrix, x in formula_kMeet constraint condition as shown in the formula:

Dx��d

The capable full rank matrix of s �� n constant that D is known, s is the number by constraint condition parameter, and n is the dimension of state vector, it is clear that s��n, d are known constrainted constants conditions;

In constrained optimization objective function, to the estimation work output in k momentWhen solving, in order to make full use of known measurement data information, adopting Kalman smoothing method, calculation formula is as follows:

${\hat{x}}_{s,k-1}={\hat{x}}_{k-1}+L\left({\hat{x}}_k-{\tilde{x}}_k\right)$

$L={\hat{P}}_{k-1}{F}_{k-1}{\tilde{P}}_k$

${\hat{y}}_k=h\left(f\right({\hat{x}}_{s,k-1}\left)\right)$

In formula,Representing the level and smooth state estimation vector in k-1 moment, L is Kalman smoothing gain,It it is the estimation work output in k moment.

5. as claimed in claim 1 based on the dynamic oscillation signal parameter identification method improving constraint EKF algorithm, it is characterised in that, constrained optimization problem is converted into a unconstrained optimization problem;

The calculation formula of conversion is as follows:

F (x)=f (x)+h (k) H (x), x �� Rⁿ

In formula, f (x) is former objective function, and h (k) is the penalty value dynamically updated, generallyOrK is current iteration number of times. H (x) is the punishment factor, and calculation formula is as follows:

$H\left(x\right)=\underset{i=1}{\overset{m}{\&Sigma;}}\mathrm{\&theta;}\left(q_i\right(x\left)\right)q_i{\left(x\right)}^{\mathrm{\&gamma;}\left(q_i\right(x\left)\right)}$

In formula, �� (q_i(x)) it is multistage partition function, q_iX () is the function relevant with violating constraint condition, q_i(x)=max{0, g_i(x) }, i=1 ..., s, wherein g_i(x)=Dx-d; �� (q_i(x)) represent penalty function effect. The rule that function value is followed is:

$\mathrm{\&gamma;}\left(q_i\right(x\left)\right)=\left\{\begin{array}{cc}1& q_i\left(x\right)\&le;1\\ 2& q_i\left(x\right)>1\end{array}\right.$

$\mathrm{\&theta;}\left(q_i\left(x\right)\right)=\left\{\begin{array}{cc}10& q_i\left(x\right)\&le;0.001\\ 20& 0.001<q_i\left(x\right)\&le;0.1\\ 100& 0.1<q_i\left(x\right)\&le;1\\ 300& q_i\left(x\right)>1\end{array}\right..$

6. as claimed in claim 1 based on the dynamic oscillation signal parameter identification method improving constraint EKF algorithm, it is characterised in that, the rule that the population speed of improvement and location updating are followed is as follows:

$V_i^{iter+1}=w\left(L\right)\&times;V_i^{iter}+c_1{r}_{1,i}({\mathrm{Pbest}}_i^{iter}-X_i^{iter})+c_2{r}_{2,i}\left({\mathrm{Gbest}}^{iter}-X_i^{iter}\right))$

$X_i^{iter+1}=X_i^{iter}+\mathrm{\&alpha;}\&CenterDot;V_i^{iter+1}$

If search space is D dimension, population comprises N number of particle. In formulaWithRepresent i-th particle respectively at the speed in moment iter and position vector. �� is the contraction factor for controlling and limit speed, c₁And c₂It is respectively cognitive and society's coefficient, r_1,iAnd r_2,iIt is two independent randomized numbers of value in [0,1] scope;Represent i-th particle by the iter moment to historical position optimum value; Gbest^iterRepresent in all particles by the historical position optimum value to the iter moment; W (L) is dynamic inertia weight, and for regulating search spatial dimension, optimum to obtaining the overall situation, its calculation formula is:

W (L)=w_start��(w_start-w_end)��(T_max-L)/T_max

In formula, w_startIt is initial inertia weight, w_endBe iteration to inertia weight during maximum times, L is the iteration number of times of the current optimizing of population, T_maxFor the iteration number of times upper limit of population optimizing.