CN107402903B - Nonlinear system state deviation evolution method based on differential algebra and Gaussian sum - Google Patents

Nonlinear system state deviation evolution method based on differential algebra and Gaussian sum Download PDF

Info

Publication number
CN107402903B
CN107402903B CN201710551107.9A CN201710551107A CN107402903B CN 107402903 B CN107402903 B CN 107402903B CN 201710551107 A CN201710551107 A CN 201710551107A CN 107402903 B CN107402903 B CN 107402903B
Authority
CN
China
Prior art keywords
gaussian distribution
sub
deviation
matrix
gaussian
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201710551107.9A
Other languages
Chinese (zh)
Other versions
CN107402903A (en
Inventor
罗亚中
孙振江
张进
杨震
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
National University of Defense Technology
Original Assignee
National University of Defense Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by National University of Defense Technology filed Critical National University of Defense Technology
Priority to CN201710551107.9A priority Critical patent/CN107402903B/en
Publication of CN107402903A publication Critical patent/CN107402903A/en
Application granted granted Critical
Publication of CN107402903B publication Critical patent/CN107402903B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/13—Differential equations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Operations Research (AREA)
  • Computing Systems (AREA)
  • Other Investigation Or Analysis Of Materials By Electrical Means (AREA)
  • Complex Calculations (AREA)

Abstract

The invention discloses a nonlinear system state deviation evolution method based on differential algebra and Gaussian sum, comprising the following steps: forecasting the terminal state of the nonlinear system by adopting a differential algebra method, and expressing the terminal state as a high-order Taylor expansion polynomial about the initial state deviation; determining a sub-Gaussian distribution covariance matrix, and fitting a Gaussian sum model aiming at each sub-Gaussian distribution target practice; calculating a sub-Gaussian distribution high-order central moment; and determining the mean value and covariance matrix of each sub-Gaussian distribution at the terminal moment, and giving a Gaussian sum form terminal state deviation distribution probability density function. The method can be automatically expanded to any specified order of deviation evolution precision, does not need to manually deduce the high-order partial derivative of a kinetic equation, is suitable for the nonlinear system deviation evolution analysis problem with strong nonlinearity and long-time prediction, can still keep the remarkable efficiency advantage compared with the Monte Carlo simulation method while accurately describing the deviation distribution and the non-Gaussian property thereof, and has the advantages of convenience in use and high calculation precision.

Description

Nonlinear system state deviation evolution method based on differential algebra and Gaussian sum
Technical Field
The invention relates to the field of spacecraft orbit dynamics, in particular to a nonlinear system state deviation evolution method based on differential algebra and Gaussian sum, which is used for predicting the state of a nonlinear system and the deviation thereof.
Background
The spacecraft trajectory deviation evolution technology is an important basic technology for developing space missions such as space collision early warning, spacecraft orbit determination, navigation filtering, spacecraft trajectory safety analysis and design and the like. The forecasting precision of the track deviation directly influences the credibility of early warning, orbit determination and navigation results, and even determines the success or failure of the space mission. Because the essence of the spacecraft orbit dynamics model is a nonlinear system in the form of a system of ordinary differential equations or a state transition equation, the spacecraft orbit deviation evolution problem is the nonlinear system state deviation forecasting problem.
Classical nonlinear system bias evolution methods are linearized covariance analysis methods and MonteCarlo simulation methods. The linear covariance analysis method is characterized in that a nonlinear dynamical system is linearized, and then the linear covariance analysis method is adopted for deviation forecasting calculation. According to the MonteCarlo simulation method, high-precision state deviation evolution results can be obtained by counting multiple target practice test results, but a large amount of sampling simulation is required by using the MonteCarlo simulation method, the calculation amount is large, and the calculation efficiency is difficult to meet the task requirements under many conditions.
The differential algebra method is a method for calculating the high-order partial derivatives of the nonlinear mapping with respect to the independent variables in a numerical environment, and taylor expansion of any order can be performed on the nonlinear system by using the high-order partial derivatives, so that the state variables of the system are predicted by using high-order polynomials. The method was originally proposed by MartinBerz, professor of particle physics, university of Michigan, USA, and the basic principle can be found in the references: berz M.Differencential algebraic description of beam dynamics to very high orders [ J ]. part.Accel.,1988,24(SSC-152): 109-. At present, no research or application result for applying the method to nonlinear dynamical system deviation prediction exists in China.
The Gaussian sum model is a method for describing the state deviation distribution of the system in the form of weighted sum of a plurality of sub-Gaussian distributions, and can be used for accurately describing the state deviation probability density function of non-Gaussian distribution. The invention is oriented to a general nonlinear dynamical system, a system dynamical model can be in a state transfer equation form or a normal differential equation form, a high-order polynomial of the terminal state about the initial state deviation is solved by using a differential algebra method, further Gaussian sum models are adopted to fit the initial state deviation distribution, and each sub-Gaussian distribution of the terminal time is predicted based on the high-order polynomial, so that the non-Gaussian distribution characteristic of the terminal state is accurately described.
Disclosure of Invention
The technical problems to be solved by the invention are as follows: aiming at the problems in the prior art, the nonlinear system state deviation evolution method based on differential algebra and Gaussian sum can be automatically expanded to any specified order of deviation evolution precision, does not need to manually derive the high-order partial derivative of a kinetic equation, is suitable for nonlinear system deviation evolution analysis problems with strong nonlinearity and long-time prediction, can accurately describe deviation distribution and non-Gaussian property thereof, can keep the remarkable efficiency advantage relative to a MonteCarlo simulation method, is convenient to use, and has high calculation precision.
In order to solve the technical problems, the invention adopts the technical scheme that:
a nonlinear system state deviation evolution method based on differential algebra and Gaussian sum is used for spacecraft trajectory deviation evolution, and comprises the following steps:
1) inputting a kinetic equation, a deviation evolution order N and a relative nominal value of an initial state of a nonlinear system
Figure GDA0002884026650000026
And initial bias covariance matrix P0；
2) Forecasting the terminal state of the nonlinear system by adopting a differential algebra method according to the kinetic equation of the nonlinear system and expressing the terminal state as the deviation delta x of the initial state0A high order taylor expansion polynomial of (1);
3) for initial bias covariance matrix P0Determining a sub-Gaussian distribution covariance matrix PiFor sub-Gaussian distribution covariance matrix PiPerforming target fitting on the initial deviation distribution by the multiple sub-Gaussian distributions to generate a Gaussian sum model;
4) calculating a sub-Gaussian distribution high-order central moment of a given order aiming at the generated Gaussian sum model;
5) and determining the mean value and covariance matrix of each sub-Gaussian distribution at the terminal moment based on the high-order polynomial of the terminal state relative to the initial state deviation and each order central moment of the sub-Gaussian distribution, and further providing a terminal state deviation distribution probability density function in a Gaussian sum form.
Preferably, the kinetic equation of the nonlinear system in the step 1) is in the form of a nonlinear mapping xf＝f(x0) Or form of ordinary differential equation
Figure GDA0002884026650000021
And the function f (-) in the form of the nonlinear mapping or the function g (-) in the form of the ordinary differential equation is an arbitrary n-dimensional real number domain RnTo RnNon-linear mapping of (2).
Preferably, the detailed steps of step 2) include:
2.1) initial state x0The initialization is represented by formula (1) and includes the relative nominal value of the initial state
Figure GDA0002884026650000022
The form of deviation of (1);
Figure GDA0002884026650000023
in the formula (1), x0In the initial state, the state of the device is as follows,
Figure GDA0002884026650000024
relative to nominal value for initial state, δ x0Is an initial state deviation;
2.2) initial state x0Substituting into the dynamic equation of the nonlinear system, and applying a differential algebra method to obtain the deviation deltax about the initial state0The terminal state described by the high-order Taylor expansion polynomial is shown as the formula (2);
Figure GDA0002884026650000025
in the formula (2), xfFor the terminal state, T is a high-order Taylor expansion polynomial, the superscript N of T is a deviation evolution order, and the subscript x of TfThe contents of the polynomial description for the higher order Taylor expansion are terminal states, δ x0Is the initial state deviation.
Preferably, the detailed steps of step 3) include:
3.1) covariance matrix P for initial bias0Inverse of (3) performing Cholesky decomposition to obtain square root information matrix R0And satisfy P0＝R0 -1R0 -T；
3.2) lower bound R of the square root information matrix given a Gaussian distribution of the childrenminAnd obtaining a sub-Gaussian distribution covariance matrix Pi；
3.3) target multiple sub-Gaussian distributions to generate Gaussian sums to fit the initial Gaussian distribution.
Preferably, the sub-Gaussian distribution covariance matrix P is solved in the step 3.2)iThe detailed steps comprise:
3.2.1) computing the matrix using the functional expression shown in equation (3)
Figure GDA0002884026650000031
Singular value decomposition of (1), wherein R0Is a square root information matrix, RminThe lower bound of the square root information matrix of sub-Gaussian distribution;
Figure GDA0002884026650000032
in the formula (3), UiAnd ViAs orthonormal, diagonal, Si＝diag(σi1,...,σin) Wherein the diagonal element σi1,...,σinIs a positive singular value, R0Is a square root information matrix, RminThe lower bound of the square root information matrix of sub-Gaussian distribution;
3.2.2) determining the diagonal matrix S for any k ═ 1iOf (a) the kth diagonal element σikIf 1 or more is true, the square root information matrix R of the sub-Gaussian distribution is obtainediIs a square root information matrix R0Jump to execute step 3.2.5); otherwise, skipping to execute the next step;
3.2.3) establishing an n-dimensional diagonal matrix shown in the vertical type (4), deleting all 0 rows of the n-dimensional diagonal matrix to obtain a singular value change matrix delta Si；
Figure GDA0002884026650000033
In the formula (4), δ Si_fullFor an n-dimensional diagonal matrix, σi1Is a diagonal matrix SiThe 1 st diagonal element of (1), σinIs a diagonal matrix SiThe nth diagonal element of (1);
3.2.4) determining the information variation matrix δ R according to equation (5)iAnd solving a square root information matrix R of sub-Gaussian distribution according to the formula (6)i；
δRi＝δSiVi TRmin (5)
In the formula (5), δ RiFor information change matrix, δ SiAs a matrix of singular value variations, ViIs an orthonormal matrix, RminThe lower bound of the square root information matrix of sub-Gaussian distribution;
Figure GDA0002884026650000041
in the formula (6), QiIs an orthonormal matrix, RiIs a square root information matrix of sub-Gaussian distribution, δ RiFor information change matrix, R0Is a square root information matrix.
3.2.5) calculating a reduction of the sub-Gaussian covariance compared to the initial Gaussian covariance according to equation (7);
δPi＝P0-Pi＝R0 -1R0 -T-Ri -1Ri -T＝δYi(δYi)T (7)
in the formula (7), δ PiIs the amount of reduction of the sub-Gaussian distribution covariance compared to the initial Gaussian distribution covariance, δ YiIs the square root of the covariance reduction, deltaPiSemi-positive and delta YiIs the square root of its matrix, δ YiThe functional expression of (a) is represented by the formula (8);
Figure GDA0002884026650000042
in the formula (8), δ YiIs the square root of the covariance reduction, R0Is a square root information matrix, δ RiFor information change matrix, matrix RdiObtained by QR decomposition as shown in a formula (9);
Figure GDA0002884026650000043
in the formula (9), QdiIs an orthonormal matrix, RdiIs an upper triangular square matrix, R0Is a square root information matrix, δ RiIs an information change matrix, and I is an n-dimensional unit matrix;
3.2.6) covariance matrix P based on initial bias0The amount of reduction of the sub-Gaussian distribution covariance compared to the initial Gaussian distribution covariance, δ PiDetermining a sub-Gaussian distribution covariance matrix Pi。
Preferably, the detailed steps of step 3.3) include:
3.3.1) note nδYMatrix δ Y being the square root of covariance reductioniColumn number of (d), singular value change matrix δ SiSetting the number of target shooting times M as the number of sub-Gaussian distributions, and initializing i;
3.3.2) recording the current target shooting frequency as the ith time;
3.3.3) generating nδYVector of dimension ηiThe components thereof obey standard Gaussian distribution and are independent of each other;
3.3.4) according to μi＝μ0+δYiηiCalculating the mean value mu of the ith sub-Gaussian distributioniIn which μ0Is the mean of the initial states, δ YiIs the square root, η, of the covariance reductioniIs nδYA dimension column vector;
3.3.5) determining the square root information matrix R corresponding to the ith sub-Gaussian distributioniAnd a covariance matrix Pi；
3.3.6) determining the ith sub-Gaussian distributionWeight of w i1/M, wherein M is the number of times of target shooting;
3.3.7) adding 1 to the variable i, judging whether i is smaller than N, and if so, skipping to execute the step 3.3.2); otherwise, skipping to execute the next step;
3.3.8) generating an initial deviation probability density function of the form of a sum of gaussians as shown in equation (10);
Figure GDA0002884026650000051
in the formula (10), p (t, x)0) Initial deviation probability density function in the form of a sum of gaussians, wiThe weight corresponding to the ith sub-Gaussian distribution, M is the number of times of target shooting, pg(x0；μi,Pi) Is the probability density function of the ith sub-Gaussian distribution and the probability density function p of the ith sub-Gaussian distributiong(x0；μi,Pi) The functional expression of (a) is represented by the formula (11);
Figure GDA0002884026650000052
in the formula (11), n is the system dimension, PiCovariance matrix, x, corresponding to ith sub-Gaussian distribution0Is in an initial state, muiIs the mean of the ith sub-gaussian distribution.
Preferably, the detailed steps of step 4) include:
4.1) determining that any odd-order central moment of sub-Gaussian distribution is 0 according to the characteristics of Gaussian distribution;
4.2) for the ith sub-Gaussian distribution, let
Figure GDA0002884026650000053
Wherein a is 1,2aIs the a-th element of the state variable x,
Figure GDA0002884026650000054
is the mean value of the states muiThe a-th element of (1); according to formula (12)Calculating a second-order central moment of the operator Gaussian distribution deviation;
E{δxaδxb}＝Pi ab,a,b∈[1,2,...,n] (12)
in the formula (12), E {. cndot } represents an expected value of a physical quantity in parentheses, E { δ x { [aδxbIs the second central moment of the sub-Gaussian distribution deviation, δ xaIs a state variable x and a state mean value muiDeviation of a-th element, δ xbIs a state variable x and a state mean value muiDeviation of the b-th element, Pi abRepresenting the covariance matrix P corresponding to the ith sub-Gaussian distributioniRow a, column b elements of (a);
4.3) on the basis of the second-order central moment of the known sub-Gaussian distribution deviation, making the m initial value be 2 and sequentially increasing by taking 2 as increment, and obtaining the m-th-order central moment and the covariance matrix PiDeducing and obtaining the m + 2-order central moment of sub-Gaussian distribution according to a formula (13) until the first 2N-order central moment of the sub-Gaussian distribution is calculated, wherein N is a deviation evolution order;
Figure GDA0002884026650000055
in the formula (13), E {. cndot } represents an expected value of a physical quantity in parentheses, E { δ x { }aδxbδxc...δxdA (a, b, c., d) th element which is the m +2 order central moment of the sub-gaussian distribution bias, a, b, c., d being the index of m +2 elements to the state variable x, (a, b, c., d) being the index of the element to the m +2 order central moment; pi abCovariance matrix P representing ith sub-Gaussian distributioniRow a, column b, elements E { δ xc...δxdThe (c, d) th element, which is the m-th order central moment of sub-gaussian distribution bias; pi acCovariance matrix P representing ith sub-Gaussian distributioniRow a, column c, element E { δ x [)b...δxdThe (b, d) th element, which is the m-th order central moment of sub-gaussian distribution bias; and so on until Pi adCovariance representing ith sub-Gaussian distributionDifference matrix PiRow a and column d elements of (1), E { δ xbδxc... } is the (b, c.,) th element of the m-th order central moment of the sub-gaussian distribution bias.
Preferably, the detailed steps of step 5) include:
5.1) traversing and selecting one sub-Gaussian distribution as the current ith sub-Gaussian distribution;
5.2) for the ith sub-Gaussian distribution, the initial state is compared with the nominal value according to the formula (14)
Figure GDA0002884026650000064
Initial state deviation δ x of0Transformed into the mean μ of the ith sub-Gaussian distributioniAnd corresponding deviation δ xi；
δx0＝μi+δxi (14)
In the formula (14), δ x0Is the deviation of the initial state, muiIs the mean value of the ith sub-Gaussian distribution, δ xiIs the deviation of the ith sub-Gaussian distribution;
5.3) mean value μ to be transformed into ith sub-Gaussian distributioniAnd corresponding deviation δ xiInitial state deviation δ x of0Substituting the terminal state with respect to the initial state deviation deltax0Obtaining a terminal state shown in a formula (15) in the high-order Taylor expansion polynomial;
Figure GDA0002884026650000061
in the formula (15), xfFor the terminal state, T is a high-order Taylor expansion polynomial, the superscript N of T is a deviation evolution order, and the subscript x of TfAt the end state, δ xiIs the deviation of the ith sub-Gaussian distribution;
5.4) calculating the mean value mu of the ith sub-Gaussian distribution at the terminal moment according to the formula (16)i_fAs the expected value of the terminal state shown in equation (15);
Figure GDA0002884026650000062
in the formula (16), mui_fIs the mean value, x, of the ith sub-Gaussian distribution at the terminal timefFor the terminal state, T is a high-order Taylor expansion polynomial, the superscript N of T is a deviation evolution order, and the subscript x of TfAt the end state, δ xiIs the deviation of the ith sub-Gaussian distribution;
5.5) covariance matrix P of ith sub-Gaussian distribution of terminal time according to equation (17)i_fExpressed as deviation deltax with respect to the initial sub-gaussian distributioniThe 2N-order Taylor expansion and the expected value corresponding to the ith sub-Gaussian distribution are solved;
Figure GDA0002884026650000063
in the formula (17), Pi_fCovariance matrix, x, representing the ith sub-Gaussian distribution at the terminal timefIn the end state, μi_fThe mean value of the ith sub-Gaussian distribution at the terminal time, T is a high-order Taylor expansion polynomial, N in the superscript of T is a deviation evolution order, and deltaxiIs the deviation of the ith sub-Gaussian distribution;
5.6) judging whether all sub-Gaussian distributions are traversed or not, and if not, skipping to execute the step 5.2); otherwise, skipping to execute the next step;
5.7) obtaining the mean value mu of all sub-Gaussian distributions at the terminal momenti_fSum covariance matrix Pi_fOn the basis, an accurate probability density function of the system state distribution at the terminal moment is obtained according to the formula (18) and output;
Figure GDA0002884026650000071
in formula (18), p (t, x)f) An accurate probability density function representing the state distribution of the system at the terminal time, t being the terminal time, xfIs in the terminal state, wiThe weight corresponding to the ith sub-Gaussian distribution, M is the number of times of target shooting, pg(xf；μi_f,Pi_f) Is the probability density function of the ith sub-Gaussian distribution, mui_fIs the mean value of the ith sub-Gaussian distribution at the terminal time, Pi_fAnd the covariance matrix represents the ith sub-Gaussian distribution of the terminal moment.
The nonlinear system state deviation evolution method based on differential algebra and Gaussian sum has the following advantages:
1. the method is based on a differential algebra method, can forecast the state of the nonlinear system with any specified order precision, and adopts a recursion algorithm to calculate the Gaussian distribution central moment of any order, so the nonlinear system state deviation evolution method based on differential algebra and Gaussian sum can be automatically expanded to any specified order deviation evolution precision, the high-order partial derivative of a kinetic equation does not need to be manually deduced, and the method has the advantages of convenient use and high calculation precision, and is suitable for the nonlinear system deviation evolution analysis problem with strong nonlinearity and long-time forecast.
2. The invention adopts Gauss and a model to represent the probability density function of the state deviation distribution of the system, has the accurate description capacity of the nonlinear system deviation non-Gauss, and each sub-Gauss distribution forecast is based on the high-order Taylor expansion polynomial of the same terminal state relative to the initial state deviation, and does not need to forecast the terminal state for each sub-Gauss distribution, so that the nonlinear system state deviation evolution method based on differential algebra and Gauss and can still keep the obvious efficiency advantage relative to the Monte Carlo simulation method while accurately describing the deviation distribution and the non-Gauss thereof.
3. In the field of spacecraft orbit dynamics, the method can be applied to spacecraft absolute or relative trajectory prediction and deviation evolution analysis thereof, orbit determination error analysis and the like, and the obtained deviation evolution information can be further used for space collision early warning, spacecraft trajectory safety analysis and design and the like.
Drawings
FIG. 1 is a schematic diagram of a basic process of an embodiment of the present invention.
FIG. 2 is a diagram of a method for obtaining a sub-Gaussian distribution covariance matrix P according to an embodiment of the present inventioniIs a schematic flow diagram.
Fig. 3 is a diagram illustrating a comparison of distribution of terminal state deviations according to a first embodiment of the present invention.
Fig. 4 is a comparison graph of the terminal probability density function according to the first embodiment of the present invention.
Fig. 5 is a diagram comparing distribution of terminal state deviations according to the second embodiment of the present invention.
Fig. 6 is a comparison graph of the terminal probability density function according to the second embodiment of the present invention.
Detailed Description
The first embodiment is as follows:
as shown in fig. 1, the method for evolving state deviation of a nonlinear system based on differential algebra and gaussian sum in this embodiment includes the following steps:
1) inputting a kinetic equation, a deviation evolution order N and a relative nominal value of an initial state of a nonlinear system
Figure GDA0002884026650000082
And initial bias covariance matrix P0；
2) Forecasting the terminal state of the nonlinear system by adopting a differential algebra method according to the kinetic equation of the nonlinear system and expressing the terminal state as the deviation delta x of the initial state0A high order taylor expansion polynomial of (1);
3) for initial bias covariance matrix P0Determining a sub-Gaussian distribution covariance matrix PiFor sub-Gaussian distribution covariance matrix PiPerforming target fitting on the initial deviation distribution by the multiple sub-Gaussian distributions to generate a Gaussian sum model;
4) calculating a sub-Gaussian distribution high-order central moment of a given order aiming at the generated Gaussian sum model;
5) and determining the mean value and covariance matrix of each sub-Gaussian distribution at the terminal moment based on the high-order polynomial of the terminal state relative to the initial state deviation and each order central moment of the sub-Gaussian distribution, and further providing a terminal state deviation distribution probability density function in a Gaussian sum form.
In the present embodiment, a differential algebra method is adopted in a nonlinear system state deviation evolution method based on differential algebra and gaussian sum to predict a terminal state of a nonlinear system, and the terminal state is expressed as a high-order polynomial about initial state deviation; fitting the initial deviation distribution in the form of weighted sum of multiple sub-Gaussian distributions by adopting a Gaussian sum model; calculating a sub-Gaussian distribution high-order central moment of a given order by using a covariance matrix of sub-Gaussian distribution; and determining the mean value and the covariance matrix of each sub-Gaussian distribution at the terminal moment based on the high-order polynomial of the terminal state deviation relative to the initial state and each order moment of the sub-Gaussian distribution. The invention is oriented to a general nonlinear dynamical system, a dynamical model can be in a state transition equation form or an ordinary differential equation form, and any high-order state prediction and deviation evolution analysis can be automatically carried out. In the field of spacecraft orbit dynamics, the method can be particularly applied to spacecraft absolute or relative trajectory prediction and deviation evolution analysis thereof, orbit determination error analysis and the like, and the obtained deviation evolution information can be further used for space collision early warning, spacecraft trajectory safety analysis and design and the like.
The nonlinear system state deviation evolution method based on differential algebra and gaussian sum in the embodiment is oriented to a general nonlinear dynamical system, the dynamical model can be in a state transition equation form or an ordinary differential equation form, and any high-order state prediction and deviation evolution analysis can be automatically performed. In this embodiment, the kinetic equation of the nonlinear system in step 1) is in the form of nonlinear mapping xf＝f(x0) Or form of ordinary differential equation
Figure GDA0002884026650000081
And the function f (-) in the form of the nonlinear mapping or the function g (-) in the form of the ordinary differential equation is an arbitrary n-dimensional real number domain RnTo RnNon-linear mapping of (2).
In this embodiment, the nonlinear system is a kepler orbit of a near-earth spacecraft, the system dimension n is 6, and the state variable x is the sum of the positions of the spacecraft in an earth inertial coordinate systemVelocity vector, i.e. x ═ rT vT]TWhere r is the position vector and v is the velocity vector. The dynamic equation is a Kepler orbit analytic solution, the specific solving form of the Kepler orbit analytic solution can be found in the existing spacecraft orbit dynamics monograph, and the spacecraft terminal state x can be obtainedfViewed as initial state x0And the track forecast time tfOf non-linear functions, i.e. like xf＝f(tf,x0) The nominal value of the classical orbit number of the spacecraft at the given initial moment is as follows:
E0＝[a0,e0,i0,Ω0,ω0,f0]＝[6778138,0.001,0.001,0.001,0.001,0.001]
wherein a is0Is a semi-major axis, e0As eccentricity, i0For track inclination, omega0Is the ascending crossing right ascension, omega0Is an angular distance of a near arch point, f0Is the true proximal angle, with the length unit being m and the angle unit being rad. According to the initial nominal number of tracks E0The nominal state of the spacecraft at the initial moment can be calculated
Figure GDA0002884026650000091
Given track forecast time tfIs one track period, i.e. tf＝5553.63s。
In addition, in the embodiment, given the evolution order N of the bias 4, the covariance matrix P of the initial bias is given0Comprises the following steps:
Figure GDA0002884026650000092
in this embodiment, the detailed steps of step 2) include:
2.1) initial state x0The initialization is represented by formula (1) and includes the relative nominal value of the initial state
Figure GDA0002884026650000093
The form of deviation of (1);
Figure GDA0002884026650000094
in the formula (1), x0In the initial state, the state of the device is as follows,
Figure GDA0002884026650000095
relative to nominal value for initial state, δ x0Is an initial state deviation;
2.2) initial state x0Substituting into the dynamic equation of the nonlinear system, and applying a differential algebra method to obtain the deviation deltax about the initial state0The terminal state described by the high-order Taylor expansion polynomial is shown as the formula (2);
Figure GDA0002884026650000096
in the formula (2), xfFor the terminal state, T is a high-order Taylor expansion polynomial, the superscript N of T is a deviation evolution order, and the subscript x of TfThe contents of the polynomial description for the higher order Taylor expansion are terminal states, δ x0Is the initial state deviation.
In this embodiment, the detailed steps of step 3) include:
3.1) covariance matrix P for initial bias0Inverse of (3) performing Cholesky decomposition to obtain square root information matrix R0And satisfy P0＝R0 -1R0 -T；
3.2) lower bound R of the square root information matrix given a Gaussian distribution of the childrenminAnd obtaining a sub-Gaussian distribution covariance matrix Pi；
3.3) target multiple sub-Gaussian distributions to generate Gaussian sums to fit the initial Gaussian distribution.
As shown in fig. 2, the sub-gaussian distribution covariance matrix P is solved in step 3.2)iThe detailed steps comprise:
3.2.1) computing the matrix using the functional expression shown in equation (3)
Figure GDA0002884026650000101
Singular value ofDecomposition of wherein R0Is a square root information matrix, RminThe lower bound of the square root information matrix of sub-Gaussian distribution;
Figure GDA0002884026650000102
in the formula (3), UiAnd ViAs orthonormal, diagonal, Si＝diag(σi1,...,σin) Wherein the diagonal element σi1,...,σinIs a positive singular value, R0Is a square root information matrix, RminThe lower bound of the square root information matrix of sub-Gaussian distribution; in this embodiment, the lower bound of the sub-gaussian square root information matrix is set as Rmin＝5R0；
3.2.2) determining the diagonal matrix S for any k ═ 1iOf (a) the kth diagonal element σikIf 1 or more is true, the square root information matrix R of the sub-Gaussian distribution is obtainediIs a square root information matrix R0Jump to execute step 3.2.5); otherwise, skipping to execute the next step;
3.2.3) establishing an n-dimensional diagonal matrix shown in the vertical type (4), deleting all 0 rows of the n-dimensional diagonal matrix to obtain a singular value change matrix delta Si；
Figure GDA0002884026650000103
In the formula (4), δ Si_fullFor an n-dimensional diagonal matrix, σi1Is a diagonal matrix SiThe 1 st diagonal element of (1), σinIs a diagonal matrix SiThe nth diagonal element of (1);
3.2.4) determining the information variation matrix δ R according to equation (5)iAnd solving a square root information matrix R of sub-Gaussian distribution according to the formula (6)i；
δRi＝δSiVi TRmin (5)
In the formula (5), δ RiFor information change matrix, δ SiAs a matrix of singular value variations, ViIs an orthonormal matrix, RminThe lower bound of the square root information matrix of sub-Gaussian distribution;
Figure GDA0002884026650000104
in the formula (6), QiIs an orthonormal matrix, RiIs a square root information matrix of sub-Gaussian distribution, δ RiFor information change matrix, R0Is a square root information matrix.
3.2.5) calculating a reduction of the sub-Gaussian covariance compared to the initial Gaussian covariance according to equation (7);
δPi＝P0-Pi＝R0 -1R0 -T-Ri -1Ri -T＝δYi(δYi)T (7)
in the formula (7), δ PiIs the amount of reduction of the sub-Gaussian distribution covariance compared to the initial Gaussian distribution covariance, δ YiIs the square root of the covariance reduction, deltaPiSemi-positive and delta YiIs the square root of its matrix, δ YiThe functional expression of (a) is represented by the formula (8);
Figure GDA0002884026650000111
in the formula (8), δ YiIs the square root of the covariance reduction, R0Is a square root information matrix, δ RiFor information change matrix, matrix RdiObtained by QR decomposition as shown in a formula (9);
Figure GDA0002884026650000112
in the formula (9), QdiIs an orthonormal matrix, RdiIs an upper triangular square matrix, R0Is the square root moment of informationArray, δ RiIs an information change matrix, and I is an n-dimensional unit matrix;
3.2.6) covariance matrix P based on initial bias0The amount of reduction of the sub-Gaussian distribution covariance compared to the initial Gaussian distribution covariance, δ PiDetermining a sub-Gaussian distribution covariance matrix Pi。
In this embodiment, the detailed steps of step 3.3) include:
3.3.1) note nδYMatrix δ Y being the square root of covariance reductioniColumn number of (d), singular value change matrix δ SiSetting the number of target shooting times M as the number of sub-Gaussian distributions, and initializing i; in this embodiment, the number M of sub-gaussian distributions is 5000;
3.3.2) recording the current target shooting frequency as the ith time;
3.3.3) generating nδYVector of dimension ηiThe components thereof obey standard Gaussian distribution and are independent of each other;
3.3.4) according to μi＝μ0+δYiηiCalculating the mean value mu of the ith sub-Gaussian distributioniIn which μ0Is the mean of the initial states, δ YiIs the square root, η, of the covariance reductioniIs nδYA dimension column vector;
3.3.5) determining the square root information matrix R corresponding to the ith sub-Gaussian distributioniAnd a covariance matrix Pi；
3.3.6) determining the weight corresponding to the ith sub-Gaussian distribution as w i1/M, wherein M is the number of times of target shooting;
3.3.7) adding 1 to the variable i, judging whether i is smaller than N, and if so, skipping to execute the step 3.3.2); otherwise, skipping to execute the next step;
3.3.8) generating an initial deviation probability density function of the form of a sum of gaussians as shown in equation (10);
Figure GDA0002884026650000113
in the formula (10), p (t, x)0) Initial bias profile in the form of a sum of gaussiansRate density function, wiThe weight corresponding to the ith sub-Gaussian distribution, M is the number of times of target shooting, pg(x0；μi,Pi) Is the probability density function of the ith sub-Gaussian distribution and the probability density function p of the ith sub-Gaussian distributiong(x0；μi,Pi) The functional expression of (a) is represented by the formula (11);
Figure GDA0002884026650000121
in the formula (11), n is the system dimension, PiCovariance matrix, x, corresponding to ith sub-Gaussian distribution0Is in an initial state, muiIs the mean of the ith sub-gaussian distribution.
In this embodiment, the detailed steps of step 4) include:
4.1) determining that any odd-order central moment of sub-Gaussian distribution is 0 according to the characteristics of Gaussian distribution;
for example, the first order central moment of a sub-gaussian distribution is:
E{δxa}＝0
wherein E {. cndot } represents an expected value of a physical quantity in parentheses, E { δ x { [aIs the first central moment of the sub-Gaussian distribution deviation, xaIs the a-th element of the state variable x, δ xaIs a state variable x and a state mean value muiDeviation of the a-th element.
For example, the third central moment of a sub-gaussian distribution is:
E{δxaδxbδxc}＝0,a,b,c∈[1,2,...,n]
wherein E {. cndot } represents an expected value of a physical quantity in parentheses, E { δ x { [aδxbδxcIs the first central moment of the sub-Gaussian distribution deviation, xaIs the a-th element of the state variable x, δ xaIs a state variable x and a state mean value muiDeviation of a-th element, xbIs the b-th element of the state variable x, δ xbIs a state variable x and a state mean value muiDeviation of the b-th element, xcIs the c-th element of the state variable x, δ xcIs a state variable x and a state mean value muiDeviation of the c-th element. The higher odd-order central moments are similar and are all 0.
4.2) for the ith sub-Gaussian distribution, let
Figure GDA0002884026650000122
Wherein a is 1,2aIs the a-th element of the state variable x,
Figure GDA0002884026650000123
is the mean value of the states muiThe a-th element of (1); calculating a second central moment of the sub-Gaussian distribution deviation according to the formula (12);
E{δxaδxb}＝Pi ab,a,b∈[1,2,...,n] (12)
in the formula (12), E {. cndot } represents an expected value of a physical quantity in parentheses, E { δ x { [aδxbIs the second central moment of the sub-Gaussian distribution deviation, δ xaIs a state variable x and a state mean value muiDeviation of a-th element, δ xbIs a state variable x and a state mean value muiDeviation of the b-th element, Pi abRepresenting the covariance matrix P corresponding to the ith sub-Gaussian distributioniRow a, column b elements of (a);
4.3) on the basis of the second-order central moment of the known sub-Gaussian distribution deviation, making the m initial value be 2 and sequentially increasing by 2 as increment (m is 2,4, 6.. multidot.), and obtaining the m-th central moment and the covariance matrix PiDeducing and obtaining the m + 2-order central moment of sub-Gaussian distribution according to a formula (13) until the first 2N-order central moment of the sub-Gaussian distribution is calculated, wherein N is a deviation evolution order;
Figure GDA0002884026650000131
in the formula (13), E {. cndot } represents an expected value of a physical quantity in parentheses, E { δ x { }aδxbδxc...δxdThe m +2 central moments of sub-Gaussian distribution deviation(a, b, c., d) elements, a, b, c., d being the index of m +2 elements to the state variable x, (a, b, c., d) being the index of the elements to the m +2 central moments; pi abCovariance matrix P representing ith sub-Gaussian distributioniRow a, column b, elements E { δ xc...δxdThe (c, d) th element, which is the m-th order central moment of sub-gaussian distribution bias; pi acCovariance matrix P representing ith sub-Gaussian distributioniRow a, column c, element E { δ x [)b...δxdThe (b, d) th element, which is the m-th order central moment of sub-gaussian distribution bias; and so on until Pi adCovariance matrix P representing ith sub-Gaussian distributioniRow a and column d elements of (1), E { δ xbδxc... } is the (b, c.,) th element of the m-th order central moment of the sub-gaussian distribution bias.
In this embodiment, the detailed steps of step 5) include:
5.1) traversing and selecting one sub-Gaussian distribution as the current ith sub-Gaussian distribution;
5.2) for the ith sub-Gaussian distribution, the initial state is compared with the nominal value according to the formula (14)
Figure GDA0002884026650000133
Initial state deviation δ x of0Transformed into the mean μ of the ith sub-Gaussian distributioniAnd corresponding deviation δ xi；
δx0＝μi+δxi (14)
In the formula (14), δ x0Is the deviation of the initial state, muiIs the mean value of the ith sub-Gaussian distribution, δ xiIs the deviation of the ith sub-Gaussian distribution;
5.3) mean value μ to be transformed into ith sub-Gaussian distributioniAnd corresponding deviation δ xiInitial state deviation δ x of0Substituting the terminal state with respect to the initial state deviation deltax0Obtaining a terminal state shown in a formula (15) in the high-order Taylor expansion polynomial;
Figure GDA0002884026650000132
in the formula (15), xfFor the terminal state, T is a high-order Taylor expansion polynomial, the superscript N of T is a deviation evolution order, and the subscript x of TfAt the end state, δ xiIs the deviation of the ith sub-Gaussian distribution;
5.4) calculating the mean value mu of the ith sub-Gaussian distribution at the terminal moment according to the formula (16)i_fAs the expected value of the terminal state shown in equation (15);
Figure GDA0002884026650000141
in the formula (16), mui_fIs the mean value, x, of the ith sub-Gaussian distribution at the terminal timefFor the terminal state, T is a high-order Taylor expansion polynomial, the superscript N of T is a deviation evolution order, and the subscript x of TfAt the end state, δ xiIs the deviation of the ith sub-Gaussian distribution;
5.5) covariance matrix P of ith sub-Gaussian distribution of terminal time according to equation (17)i_fExpressed as deviation deltax with respect to the initial sub-gaussian distributioniThe 2N-order Taylor expansion and the expected value corresponding to the ith sub-Gaussian distribution are solved;
Figure GDA0002884026650000142
in the formula (17), Pi_fCovariance matrix, x, representing the ith sub-Gaussian distribution at the terminal timefIn the end state, μi_fThe mean value of the ith sub-Gaussian distribution at the terminal time, T is a high-order Taylor expansion polynomial, N in the superscript of T is a deviation evolution order, and deltaxiIs the deviation of the ith sub-Gaussian distribution;
5.6) judging whether all sub-Gaussian distributions are traversed or not, and if not, skipping to execute the step 5.2); otherwise, skipping to execute the next step;
5.7) obtaining all sub-Gaussian distribution at the terminalMean value of moment mui_fSum covariance matrix Pi_fOn the basis, an accurate probability density function of the system state distribution at the terminal moment is obtained according to the formula (18) and output;
Figure GDA0002884026650000143
in formula (18), p (t, x)f) An accurate probability density function representing the state distribution of the system at the terminal time, t being the terminal time, xfIs in the terminal state, wiThe weight corresponding to the ith sub-Gaussian distribution, M is the number of times of target shooting, pg(xf；μi_f,Pi_f) Is the probability density function of the ith sub-Gaussian distribution, mui_fIs the mean value of the ith sub-Gaussian distribution at the terminal time, Pi_fAnd the covariance matrix represents the ith sub-Gaussian distribution of the terminal moment.
Referring to a comparison graph of terminal state deviation distribution shown in fig. 3, sub-gaussian distribution represents the sub-gaussian distribution situation at the terminal time obtained by prediction of the present invention, MonteCarlo targeting represents the distribution situation of terminal time targeting points obtained by a conventional MonteCarlo simulation method, MonteCarlo covariance represents covariance 3-sigma ellipsoid counted according to MonteCarlo targeting points, MonteCarlo mean represents the distribution mean counted according to MonteCarlo targeting points, linearized covariance represents the terminal time covariance 3-sigma ellipsoid obtained by a conventional linearization method, and linearized mean represents the terminal time distribution mean obtained by a conventional linearization method. As can be seen from fig. 3, compared with the accurate MonteCarlo simulation, the terminal mean and covariance matrix predicted by the conventional linearization method have large errors, and the present embodiment well fits the distribution of MonteCarlo simulation sample points by using the gaussian sum model based on the nonlinear system state deviation evolution method of differential algebra and gaussian sum, and verifies the correctness and validity of the method.
Referring to the comparison graph of the terminal probability density function shown in fig. 4, the gaussian sum represents the terminal gaussian sum model probability density function obtained by the present invention, MonteCarlo represents the probability density function calculated according to MonteCarlo targeting, and the linearization method represents the terminal time probability density function obtained by the conventional linearization method. As can be further seen from fig. 4, the probability density function calculated by the nonlinear system state deviation evolution method based on differential algebra and sum of gaussians in the present embodiment is consistent with the MonteCarlo simulation statistical result, and the obvious probability density function no longer presents in the form of gaussian distribution, and meanwhile, it can be found that the traditional linearization method has obvious errors. Therefore, the method can accurately forecast the non-Gaussian deviation evolution condition of the strong nonlinear system.
In summary, the nonlinear system state deviation evolution method based on differential algebra and gaussian sum in this embodiment is a method for predicting the nonlinear system state deviation evolution process with high accuracy and high efficiency, based on the terminal state of the differential algebra prediction system, a recursive algorithm is applied to calculate the central moment of any order of gaussian distribution, which can be automatically expanded to the deviation evolution prediction of any order of precision, and gaussian sum model is used to describe the non-gaussian distribution characteristic of state deviation, so that the method can maintain the significant efficiency advantage relative to the MonteCarlo simulation method while accurately describing the deviation distribution and the non-gaussian characteristic thereof, and is suitable for strong nonlinear systems and long-time accurate prediction problems. The dynamic equation may be in any nonlinear mapping form or nonlinear very differential equation form, and the dynamic equation of this embodiment is a kepler orbit analytic solution in a nonlinear mapping form. The nonlinear system state deviation evolution method based on differential algebra and Gaussian sum in the embodiment is oriented to a general nonlinear dynamical system, can be automatically expanded to any high-order precision based on the differential algebra and Gaussian sum method, and realizes the precise evolution of the state deviation of the general nonlinear system.
Example two:
the process of this embodiment is basically the same as that of the first embodiment, and the main differences are as follows:
in step 2.1), the nonlinear system is a spacecraft nonlinear relative motion track of which the earth stationary orbit considers sunlight pressure, the system dimension n is 6, and the state variable x is a relative position and a velocity vector of the slave spacecraft in the orbit coordinate system of the master spacecraft, namely x is [ r ═ rT vT]TWhich isWherein r is [ x, y, z ]]TIs a relative position vector, v ═ vx,vy,vz]TIs a relative velocity vector.
In this embodiment, the kinetic equation of the nonlinear system is in the form of ordinary differential equation, and the terminal state is obtained by a variable-step length longge stota integration method.
In this embodiment, the specific expression of the dynamical equation of the nonlinear system is as follows:
Figure GDA0002884026650000161
where a is the major axis of the orbit of the main spacecraft, i.e. the radius a of the geostationary orbit of the earth is 42164200m in this embodiment, and the orbital angular velocity of the stationary orbit spacecraft is n 7.2921 × 10-5rad/s, constant of earth's gravity of 3.986 × 1014m3/s2The perturbation acceleration of sunlight pressure is gamma ═ gammax,γy,γz]TThe specific expression form can be referred to as the following documents: fehse, W., "Disturbance by Solar Pressure of Rendezvous traps in GEO," 5th ICATT, Noordwijk, Netherlands, May 2012.
The nominal value of the number of classical orbits of the spacecraft at the initial moment is given as follows:
Figure GDA0002884026650000162
given track forecast time tfIs one track period, i.e. t in this embodimentf＝86164s。
In step 2.1), an initial state covariance matrix P is given0Comprises the following steps:
Figure GDA0002884026650000163
compared with the first embodiment, the nonlinear system of the embodiment is changed into the nonlinear relative motion track of the spacecraft with the earth stationary orbit and the sunlight pressure considered, the kinetic equation is changed into a form of ordinary differential equation, and the terminal state is obtained by a variable-step length Runge Kutta integration method.
Referring to the comparison graph of terminal state deviation shown in fig. 5, sub-gaussian distribution represents the sub-gaussian distribution situation at the terminal time obtained by the prediction of the present invention, MonteCarlo targeting represents the distribution situation of terminal time targeting points obtained by the conventional MonteCarlo simulation method, MonteCarlo covariance represents the covariance 3-sigma ellipsoid counted according to MonteCarlo targeting points, MonteCarlo mean represents the distribution mean counted according to MonteCarlo targeting points, linearized covariance represents the terminal time covariance 3-sigma ellipsoid obtained by the conventional linearization method, and linearized mean represents the terminal time distribution mean obtained by the conventional linearization method. Referring to the comparison graph of the terminal probability density function shown in fig. 6, the gaussian sum represents the terminal gaussian sum model probability density function obtained by the present invention, MonteCarlo represents the probability density function calculated according to MonteCarlo targeting, and the linearization method represents the terminal time probability density function obtained by the conventional linearization method. Referring to the comparison between the terminal state deviation distribution shown in fig. 5 and the terminal probability density function shown in fig. 6, the nonlinear system state deviation evolution method based on differential algebra and gaussian sum in the embodiment can still well fit the MonteCarlo simulation result, effectively capture the non-gaussian property of the spacecraft relative trajectory deviation distribution, and accurately forecast the probability density function of the terminal time deviation distribution. The embodiment verifies the universality of the nonlinear system state deviation evolution method based on differential algebra and Gaussian sum, and for a general nonlinear system, the kinetic equation can be in a nonlinear mapping form or a normal differential equation form, and the evolution condition of the state deviation can be accurately forecasted.
The above description is only a preferred embodiment of the present invention, and the protection scope of the present invention is not limited to the above embodiments, and all technical solutions belonging to the idea of the present invention belong to the protection scope of the present invention. It should be noted that modifications and embellishments within the scope of the invention may occur to those skilled in the art without departing from the principle of the invention, and are considered to be within the scope of the invention.

Claims (6)

1. A nonlinear system state deviation evolution method based on differential algebra and Gaussian sum is used for spacecraft trajectory deviation evolution, and is characterized in that the method comprises the following steps:
1) inputting a kinetic equation, a deviation evolution order N and a relative nominal value of an initial state of a nonlinear system
Figure FDA0002884026640000011
And initial bias covariance matrix P0；
2) Forecasting the terminal state of the nonlinear system by adopting a differential algebra method according to the kinetic equation of the nonlinear system and expressing the terminal state as the deviation delta x of the initial state0A high order taylor expansion polynomial of (1);
3) for initial bias covariance matrix P0Determining a sub-Gaussian distribution covariance matrix PiFor sub-Gaussian distribution covariance matrix PiPerforming target fitting on the initial deviation distribution by the multiple sub-Gaussian distributions to generate a Gaussian sum model;
4) calculating a sub-Gaussian distribution high-order central moment of a given order aiming at the generated Gaussian sum model;
5) determining a mean value and a covariance matrix of each sub-Gaussian distribution at a terminal moment based on a high-order polynomial of the terminal state about the initial state deviation and each order central moment of the sub-Gaussian distribution, and further providing a terminal state deviation distribution probability density function in a Gaussian sum form;
the detailed steps of the step 4) comprise:
4.1) determining that any odd-order central moment of sub-Gaussian distribution is 0 according to the characteristics of Gaussian distribution;
4.2) for the ith sub-Gaussian distribution, let
Figure FDA0002884026640000012
Wherein a is 1,2aIs the a-th element of the state variable x,
Figure FDA0002884026640000013
is the mean value of the states muiThe a-th element of (1); calculating a second central moment of the sub-Gaussian distribution deviation according to the formula (12);
E{δxaδxb}＝Pi ab,a,b∈[1,2,...,n] (12)
in the formula (12), E {. cndot } represents an expected value of a physical quantity in parentheses, E { δ x { [aδxbIs the second central moment of the sub-Gaussian distribution deviation, δ xaIs a state variable x and a state mean value muiDeviation of a-th element, δ xbIs a state variable x and a state mean value muiDeviation of the b-th element, Pi abRepresenting the covariance matrix P corresponding to the ith sub-Gaussian distributioniRow a, column b elements of (a);
4.3) on the basis of the second-order central moment of the known sub-Gaussian distribution deviation, making the m initial value be 2 and sequentially increasing by taking 2 as increment, and obtaining the m-th-order central moment and the covariance matrix PiDeducing and obtaining the m + 2-order central moment of sub-Gaussian distribution according to a formula (13) until the first 2N-order central moment of the sub-Gaussian distribution is calculated, wherein N is a deviation evolution order;
Figure FDA0002884026640000014
in the formula (13), E {. cndot } represents an expected value of a physical quantity in parentheses, E { δ x { }aδxbδxc...δxdThe (a, b, c.., d) th element which is the m +2 order central moment of the sub-gaussian distribution deviation, (a, b, c.., d) is the index to the element of the m +2 order central moment; pi abCovariance matrix P representing ith sub-Gaussian distributioniRow a, column b, elements E { δ xc...δxdThe (c, d) th element, which is the m-th order central moment of sub-gaussian distribution bias; pi acCovariance matrix P representing ith sub-Gaussian distributioniRow a, column c, element E { δ x [)b...δxdThe (b, d) th element, which is the m-th order central moment of sub-gaussian distribution bias; and so on until Pi adCovariance matrix P representing ith sub-Gaussian distributioniRow a and column d elements of (1), E { δ xbδxc... } is the (b, c.,) th element of the m-th order central moment of the sub-gaussian distribution bias;
the detailed steps of the step 5) comprise:
5.1) traversing and selecting one sub-Gaussian distribution as the current ith sub-Gaussian distribution;
5.2) for the ith sub-Gaussian distribution, the initial state is compared with the nominal value according to the formula (14)
Figure FDA0002884026640000023
Initial state deviation δ x of0Transformed into the mean μ of the ith sub-Gaussian distributioniAnd corresponding deviation δ xi；
δx0＝μi+δxi (14)
In the formula (14), δ x0Is the deviation of the initial state, muiIs the mean value of the ith sub-Gaussian distribution, δ xiIs the deviation of the ith sub-Gaussian distribution;
5.3) mean value μ to be transformed into ith sub-Gaussian distributioniAnd corresponding deviation δ xiInitial state deviation δ x of0Substituting the terminal state with respect to the initial state deviation deltax0Obtaining a terminal state shown in a formula (15) in the high-order Taylor expansion polynomial;
Figure FDA0002884026640000021
in the formula (15), xfFor the terminal state, T is a high-order Taylor expansion polynomial, the superscript N of T is a deviation evolution order, and the subscript x of TfAt the end state, δ xiIs the deviation of the ith sub-Gaussian distribution;
5.4) calculating the mean value mu of the ith sub-Gaussian distribution at the terminal moment according to the formula (16)i_fAs the expected value of the terminal state shown in equation (15);
Figure FDA0002884026640000022
in the formula (16), mui_fIs the mean value, x, of the ith sub-Gaussian distribution at the terminal timefFor the terminal state, T is a high-order Taylor expansion polynomial, the superscript N of T is a deviation evolution order, and the subscript x of TfAt the end state, δ xiIs the deviation of the ith sub-Gaussian distribution;
5.5) covariance matrix P of ith sub-Gaussian distribution of terminal time according to equation (17)i_fExpressed as deviation deltax with respect to the initial sub-gaussian distributioniThe 2N-order Taylor expansion and the expected value corresponding to the ith sub-Gaussian distribution are solved;
Figure FDA0002884026640000031
in the formula (17), Pi_fCovariance matrix, x, representing the ith sub-Gaussian distribution at the terminal timefIn the end state, μi_fThe mean value of the ith sub-Gaussian distribution at the terminal time, T is a high-order Taylor expansion polynomial, N in the superscript of T is a deviation evolution order, and deltaxiIs the deviation of the ith sub-Gaussian distribution;
5.6) judging whether all sub-Gaussian distributions are traversed or not, and if not, skipping to execute the step 5.2); otherwise, skipping to execute the next step;
5.7) obtaining the mean value mu of all sub-Gaussian distributions at the terminal momenti_fSum covariance matrix Pi_fOn the basis, an accurate probability density function of the system state distribution at the terminal moment is obtained according to the formula (18) and output;
Figure FDA0002884026640000032
in formula (18), p (t, x)f) An accurate probability density function representing the state distribution of the system at the terminal time, t being the terminal time, xfIs in the terminal state, wiIs the ithWeight corresponding to sub-Gaussian distribution, M is the number of times of target shooting, pg(xf；μi_f,Pi_f) Is the probability density function of the ith sub-Gaussian distribution, mui_fIs the mean value of the ith sub-Gaussian distribution at the terminal time, Pi_fAnd the covariance matrix represents the ith sub-Gaussian distribution of the terminal moment.
2. The differential algebra and Gaussian sum based nonlinear system state deviation evolution method as claimed in claim 1, wherein the dynamical equation of the nonlinear system in the step 1) is a nonlinear mapping form xf＝f(x0) Or form of ordinary differential equation
Figure FDA0002884026640000033
And the function f (-) in the form of the nonlinear mapping or the function g (-) in the form of the ordinary differential equation is an arbitrary n-dimensional real number domain RnTo RnNon-linear mapping of (2), where x0Is in an initial state, xfIs the terminal state, x is the state variable.
3. The differential algebra and gaussian sum based nonlinear system state deviation evolution method as claimed in claim 1, wherein the detailed steps of step 2) comprise:
2.1) initial state x0The initialization is represented by formula (1) and includes the relative nominal value of the initial state
Figure FDA0002884026640000034
The form of deviation of (1);
Figure FDA0002884026640000035
in the formula (1), x0In the initial state, the state of the device is as follows,
Figure FDA0002884026640000036
relative to nominal value for initial state, δ x0Is an initial state deviation;
2.2) initial state x0Substituting into the dynamic equation of the nonlinear system, and applying a differential algebra method to obtain the deviation deltax about the initial state0The terminal state described by the high-order Taylor expansion polynomial is shown as the formula (2);
Figure FDA0002884026640000037
in the formula (2), xfFor the terminal state, T is a high-order Taylor expansion polynomial, the superscript N of T is a deviation evolution order, and the subscript x of TfThe contents of the polynomial description for the higher order Taylor expansion are terminal states, δ x0Is the initial state deviation.
4. The differential algebra and gaussian sum based nonlinear system state deviation evolution method as claimed in claim 1, wherein the detailed step of step 3) comprises:
3.1) covariance matrix P for initial bias0Inverse of (3) performing Cholesky decomposition to obtain square root information matrix R0And satisfy P0＝R0 -1R0 -T；
3.2) lower bound R of the square root information matrix given a Gaussian distribution of the childrenminAnd obtaining a sub-Gaussian distribution covariance matrix Pi；
3.3) target multiple sub-Gaussian distributions to generate Gaussian sums to fit the initial Gaussian distribution.
5. The method for evolving state deviation of nonlinear system based on differential algebra and Gaussian sum as claimed in claim 4, wherein the sub-Gaussian distribution covariance matrix P in step 3.2) is obtainediThe detailed steps comprise:
3.2.1) computing the matrix using the functional expression shown in equation (3)
Figure FDA0002884026640000041
Singular value decomposition of (1), wherein R0Is square root informationMatrix, RminThe lower bound of the square root information matrix of sub-Gaussian distribution;
Figure FDA0002884026640000042
in the formula (3), UiAnd ViAs orthonormal, diagonal, Si＝diag(σi1,...,σin) Wherein the diagonal element σi1,...,σinIs a positive singular value, R0Is a square root information matrix, RminThe lower bound of the square root information matrix of sub-Gaussian distribution;
3.2.2) determining the diagonal matrix S for any k ═ 1iOf (a) the kth diagonal element σikIf 1 or more is true, the square root information matrix R of the sub-Gaussian distribution is obtainediIs a square root information matrix R0Jump to execute step 3.2.5); otherwise, skipping to execute the next step;
3.2.3) establishing an n-dimensional diagonal matrix shown in the vertical type (4), deleting all 0 rows of the n-dimensional diagonal matrix to obtain a singular value change matrix delta Si；
Figure FDA0002884026640000043
In the formula (4), δ Si_fullFor an n-dimensional diagonal matrix, σi1Is a diagonal matrix SiThe 1 st diagonal element of (1), σinIs a diagonal matrix SiThe nth diagonal element of (1);
3.2.4) determining the information variation matrix δ R according to equation (5)iAnd solving a square root information matrix R of sub-Gaussian distribution according to the formula (6)i；
δRi＝δSiVi TRmin (5)
In the formula (5), δ RiFor information change matrix, δ SiAs a matrix of singular value variations, ViIs an orthonormal matrix, RminIs a height ofThe lower bound of the square root information matrix of the gaussian distribution;
Figure FDA0002884026640000051
in the formula (6), QiIs an orthonormal matrix, RiIs a square root information matrix of sub-Gaussian distribution, δ RiFor information change matrix, R0Is a square root information matrix;
3.2.5) calculating a reduction of the sub-Gaussian covariance compared to the initial Gaussian covariance according to equation (7);
δPi＝P0-Pi＝R0 -1R0 -T-Ri -1Ri -T＝δYi(δYi)T (7)
in the formula (7), δ PiIs the amount of reduction of the sub-Gaussian distribution covariance compared to the initial Gaussian distribution covariance, δ YiIs the square root of the covariance reduction, deltaPiSemi-positive and delta YiIs δ PiSquare root of matrix of (d), δ YiThe functional expression of (a) is represented by the formula (8);
Figure FDA0002884026640000052
in the formula (8), δ YiIs the square root of the covariance reduction, R0Is a square root information matrix, δ RiFor information change matrix, matrix RdiObtained by QR decomposition as shown in a formula (9);
Figure FDA0002884026640000053
in the formula (9), QdiIs an orthonormal matrix, RdiIs an upper triangular square matrix, R0Is a square root information matrix, δ RiIs an information change matrix, and I is an n-dimensional unit matrix;
3.2.6) covariance matrix P based on initial bias0The amount of reduction of the sub-Gaussian distribution covariance compared to the initial Gaussian distribution covariance, δ PiDetermining a sub-Gaussian distribution covariance matrix Pi。
6. The differential algebra and gaussian sum based nonlinear system state deviation evolution method as claimed in claim 5, wherein the detailed step of step 3.3) comprises:
3.3.1) note nδYMatrix δ Y being the square root of covariance reductioniColumn number of (d), singular value change matrix δ SiSetting the number of target shooting times M as the number of sub-Gaussian distributions, and initializing i;
3.3.2) recording the current target shooting frequency as the ith time;
3.3.3) generating nδYVector of dimension ηiThe components thereof obey standard Gaussian distribution and are independent of each other;
3.3.4) according to μi＝μ0+δYiηiCalculating the mean value mu of the ith sub-Gaussian distributioniIn which μ0Is the mean of the initial states, δ YiIs the square root, η, of the covariance reductioniIs nδYA dimension column vector;
3.3.5) determining the square root information matrix R corresponding to the ith sub-Gaussian distributioniAnd a covariance matrix Pi；
3.3.6) determining the weight corresponding to the ith sub-Gaussian distribution as wi1/M, wherein M is the number of times of target shooting;
3.3.7) adding 1 to the variable i, judging whether i is smaller than N, and if so, skipping to execute the step 3.3.2); otherwise, skipping to execute the next step;
3.3.8) generating an initial deviation probability density function of the form of a sum of gaussians as shown in equation (10);
Figure FDA0002884026640000061
in the formula (10), p (t, x)0) In the form of a sum of gaussiansOf the initial deviation probability density function, wiThe weight corresponding to the ith sub-Gaussian distribution, M is the number of times of target shooting, pg(x0；μi,Pi) Is the probability density function of the ith sub-Gaussian distribution and the probability density function p of the ith sub-Gaussian distributiong(x0；μi,Pi) The functional expression of (a) is represented by the formula (11);
Figure FDA0002884026640000062
in the formula (11), n is the system dimension, PiCovariance matrix, x, corresponding to ith sub-Gaussian distribution0Is in an initial state, muiIs the mean of the ith sub-gaussian distribution.
CN201710551107.9A 2017-07-07 2017-07-07 Nonlinear system state deviation evolution method based on differential algebra and Gaussian sum Active CN107402903B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710551107.9A CN107402903B (en) 2017-07-07 2017-07-07 Nonlinear system state deviation evolution method based on differential algebra and Gaussian sum

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710551107.9A CN107402903B (en) 2017-07-07 2017-07-07 Nonlinear system state deviation evolution method based on differential algebra and Gaussian sum

Publications (2)

Publication Number Publication Date
CN107402903A CN107402903A (en) 2017-11-28
CN107402903B true CN107402903B (en) 2021-02-26

Family

ID=60404964

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710551107.9A Active CN107402903B (en) 2017-07-07 2017-07-07 Nonlinear system state deviation evolution method based on differential algebra and Gaussian sum

Country Status (1)

Country Link
CN (1) CN107402903B (en)

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108490966B (en) * 2018-01-31 2021-02-05 中国人民解放军国防科技大学 Static orbit perturbation relative trajectory high-order guidance method based on differential algebra
CN109032176B (en) * 2018-07-25 2021-06-22 西北工业大学 Geosynchronous orbit determination and parameter determination method based on differential algebra
CN109508482A (en) * 2018-10-26 2019-03-22 天津大学 A kind of calculation method for complex-curved surface profile degree error uncertainty
CN110442831B (en) * 2019-07-31 2023-03-24 中国人民解放军国防科技大学 Space non-cooperative target space-based search method based on nonlinear deviation evolution
CN111899297B (en) * 2020-08-06 2024-01-23 中国铁建重工集团股份有限公司 Method for extracting center of light stripe of line structure
CN113836697B (en) * 2021-08-26 2024-01-23 北京理工大学 Spatial gravitational wave detection configuration stability evolution method based on dynamic Gaussian mixture
CN114329887A (en) * 2021-11-22 2022-04-12 北京理工大学 Data-driven dynamic Gaussian mixture spacecraft orbit error evolution method

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106697333A (en) * 2017-01-12 2017-05-24 北京理工大学 Robustness analysis method for spacecraft orbit control strategy

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101770024B (en) * 2010-01-25 2011-08-31 上海交通大学 Multi-target tracking method
CN102040008A (en) * 2010-12-13 2011-05-04 北京航空航天大学 Anti-collision control method for safety of in-obit operation of formation-flying satellites
CN105303052A (en) * 2015-11-11 2016-02-03 中国人民解放军国防科学技术大学 Low-speed approaching spacecraft track safety assessment method

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106697333A (en) * 2017-01-12 2017-05-24 北京理工大学 Robustness analysis method for spacecraft orbit control strategy

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Nonlinear Mapping of Uncertainties in Celestial Mechanics;M. Valli et al.;《Journal of Guidance Control and Dynamics》;20130228;第36卷(第1期);第48-63页 *
航天器轨迹偏差的非线性非高斯预报方法;孙振江 等;《宇航学报》;20161130;第37卷(第11期);第1289-1297页 *

Also Published As

Publication number Publication date
CN107402903A (en) 2017-11-28

Similar Documents

Publication Publication Date Title
CN107402903B (en) Nonlinear system state deviation evolution method based on differential algebra and Gaussian sum
Gross et al. A comparison of extended Kalman filter, sigma-point Kalman filter, and particle filter in GPS/INS sensor fusion
da Silva et al. Ensemble-based state estimator for aerodynamic flows
CN106970643B (en) Analytic satellite nonlinear relative motion deviation propagation analysis method
CN110378012B (en) Strict regression orbit design method, system and medium considering high-order gravity field
CN110059439B (en) Spacecraft orbit determination method based on data driving
CN111474960A (en) Critical altitude supersonic speed target tracking method based on control quantity characteristic assistance
CN114936471A (en) Spacecraft collision early warning layered rapid screening method based on parallel computing
Huang et al. Reliability-based trajectory optimization using nonintrusive polynomial chaos for Mars entry mission
Xue et al. A comparison of several nonlinear filters for mobile robot pose estimation
Brock et al. Dynamic CFD simulations of the supersonic inflatable aerodynamic decelerator (SAID) ballistic range tests
CN110231619B (en) Radar handover time forecasting method and device based on Enk method
Aeschliman et al. A proposed methodology for computational fluid dynamics code verification, calibration, and validation
CN114462256B (en) Method, device, equipment and medium for determining non-cooperative low-thrust maneuvering target track
CN111428912A (en) Mars detector orbit prediction method and system based on support vector machine
CN111259337A (en) Heavy debris real-time drop point forecasting method based on statistics
Zhou et al. Inverse simulation system for manual-controlled rendezvous and docking based on artificial neural network
Huang et al. Uncertainty analysis of reachable set for planetary entry using polynomial chaos
Morris et al. Real-time system identification of quadrotor dynamics
CN112327333B (en) Satellite position calculation method and device
Qi et al. Contact stiffness and damping identification for hardware-in-the-loop contact simulator with measurement delay compensation
CN109709805B (en) Spacecraft robust intersection trajectory design method considering uncertainty factors
Li et al. Efficient and accurate error propagation in the semi-analytic orbit dynamics system for space debris
Qi et al. Hybrid Smith predictor and phase lead based divergence compensation for hardware-in-the-loop contact simulation with measurement delay
CN109583007B (en) Mars entering flight state uncertainty quantification method

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant