CN109343550B - Spacecraft angular velocity estimation method based on rolling time domain estimation - Google Patents

Spacecraft angular velocity estimation method based on rolling time domain estimation Download PDF

Info

Publication number
CN109343550B
CN109343550B CN201811193664.9A CN201811193664A CN109343550B CN 109343550 B CN109343550 B CN 109343550B CN 201811193664 A CN201811193664 A CN 201811193664A CN 109343550 B CN109343550 B CN 109343550B
Authority
CN
China
Prior art keywords
spacecraft
angular velocity
estimation
constraint
time domain
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
CN201811193664.9A
Other languages
Chinese (zh)
Other versions
CN109343550A (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.)
Beihang University
Original Assignee
Beihang University
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 Beihang University filed Critical Beihang University
Priority to CN201811193664.9A priority Critical patent/CN109343550B/en
Publication of CN109343550A publication Critical patent/CN109343550A/en
Application granted granted Critical
Publication of CN109343550B publication Critical patent/CN109343550B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Abstract

The invention relates to a rolling time domain estimation-based spacecraft angular velocity estimation method, which comprises the following steps of: constructing a state equation and an observation equation of the angular velocity of the spacecraft based on a spacecraft dynamics model; simplifying four paths of twelve-dimensional angular velocity measurement information of the redundant fiber optic gyroscope combination to obtain three-axis angular velocity measurement data of the spacecraft; based on the constructed state equation and observation equation and the obtained spacecraft triaxial angular velocity measurement data, spacecraft angular velocity constraint and control moment constraint are considered at the same time, online estimation is carried out by using a rolling time domain estimation strategy and a multiple target shooting method, and the estimated value of the spacecraft optimal angular velocity at the current moment is obtained. The method has the advantages of high calculation efficiency, high reliability, good operation effect, high precision and the like, simultaneously fully considers the influences of the physical constraint and the control moment constraint of the spacecraft, is an effective optimal online state estimation method, and is suitable for the angular velocity estimation problem of the constrained spacecraft.

Description

Spacecraft angular velocity estimation method based on rolling time domain estimation
Technical Field
The invention relates to a rolling time domain estimation-based spacecraft angular velocity estimation method, which is mainly applied to the problem of spacecraft angular velocity on-line estimation considering angular velocity constraint and control moment constraint and belongs to the field of spacecraft attitude determination.
Background
The development of the aerospace technology is based on the control and navigation technology of the spacecraft, and how to improve the control and navigation precision of the spacecraft on the basis of the existing system hardware equipment has important significance. In recent years, with the trend of diversification and complication of flight missions of spacecrafts, for the spacecrafts performing the flight missions in orbit, the accuracy of determining attitude parameters of the spacecrafts determines the accuracy of attitude control, and further influences the function realization of the whole flight missions.
The emergence of new complex space missions puts higher requirements on the control precision of the spacecraft, so that the research of a high-precision angular velocity estimation method is particularly important. Most of existing spacecraft angular velocity estimation methods are researched around Kalman filtering algorithms, Kalman filtering has many advantages, meanwhile, conclusions about linear Kalman filtering are very mature, but Kalman filtering cannot do constraint processing and system model uncertainty, so that very useful information existing in reality in an actual system is ignored, and the practicability of Kalman filtering is limited. The rolling time domain estimation method is designed based on a model, only a fixed amount of measurement data before the current moment is considered, and the influence of the measurement data at other moments on estimation is described by an approximate method, so that the problems of 'infinite memory' and 'data explosion' are well solved, known information about system interference and system state in a time domain constraint form can be fully utilized, the time domain constraint and nonlinear characteristics of the system are directly expressed in an optimization problem, and the system state is accurately estimated, and the reasonability and accuracy of estimation are improved.
For the problem of spacecraft angular velocity estimation, patent CN201711456224.3 selects a Singer acceleration probability model taking spacecraft angular velocity and angular acceleration as state quantities as a state model, performs image processing on a multi-frame fuzzy star map, extracts optical flow information as quantity measurement, establishes a measurement model according to the relation between star point optical flow and spacecraft angular velocity, and finally completes estimation of spacecraft triaxial angular velocity through kalman filtering, but the state model, the measurement model and the estimation algorithm in the method are all different from the method; the patent CN201610180095.9 provides a method for estimating the constrained angular velocity of a gyroscope-free inertial measurement system, which has the characteristics of clear estimation steps and estimation error which does not change along with time, can be used for estimating the angular velocity under the condition of configuring a redundant accelerometer of the gyroscope-free inertial measurement system, and can effectively improve the estimation precision of the angular velocity, but the method does not consider the constraint of the angular velocity;
for an application example of a rolling time domain estimation algorithm, patent CN201210102226.3 firstly identifies parameters of a power battery model, and then estimates the state of charge of the power battery by using a rolling time domain estimation method on the basis of the determined model; the patent CN201610201135.3 estimates the wave frequency model parameters and wave frequency motion of the dynamic positioning ship by using a rolling time domain estimation method; the patent CN201510235759.2 considers the state constraint condition for a linear model in radar tracking, and adopts a linear rolling time domain method to estimate the position of a tracking target; the patent CN201711466618.7 firstly carries out sensitivity analysis based on a circuit parameter function to establish an augmented nonlinear state space equation, then establishes an SOC and model parameter adaptive joint estimation model based on a rolling time domain estimation strategy, sets each algorithm parameter, and finally carries out SOC and model parameter joint estimation by utilizing the online estimation model based on detection voltage and current. The system models used in the above patents are not spacecraft models, so the above methods are not suitable for spacecraft.
Disclosure of Invention
The technical problem to be solved by the invention is as follows: because the spacecraft is often influenced by various external disturbance moments such as gravity gradient moment, aerodynamic resistance, sunlight pressure and the like in the process of attitude control in space, a spacecraft dynamics model has certain uncertainty, and meanwhile, a measurement error exists in the measurement of the angular speed of an optical fiber gyroscope on the spacecraft. In order to obtain effective angular velocity information with higher precision, the invention provides a rolling time domain estimation-based spacecraft angular velocity estimation method, which is a state estimation method which is designed based on a model, considers state time domain constraint, considers fixed time domain measurement data, is reasonable and has higher estimation precision.
According to an aspect of the invention, a method for estimating the angular velocity of a spacecraft based on rolling time domain estimation is provided, which comprises the following steps:
s1: constructing a state equation and an observation equation of the angular velocity of the spacecraft based on a spacecraft dynamics model in a special orthogonal group SO (3) description mode;
s2: simplifying four-way twelve-dimensional angular rate measurement information combined by redundant fiber optic gyroscopes by using a new weighted data fusion method for distributing weights based on failure factors to obtain spacecraft three-axis angular rate measurement data;
s3: based on the state equation and the observation equation constructed in the step S1 and the spacecraft triaxial angular rate measurement data obtained in the step S2, spacecraft angular rate constraint and control moment constraint are considered at the same time, online estimation is carried out by using a rolling time domain estimation strategy and a multiple target shooting method, and an estimated value of the spacecraft optimal angular rate at the current moment is obtained.
Further, the spacecraft dynamics model in the special orthogonal group SO (3) description mode in step S1 is:
Figure GDA0001934996820000031
Rk+1＝RkFk(2)
Figure GDA00019349968200000314
wherein h is a sampling time interval; superscript T represents matrix transposition; the index k denotes the kth sampling instant; matrix Rk,FkE SO (3) is a rotation matrix,
Figure GDA0001934996820000032
a rotation matrix representing the euler attitude angle description of the spacecraft,
Figure GDA0001934996820000033
a rotation matrix representing the description of the Euler attitude error angle of the spacecraft between two sampling instants, FkIs the solution of implicit equation (1);
Figure GDA0001934996820000034
which represents the angular momentum of the spacecraft,
Figure GDA0001934996820000035
is an inertial matrix of the spacecraft,
Figure GDA0001934996820000036
is the spacecraft angular velocity;
Figure GDA0001934996820000037
representing the external disturbance moment of the spacecraft caused by gravity gradient moment, aerodynamic resistance, sunlight pressure and the like；
Figure GDA0001934996820000038
Representing the control moment required by the attitude maneuver or the stability of the spacecraft; positive definite array
Figure GDA0001934996820000039
Represents a non-standard inertia matrix, which is related to the inertia matrix J;
Figure GDA00019349968200000310
is about real vector
Figure GDA00019349968200000311
For any vector
Figure GDA00019349968200000312
Its corresponding oblique antisymmetric array x×The form is as follows:
Figure GDA00019349968200000313
further, based on the spacecraft dynamics model, a state equation and an observation equation for rolling time domain estimation can be constructed as follows:
taking the system state quantity as omega ═ omegax,ωy,ωz)TThen the state equation is
Figure GDA0001934996820000041
Taking the observed quantity of the system as y ═ yx,yy,yz)TThen the observation equation is
Jyk＝CJωk+νk(5)
Wherein ω ═ ω (ω ═ ω)x,ωy,ωz)TIs the spacecraft angular velocity; y ═ yx,yy,yz)TThe output value of the angular speed measurement data after the redundant fiber optic gyroscope combined data processing,yx,yy,yzmeasuring the angular rate of each inertia main shaft of the spacecraft; fkThe epsilon SO (3) is an auxiliary variable and represents the relative attitude between two adjacent moments;
Figure GDA0001934996820000042
and
Figure GDA0001934996820000043
the system external interference and the measurement noise are respectively independent from each other, and the measurement noise is independent from the system external interference and the initial state;
Figure GDA0001934996820000044
is a measurement matrix.
Furthermore, since the observation equation needs three-dimensional measurement data, and the combination of the redundant fiber optic gyroscopes feels that the combination of the fiber optic gyroscopes is four-way twelve-dimensional measurement data, a new data fusion method for distributing weights based on failure factors needs to be designed. The new weighted data fusion method for distributing the weights based on the failure factors in step S2 is as follows:
Figure GDA0001934996820000045
wherein the content of the first and second substances,
Figure GDA0001934996820000046
measuring data of the angular velocity of the spacecraft after data processing;
Figure GDA0001934996820000047
the angular velocity information of the spacecraft sensed by the ith optical fiber gyroscope in the redundant optical fiber gyroscopes,
Figure GDA0001934996820000048
for four-way twelve-dimensional angular rate measurement information αiAnd i is 1,2,3 and 4 is a mounting deflection angle of the optical fiber gyroscope, and represents the projection of the axial direction of the ith optical fiber gyroscope on the xoy plane and the ox axisβiThe i is 1,2,3 and 4, which is the installation high-low angle of the optical fiber gyroscope and represents the included angle between the axial direction of the ith optical fiber gyroscope and the xoy plane;
Figure GDA0001934996820000051
as a weighting factor, the failure factor ρ of the ith optical fiber gyroscopei∈[0,1]Influence and need to satisfy the constraint τ1+τ2+τ3+τ4＝1。
Further, the specific steps of performing online estimation by using the rolling time domain estimation strategy and the multiple targeting method in step S3 are as follows:
1) a priori estimation of initial state assuming a known system
Figure GDA0001934996820000052
Covariance matrix of initial state error is P0The weight matrix of the external interference and the measured noise of the system is Q1,Q2And P is0,Q1,Q2All are positive definite matrixes, and the length of a rolling time domain is N;
2) at time T, the optimization problem defined by the rolling time domain estimation is as follows:
Figure GDA0001934996820000053
in the two cases, the first term in the objective function represents the cost sum of external interference of the system at each sampling moment, the second term represents the cost sum of measurement noise at each sampling moment, and the third term represents prior estimation when T is less than or equal to N
Figure GDA0001934996820000054
Degree of trust of, at T>When the N is greater than the N value,
Figure GDA0001934996820000055
to approximate the cost function of a constrained nonlinear system, i.e., to the arrival cost function of an unconstrained linear system, representing the effect of the remaining measurement data on the estimate, the cost function of a constrained nonlinear system is approximated by the arrival cost function of an unconstrained linear systemThe cost function is of the form:
Figure GDA0001934996820000056
in the formula, initial state estimation
Figure GDA0001934996820000057
Obey normal distribution
Figure GDA0001934996820000058
PT-NTo estimate the error covariance, it can be updated by the following equation (8);
Figure GDA0001934996820000059
meanwhile, the objective function needs to satisfy spacecraft dynamics, spacecraft angular velocity time domain constraint and control moment constraint, and the specific mathematical expression form is as follows:
angular velocity time domain constraint:
Figure GDA0001934996820000061
wherein the content of the first and second substances,
Figure GDA0001934996820000062
for the angular rate estimate obtained by the rolling time domain estimation,
Figure GDA0001934996820000063
the minimum and maximum angular velocities allowed by the main axes of the spacecraft are determined by the properties of physical devices such as solar sailboards, sensors and the like installed on the spacecraft. It should be noted that the angular velocity time domain constraint proposed herein is the optimal angular velocity estimated by the rolling time domain online estimation algorithm
Figure GDA0001934996820000064
Constraining rather than imposing an actual operating angular velocity ω of the spacecraftkGo on to aboutBundles, require special attention to distinguish in order to avoid ambiguity.
Controlling torque constraint:
Figure GDA0001934996820000065
and because the spacecraft system is in a discrete system form, the angular acceleration
Figure GDA0001934996820000066
Can be approximated in the following form:
Figure GDA0001934996820000067
the control torque constraint can be written approximately as:
Figure GDA0001934996820000068
simplifying the above formula, can convert control moment restraint into the form of angular velocity restraint, the concrete form is as follows:
Figure GDA0001934996820000069
wherein the content of the first and second substances,
Figure GDA00019349968200000610
respectively providing minimum and maximum control moments for the spacecraft actuating mechanism in the directions of all main shafts;
Figure GDA00019349968200000611
the optimal estimated value of the angular rate of the spacecraft at the previous moment is obtained; t issThe time interval is sampled for a discrete system. Similarly, it should be noted that the constraint equation (13) is the optimal angular velocity estimated by the rolling time domain online estimation algorithm
Figure GDA00019349968200000612
For constraining, not for spacecraftActual operating angular velocity ωkAnd (6) carrying out constraint.
Combining the angular velocity constraint (9) and the control moment constraint (13), two constraints can be simplified as:
Figure GDA00019349968200000613
Figure GDA00019349968200000614
Figure GDA0001934996820000071
at this time, the constrained rolling horizon optimization problem can be written as follows:
Figure GDA0001934996820000072
Figure GDA0001934996820000073
Figure GDA0001934996820000074
Figure GDA0001934996820000075
3) at the time T, if the redundant fiber optic gyroscope after data processing is combined with the angular speed measurement output sequence
Figure GDA0001934996820000076
Or
Figure GDA0001934996820000077
And spacecraft attitude control moment sequence
Figure GDA0001934996820000078
Or
Figure GDA0001934996820000079
Known, therein for arbitrary variables
Figure GDA00019349968200000710
Sequence of
Figure GDA00019349968200000711
To represent
Figure GDA00019349968200000712
The optimization problem (15) can be solved using a multiple targeting method to obtain a unique existing optimal solution
Figure GDA00019349968200000713
4) At time T, the specific steps for solving the optimization problem (15) by using the multiple targeting method are as follows:
firstly, the angular momentum pi of the spacecraft is used herek＝JωkReplacing spacecraft angular velocity omegakThe optimization problem (15) can be equivalently transformed into the form of the following optimization problem (16):
Figure GDA00019349968200000714
Figure GDA00019349968200000715
Figure GDA00019349968200000716
Figure GDA00019349968200000717
wherein, Yk＝JykFrom the spacecraft angular velocity measurement y for time kkCalculating the angular momentum of the spacecraft;
Figure GDA0001934996820000081
3, estimating values of angular momentum of each main shaft of the spacecraft at the moment k;
Figure GDA0001934996820000082
Figure GDA0001934996820000083
respectively the minimum and maximum angular momentum of each main shaft of the spacecraft; equation (16c) is solving for the auxiliary variables
Figure GDA0001934996820000084
An implicit equation of (c);
considering system dynamics equality constraint (16b) and angular momentum time domain inequality constraint (16d), defining an augmented objective function:
Figure GDA0001934996820000085
wherein the content of the first and second substances,
Figure GDA0001934996820000086
is a Lagrangian multiplier,. mu. > 0 and
Figure GDA0001934996820000087
in order to be a penalty parameter,
Figure GDA0001934996820000088
the penalty function for the corresponding inequality constraint (16d) is expressed in the following way:
Figure GDA0001934996820000089
when T is larger than N, the broadening cost function is subjected to variation by using a variation method, wherein the variation expression is as follows:
Figure GDA00019349968200000810
wherein the content of the first and second substances,
Figure GDA00019349968200000811
as a penalty function
Figure GDA00019349968200000812
To pair
Figure GDA00019349968200000813
First order variational sum, variational coefficient
Figure GDA00019349968200000814
As follows:
Figure GDA00019349968200000815
wherein e isiI is 1,2 and 3 are unit vectors of main axes of inertia of the spacecraft respectively;
let δ ΦTaThe first order requirement for the 0 available optimization problem (16) is as follows:
Figure GDA0001934996820000091
wherein the content of the first and second substances,
Figure GDA0001934996820000092
the following first order requirements can be met by modifying equation (21):
Figure GDA0001934996820000093
wherein the content of the first and second substances,
Figure GDA0001934996820000094
for the vectors to be solved of the first order requirement, by lagrange multiplier sequence
Figure GDA0001934996820000095
And spacecraft angular momentum estimation sequence
Figure GDA0001934996820000096
Composition is carried out;
at this time, solving the optimization problem (16) can be equivalently converted into solving the root of the first-order requirement (22), namely, the root of the first-order requirement and the optimal solution of the optimization problem are the same, and the numerical iteration solving is carried out on the formula (22) by adopting a multiple targeting method, namely, a group of specific sequences s is required to be found*So that Fs(s*)＝06N+3If the first order requirement (22) has a root and is unique, the root is the optimal solution of the optimization problem (16);
the root-finding iterative relationship in the multiple targeting method is as follows:
si+1＝si+κΔ(s)i(23)
wherein s isi+1For the i +1 th iterative solution, siFor the ith iterative solution,. kappa.. epsilon. (0, 1)]In order to iterate the step size,
Figure GDA0001934996820000101
for iterative descending direction, Δ(s)iIn (1)
Figure GDA0001934996820000102
A Jacobian matrix of formula (22);
a variational expression according to equation (22):
Figure GDA0001934996820000103
Figure GDA0001934996820000104
Figure GDA0001934996820000105
Figure GDA0001934996820000106
Figure GDA0001934996820000107
wherein, the coefficient matrix corresponding to each formula is as follows:
Figure GDA0001934996820000111
Figure GDA0001934996820000112
Figure GDA0001934996820000113
Figure GDA0001934996820000114
Figure GDA0001934996820000115
Figure GDA0001934996820000116
Figure GDA0001934996820000117
Figure GDA0001934996820000118
Figure GDA0001934996820000119
Figure GDA00019349968200001110
Figure GDA00019349968200001111
thus, the formula (2) can be obtained from the formula (24) to the formula (28)2) Jacobian matrix of
Figure GDA00019349968200001112
The following were used:
Figure GDA00019349968200001113
wherein the content of the first and second substances,
Figure GDA00019349968200001114
k∈[T-N+1,T-1]；
the iteration termination condition in the algorithm flow is as follows:
||Fs(s)||≤σ (30)
wherein σ is the acceptable iteration precision error;
to this end, the root of the first order requirement (22) can be solved by the above steps
Figure GDA00019349968200001115
Then passing through the relation lambdak＝Q1wkAn optimal solution of the optimization problem (16) is obtained
Figure GDA0001934996820000121
Then passing through relation IIk＝JωkAn optimal solution of the optimization problem (15) can be obtained
Figure GDA0001934996820000122
When T is less than or equal to N, namely when the full information rolling time domain estimation is considered, the root of the first-order requirement (22) can be solved by only making all T-N in the algorithm steps equal to 0
Figure GDA0001934996820000123
Further by the relation lambdak＝Q1wkAn optimal solution of the optimization problem (16) is obtained
Figure GDA0001934996820000124
Then passing through relation IIk＝JωkCan be optimizedOptimal solution of problem (15)
Figure GDA0001934996820000125
5) At the time T, the optimal solution of the step 4) can be obtained
Figure GDA0001934996820000126
Or
Figure GDA0001934996820000127
And attitude control moment sequence
Figure GDA0001934996820000128
Or
Figure GDA0001934996820000129
Substituting the attitude dynamics and the kinematic equation of the discrete spacecraft into iteration and recursion to obtain the optimal estimated value of the angular velocity of the spacecraft at the current moment
Figure GDA00019349968200001210
6) Calculating the next time estimation error covariance P according to the formula (8)T-N+1；
7) At the next time T +1, the data sequence is updated with the latest measured value and the latest control torque, and T: ═ T +1 is defined, returning to step 2).
The optimal online estimation of the angular velocity of the spacecraft can be completed through the steps.
Compared with the prior art, the invention has the advantages that:
(1) the method describes the attitude dynamics and the kinematics equation of the discrete spacecraft by introducing a special orthogonal group SO (3), the kinematics equation uses multiplication of two rotation matrixes to describe the rotation motion of the spacecraft in a popular and understandable way and greatly simplify the calculation amount of the attitude motion of the spacecraft, and simultaneously compared with the traditional quaternion description method and the Euler angle description method, the constraint problem in the quaternion method and the singular point problem in the Euler angle method are avoided.
(2) Different from the conventional angular velocity estimation method, compared with the classical Kalman filtering, the estimation method of the angular velocity of the spacecraft based on the rolling time domain estimation provided by the invention has the advantages that on one hand, the state estimation is carried out by fully utilizing the system information and the measurement information in the current and past periods, more information is considered, and the estimation precision is further improved; on the other hand, angular velocity constraint and control moment constraint of the spacecraft system are considered and are directly expressed in an optimization problem, and the optimization is dynamically satisfied through online rolling optimization; and the optimization problem is solved in a limited time domain by adopting the approximate arrival cost function, so that the solving time of the optimization problem is shortened. The rolling time domain estimation method is an estimation method which can solve the optimization problem including constraint conditions, and is high in reliability, accurate and effective.
(3) In the optimal problem solving algorithm, the optimization problem in a limited time domain is solved by adopting a multiple shooting method, compared with other optimal algorithms, the algorithm can process a rigid problem and can adopt a parallel framework, and the method has the characteristics of high calculation efficiency and good operation effect, so that the solving efficiency and precision are improved, and the real-time requirement of a spacecraft system is met.
Drawings
Fig. 1 is a schematic block diagram of the method for estimating the angular velocity of a spacecraft based on rolling time domain estimation according to the present invention.
Fig. 2 is a schematic diagram of the rolling time domain estimation method of the present invention.
Fig. 3 is a flow chart of the method for estimating the angular velocity of a spacecraft based on rolling time domain estimation according to the present invention.
FIG. 4 is a flow chart of the multiple targeting algorithm of the present invention to solve the optimization problem.
Detailed Description
As shown in fig. 1, the spacecraft angular velocity estimation method of the present invention is a rolling time domain estimation method considering the spacecraft angular velocity physical constraint and the control moment constraint. At each sampling moment, the optimal angular velocity estimation value obtained by the estimation method can be fed back to the controller, and the controller generates a control moment to act on the spacecraft interfered by the external moment by using the optimal estimation value to complete attitude stabilization or attitude maneuver; the redundant fiber optic gyroscope is combined with the angular velocity of the sensitive spacecraft, new weighted data fusion processing is carried out on four paths of twelve-dimensional measurement signals with measurement noise on the basis of fault failure factor distribution, and the new weighted data fusion processing and the control moment signals are input into a rolling time domain estimator, so that the real-time, reliable, accurate and efficient estimation of the triaxial angular velocity of the spacecraft is completed; and pushing the latest measured value and the latest control moment value into the data sequence at the next sampling moment to finish updating the data sequence, thereby forming a closed loop structure.
As shown in fig. 2, the rolling time domain estimation algorithm of the present invention is derived from the analysis of a rolling time domain window, each time only uses the latest first N data at the current time T, and the influence of the rest of the measurement data on the estimation can be described by an arrival cost function; at the next moment, the rolling time domain window is shifted backwards by one sampling time, and the data sequence is updated to the latest first N data at the current T +1 moment.
It should be understood that in this context, spacecraft angular velocities describe "vectors" and spacecraft angular velocities describe "scalars", i.e. angular velocities refer to (ω)x,ωy,ωz) Angular rate refers to the angular rate ω of a certain principal axis of inertiaxOr ωyOr ωzAnd the three-axis angular rate refers to ωx,ωy,ωzThree scalars.
As shown in fig. 3, the method for estimating the angular velocity of a spacecraft based on rolling time domain estimation of the present invention includes the steps of: firstly, constructing a state equation and a measurement equation which can be used for a rolling time domain estimation strategy based on a discrete spacecraft attitude dynamics and kinematics model; then, carrying out new weighted data fusion of weight values distributed based on fault failure factors on four paths of twelve-dimensional angular rate signals sensed by the combination of the redundant fiber optic gyroscopes to obtain angular rate measured values of all inertial main shafts under a spacecraft body coordinate system; and finally, according to the state equation, the measurement equation, the angular velocity measurement sequence and the controller input sequence, simultaneously considering the angular velocity constraint and the control moment constraint of the spacecraft, and performing online estimation by using a rolling time domain estimation strategy and a multiple target shooting method to obtain the optimal estimation of the angular velocity of the spacecraft. The specific operation steps are as follows:
firstly, aiming at the problem of spacecraft angular velocity estimation, a state equation and an observation equation of the spacecraft angular velocity are constructed based on a spacecraft dynamics model in a special orthogonal group SO (3) description mode. The spacecraft dynamics model under the special orthogonal group SO (3) description mode is as follows:
Figure GDA0001934996820000141
Rk+1＝RkFk(2)
Figure GDA00019349968200001412
wherein h is a sampling time interval; superscript T represents matrix transposition; the index k denotes the kth sampling instant; matrix Rk,FkE SO (3) is a rotation matrix,
Figure GDA0001934996820000142
a rotation matrix representing the euler attitude angle description of the spacecraft,
Figure GDA0001934996820000143
a rotation matrix representing the description of the Euler attitude error angle of the spacecraft between two sampling instants, FkIs the solution of implicit equation (1);
Figure GDA0001934996820000144
denotes the overall angular momentum of the spacecraft, J ═ diag ([2.5, 2.5)])kg·m2A matrix of the inertia of the spacecraft is represented,
Figure GDA0001934996820000145
representing the three-axis angular velocity of the spacecraft;
Figure GDA0001934996820000146
the method comprises the following steps of representing external disturbance moment caused by gravity gradient moment, aerodynamic resistance, sunlight pressure and the like to which a spacecraft is subjected;
Figure GDA0001934996820000147
representing the control moment required by the attitude maneuver or the stability of the spacecraft; positive definite array
Figure GDA0001934996820000148
Represents a non-standard inertia matrix, which is related to the inertia matrix J;
Figure GDA0001934996820000149
is about real vector
Figure GDA00019349968200001410
For any vector
Figure GDA00019349968200001411
Its corresponding oblique antisymmetric array x×The form is as follows:
Figure GDA0001934996820000151
based on the discrete spacecraft attitude dynamics model, a state equation and an observation equation for rolling time domain estimation can be constructed as follows:
taking the system state quantity as omega ═ omegax,ωy,ωz)TThen the state equation is
Figure GDA0001934996820000152
Wherein ω ═ ω (ω ═ ω)x,ωy,ωz)TIs the spacecraft angular velocity; fkThe epsilon SO (3) is an auxiliary variable and represents the relative attitude between two adjacent moments;
Figure GDA0001934996820000153
is system external interference, independent of measurement noise and initial state;
taking the observed quantity of the system as y ═ yx,yy,yz)TThen the observation equation is
Jyk＝CJωk+νk(5)
Wherein y ═ yx,yy,yz)TMeasuring data output value, y, for redundant fiber optic gyroscopes combined data processed angular velocityx,yy,yzMeasuring the angular rate of each inertia main shaft of the spacecraft;
Figure GDA0001934996820000154
for noise measurement, it is independent of system external interference and initial state; measurement matrix C ═ I3×3Is an identity matrix;
secondly, because the observation equation needs three-dimensional measurement data, and the combination of the redundant fiber optic gyroscopes feels that the combination of the redundant fiber optic gyroscopes is four-way twelve-dimensional measurement data, a new weighted data fusion method for distributing weights based on failure factors needs to be designed:
Figure GDA0001934996820000155
wherein the content of the first and second substances,
Figure GDA0001934996820000156
for the spacecraft angular velocity measurement data after data processing, it should be noted that y here1,y2,y3And y mentioned abovex,yy,yzBeing the same quantity, the above observation equation requires three-dimensional angular rate measurements, whereas redundant fiber optic gyroscopes are sensitive to four-way twelve-dimensional angular rate measurements, so that from twelve-dimensional to three-dimensional data processing is required for use with the above observation equation.
Figure GDA0001934996820000161
The angular velocity information of the spacecraft sensed by the ith optical fiber gyroscope in the redundant optical fiber gyroscopes is measured by four paths of twelve-dimensional angular velocity measurement information
Figure GDA0001934996820000162
α1＝45°,α2＝135°,α3＝225°,α4315 degrees are respectively the installation deflection angles of four optical fiber gyroscopes, and represent the included angle between the projection of the axial direction of the ith optical fiber gyroscope on the xoy plane and the ox axis, βi54.75 degrees, and i is 1,2,3 and 4 respectively representing the installation high and low angles of the four optical fiber gyroscopes, and representing the included angle between the axial direction of the ith optical fiber gyroscope and the xoy plane;
Figure GDA0001934996820000163
as a weighting factor, influenced by the failure factor of the ith optical fiber gyroscope, where ρ is takeni1, i is 1,2,3,4, so the weighting factor is τi＝0.25,i＝1,2,3,4τi＝0.25,i＝1,2,3,4；
And thirdly, performing optimal estimation on the angular velocity of the spacecraft by using a rolling time domain estimation algorithm on the basis of the first step and the second step, wherein the method specifically comprises the following steps:
1) a priori estimation of initial state assuming a known system
Figure GDA0001934996820000164
Covariance matrix of initial state error is P0＝diag([1,1,1]) The weight matrix of the external interference and the measured noise of the system is Q1＝diag([100,100,100]),Q2＝diag([1000,1000,1000]) And P is0,Q1,Q2All the time zones are positive definite matrixes, and the length of a rolling time domain is N-5;
2) at time T, the optimization problem defined by the rolling time domain estimation is as follows:
Figure GDA0001934996820000165
wherein the content of the first and second substances,
Figure GDA0001934996820000166
in order to achieve the cost function, the invention adopts the arrival cost function of an unconstrained linear system to approximate the cost function of a constrained nonlinear system, namely the approximate arrival cost function has the following form:
Figure GDA0001934996820000167
in the formula, initial state estimation
Figure GDA0001934996820000168
Obey normal distribution
Figure GDA0001934996820000169
PT-N is the estimation error covariance, updated by equation (8) below;
Figure GDA0001934996820000171
meanwhile, the objective function needs to satisfy spacecraft dynamics, spacecraft angular velocity time domain constraint and control moment constraint, and the specific mathematical expression form is as follows:
angular velocity time domain constraint:
Figure GDA0001934996820000172
wherein the content of the first and second substances,
Figure GDA0001934996820000173
for the angular rate estimate obtained by the rolling time domain estimation,
Figure GDA0001934996820000174
the minimum and maximum angular velocities allowed by the main axes of the spacecraft are determined by the properties of physical devices such as solar sailboards, sensors and the like installed on the spacecraft. It should be noted that the angular velocity time domain constraint proposed herein is the optimal angular velocity estimated by the rolling time domain online estimation algorithm
Figure GDA0001934996820000175
Constraining rather than imposing an actual operating angular velocity ω of the spacecraftkTo make the constraints, special care needs to be taken to distinguish in order to avoid ambiguity.
Controlling torque constraint:
Figure GDA0001934996820000176
and because the spacecraft system is in a discrete system form, the angular acceleration
Figure GDA0001934996820000177
Can be approximated in the following form:
Figure GDA0001934996820000178
the control torque constraint can be written approximately as:
Figure GDA0001934996820000179
simplifying the above formula, can convert control moment restraint into the form of angular velocity restraint, the concrete form is as follows:
Figure GDA00019349968200001710
wherein the content of the first and second substances,
Figure GDA00019349968200001711
respectively providing minimum and maximum control moments for the spacecraft actuating mechanism in the directions of all main shafts; t iss0.01s is a discrete system sampling time interval;
Figure GDA00019349968200001712
the optimal estimated value of the angular velocity of the spacecraft at the previous moment is obtained. Similarly, it should be noted that the constraint equation (13) is the optimal angular velocity estimated by the rolling time domain online estimation algorithm
Figure GDA0001934996820000181
Constraining rather than imposing an actual operating angular velocity ω of the spacecraftkAnd (6) carrying out constraint.
Combining the angular velocity constraint (9) and the control moment constraint (13), two constraints can be simplified as:
Figure GDA0001934996820000182
Figure GDA0001934996820000183
Figure GDA0001934996820000184
at this time, the constrained rolling horizon optimization problem can be written as follows:
Figure GDA0001934996820000185
Figure GDA0001934996820000186
Figure GDA0001934996820000187
Figure GDA0001934996820000188
3) at the time T, if the redundant fiber optic gyroscope after data processing is combined with the angular speed measurement output sequence
Figure GDA0001934996820000189
Or
Figure GDA00019349968200001810
And spacecraft attitude control moment sequence
Figure GDA00019349968200001811
Or
Figure GDA00019349968200001812
If known, the optimization problem (15) can be solved by using a multiple targeting method, and the only existing optimal solution can be solved
Figure GDA00019349968200001813
4) At time T, the specific steps for solving the optimization problem (15) by using the multiple targeting method are as follows:
firstly, the angular momentum pi of the spacecraft is used herek＝JωkReplacing spacecraft angular velocity omegakThe optimization problem (15) can be equivalently written in the form of the following optimization problem (16):
Figure GDA00019349968200001814
Figure GDA00019349968200001815
Figure GDA0001934996820000191
Figure GDA0001934996820000192
wherein, Yk＝JykFrom the spacecraft angular velocity measurement y for time kkThe angular momentum of the spacecraft is obtained by calculation,
Figure GDA0001934996820000193
for the estimated value of the angular momentum of each main shaft of the spacecraft at the moment k,
Figure GDA0001934996820000194
Figure GDA0001934996820000195
respectively the minimum and maximum angular momentum of each main shaft of the spacecraft, and the formula (16c) is used for solving auxiliary variables
Figure GDA0001934996820000196
Implicit equation of (2).
Considering system dynamics equality constraint (16b) and angular momentum time domain inequality constraint (16d), defining an augmented objective function:
Figure GDA0001934996820000197
wherein the content of the first and second substances,
Figure GDA0001934996820000198
is a Lagrangian multiplier,. mu. > 0 and
Figure GDA0001934996820000199
in order to be a penalty parameter,
Figure GDA00019349968200001910
the penalty function for the corresponding inequality constraint (16d) is expressed in the following way:
Figure GDA00019349968200001911
wherein h is 0.01, which is the time interval between two moments;
when T is larger than N, the broadening cost function is subjected to variation by using a variation method, wherein the variation expression is as follows:
Figure GDA00019349968200001912
wherein the content of the first and second substances,
Figure GDA00019349968200001913
as a penalty function
Figure GDA00019349968200001914
To pair
Figure GDA00019349968200001915
First order variational sum, variational coefficient
Figure GDA0001934996820000201
As follows:
Figure GDA0001934996820000202
let δ ΦTaThe first order requirement for the 0 available optimization problem (16) is as follows:
Figure GDA0001934996820000203
wherein the content of the first and second substances,
Figure GDA0001934996820000204
the following first order requirements can be met by modifying equation (21):
Figure GDA0001934996820000205
wherein the content of the first and second substances,
Figure GDA0001934996820000206
for the vectors to be solved of the first order requirement, by lagrange multiplier sequence
Figure GDA0001934996820000207
And spacecraft angular momentum estimation sequence
Figure GDA0001934996820000208
Composition is carried out;
at this point, solving the optimization problem (16) can be equivalently transformed into solving the root of the first order requirement (22), i.e., the root of the first order requirement is the same as the optimal solution of the optimization problem. The numerical iterative solution of equation (22) using multiple targeting is to find a specific set of sequences s*So that Fs(s*)＝06N+3If the first order requirement (22) has a root and is unique, the root is the optimal solution of the optimization problem (16);
the flow of the multiple targeting algorithm is shown in fig. 4, and the root-finding iterative relationship in the algorithm is as follows:
si+1＝si+κΔ(s)i(23)
wherein s isi+1For the i +1 th iterative solution, siFor the ith iterative solution,. kappa.. epsilon. (0, 1)]In order to iterate the step size,
Figure GDA0001934996820000211
for iterative descending direction, Δ(s)iIn (1)
Figure GDA0001934996820000212
A Jacobian matrix of formula (22);
a variational expression according to equation (22):
Figure GDA0001934996820000213
Figure GDA0001934996820000214
Figure GDA0001934996820000215
Figure GDA0001934996820000216
Figure GDA0001934996820000217
wherein, the coefficient matrix corresponding to each formula is as follows:
Figure GDA0001934996820000221
Figure GDA0001934996820000222
Figure GDA0001934996820000223
Figure GDA0001934996820000224
Figure GDA0001934996820000225
Figure GDA0001934996820000226
Figure GDA0001934996820000227
Figure GDA0001934996820000228
Figure GDA0001934996820000229
Figure GDA00019349968200002210
Figure GDA00019349968200002211
therefore, the Jacobian matrix of formula (22) can be obtained from formula (24) to formula (28)
Figure GDA00019349968200002212
The following;
Figure GDA00019349968200002213
wherein the content of the first and second substances,
Figure GDA00019349968200002214
k∈[T-N+1,T-1]；
the iteration termination condition in the algorithm flow is as follows:
||Fs(s)||≤σ (30)
wherein, σ is 0.00001 as iteration precision error;
to this end, the root of the first order requirement (22) can be solved by the above steps
Figure GDA00019349968200002215
Then passing through the relation lambdak＝Q1wkAn optimal solution of the optimization problem (16) is obtained
Figure GDA0001934996820000231
Then passing through relation IIk＝JωkAn optimal solution of the optimization problem (15) can be obtained
Figure GDA0001934996820000232
When T is less than or equal to N, namely when the full information rolling time domain estimation is considered, the root of the first-order requirement (22) can be solved by only making all T-N in the algorithm steps equal to 0
Figure GDA0001934996820000233
Further by the relation lambdak＝Q1wkAn optimal solution of the optimization problem (16) is obtained
Figure GDA0001934996820000234
Then passing through relation IIk＝JωkAn optimal solution of the optimization problem (15) can be obtained
Figure GDA0001934996820000235
5) At time T, the optimal solution of the steps can be obtained
Figure GDA0001934996820000236
Or
Figure GDA0001934996820000237
And attitude control moment sequenceColumn(s) of
Figure GDA0001934996820000238
Or
Figure GDA0001934996820000239
Substituting the attitude dynamics and the kinematic equation of the discrete spacecraft into iteration and recursion to obtain the optimal estimated value of the angular velocity of the spacecraft at the current moment
Figure GDA00019349968200002310
6) Calculating the next time estimation error covariance P according to the formula (8)T-N+1；
7) At the next time T +1, updating the data sequence by using the latest measured value and the latest control torque, defining T: (T + 1), and returning to the step 2);
the optimal online estimation of the angular velocity of the spacecraft can be completed through the steps.
Those skilled in the art will appreciate that the invention may be practiced without these specific details.

Claims (1)

1. A spacecraft angular velocity estimation method based on rolling time domain estimation is characterized by comprising the following steps:
s1: constructing a state equation and an observation equation of the angular velocity of the spacecraft based on a spacecraft dynamics model in a special orthogonal group SO (3) description mode;
the spacecraft dynamics model under the special orthogonal group SO (3) description mode is as follows:
Figure FDA0002388131010000011
Rk+1＝RkFk(2)
Figure FDA0002388131010000012
wherein h is a sampling time interval; superscript T represents matrix transposition;the index k denotes the kth sampling instant; matrix Rk,FkE SO (3) is a rotation matrix,
Figure FDA0002388131010000013
a rotation matrix representing the euler attitude angle description of the spacecraft,
Figure FDA0002388131010000014
a rotation matrix representing the description of the Euler attitude error angle of the spacecraft between two sampling instants, FkIs the solution of implicit equation (1);
Figure FDA0002388131010000015
which represents the angular momentum of the spacecraft,
Figure FDA0002388131010000016
Figure FDA0002388131010000017
is an inertial matrix of the spacecraft,
Figure FDA0002388131010000018
the angular velocity of the spacecraft at the kth sampling moment;
Figure FDA0002388131010000019
Figure FDA00023881310100000110
representing external disturbance moments caused by external disturbances including gravity gradient moments, aerodynamic drag, sunlight pressure, to which the spacecraft is subjected;
Figure FDA00023881310100000111
representing the control moment required by the attitude maneuver or the stability of the spacecraft; positive definite array
Figure FDA00023881310100000112
Representing non-standard inertia matricesIt is related to the inertia matrix J;
Figure FDA00023881310100000113
is about real vector
Figure FDA00023881310100000114
For any vector
Figure FDA00023881310100000115
Its corresponding oblique antisymmetric array x×The form is as follows:
Figure FDA00023881310100000116
the state equation and the observation equation for rolling time domain estimation are constructed based on a spacecraft dynamics model as follows:
taking the system state quantity as spacecraft angular velocity omega ═ omegax,ωy,ωz)TThen the state equation is
Figure FDA0002388131010000021
Taking the system observed quantity as the angular velocity measurement data output value y after the redundant fiber-optic gyroscope combined data is processed (y ═ y)x,yy,yz)TThen the observation equation is
Jyk＝CJωk+νk(5)
Wherein, yx,yy,yzMeasuring the angular rate of each inertia main shaft of the spacecraft;
Figure FDA0002388131010000022
and
Figure FDA0002388131010000023
respectively system external interference and measurement noise;
Figure FDA0002388131010000024
is a measurement matrix;
s2: simplifying four-way twelve-dimensional angular rate measurement information combined by redundant fiber optic gyroscopes by using a new weighted data fusion method for distributing weights based on failure factors to obtain spacecraft three-axis angular rate measurement data;
the new weighted data fusion method for distributing the weight value based on the failure factor comprises the following steps:
Figure FDA0002388131010000025
wherein the content of the first and second substances,
Figure FDA0002388131010000026
the angular velocity information of the spacecraft sensed by the ith optical fiber gyroscope in the redundant optical fiber gyroscopes,
Figure FDA0002388131010000027
for four-way twelve-dimensional angular rate measurement information αiWherein i is 1,2,3 and 4 is the installation deflection angle of the optical fiber gyroscope and represents the included angle between the projection of the axial direction of the ith optical fiber gyroscope on the xoy plane and the ox axis, βiThe i is 1,2,3 and 4, which is the installation high-low angle of the optical fiber gyroscope and represents the included angle between the axial direction of the ith optical fiber gyroscope and the xoy plane;
Figure FDA0002388131010000028
as a weighting factor, the failure factor ρ of the ith optical fiber gyroscopei∈[0,1]Influence and need to satisfy the constraint τ1+τ2+τ3+τ4＝1；
S3: based on the state equation and the observation equation constructed in the step S1 and the spacecraft triaxial angular rate measurement data obtained in the step S2, considering spacecraft angular velocity constraint and control moment constraint, and performing online estimation by using a rolling time domain estimation strategy and a multiple target shooting method to obtain an estimated value of the optimal angular velocity of the spacecraft at the current moment;
the specific steps of on-line estimation by using a rolling time domain estimation strategy and a multiple target practice method are as follows:
1) a priori estimation of initial state assuming a known system
Figure FDA0002388131010000031
Covariance matrix of initial state error is P0The weight matrix of the external interference and the measured noise of the system is Q1,Q2And P is0,Q1,Q2All are positive definite matrixes, and the length of a rolling time domain is N;
2) at time T, the optimization problem defined by the rolling time domain estimation is as follows:
Figure FDA0002388131010000032
in the two cases, the first term in the objective function represents the cost sum of external interference of the system at each sampling moment, the second term represents the cost sum of measurement noise at each sampling moment, and the third term represents prior estimation when T is less than or equal to N
Figure FDA0002388131010000033
Degree of trust of, at T>When the N is greater than the N value,
Figure FDA0002388131010000034
to approximate the arrival cost function, which represents the influence of the remaining measurement data on the estimation, the arrival cost function of the unconstrained linear system is used to approximate the cost function of the constrained nonlinear system, i.e., the approximate arrival cost function is of the form:
Figure FDA0002388131010000035
in the formula, initial state estimation
Figure FDA0002388131010000036
Obey normal distribution
Figure FDA0002388131010000037
PT-NTo estimate the error covariance, it can be updated by the following equation (8);
Figure FDA0002388131010000038
meanwhile, the objective function needs to satisfy spacecraft dynamics, spacecraft angular velocity time domain constraint and control moment constraint, and the specific mathematical expression form is as follows:
angular velocity time domain constraint:
Figure FDA0002388131010000039
wherein the content of the first and second substances,
Figure FDA0002388131010000041
for the angular rate estimate obtained by the rolling time domain estimation,
Figure FDA0002388131010000042
the minimum and maximum angular rates allowed by each main shaft of the spacecraft are determined by the properties of physical devices including a solar panel and a sensor which are arranged on the spacecraft;
controlling torque constraint:
Figure FDA0002388131010000043
and because the spacecraft system is in a discrete system form, the angular acceleration
Figure FDA0002388131010000044
Can be approximated in the following form:
Figure FDA0002388131010000045
the control torque constraint can be written approximately as:
Figure FDA0002388131010000046
simplifying the above formula, can convert control moment restraint into the form of angular velocity restraint, the concrete form is as follows:
Figure FDA0002388131010000047
wherein the content of the first and second substances,
Figure FDA0002388131010000048
respectively providing minimum and maximum control moments for the spacecraft actuating mechanism in the directions of all main shafts;
Figure FDA0002388131010000049
the optimal estimated value of the angular rate of the spacecraft at the previous moment is obtained; t issSampling a time interval for a discrete system;
combining the angular velocity constraint (9) and the control moment constraint (13), two constraints can be simplified as:
Figure FDA00023881310100000410
Figure FDA00023881310100000411
Figure FDA00023881310100000412
at this time, the constrained rolling horizon optimization problem can be written as follows:
Figure FDA00023881310100000413
Figure FDA0002388131010000051
Figure FDA0002388131010000052
Figure FDA0002388131010000053
3) at the time T, if the redundant fiber optic gyroscope after data processing is combined with the angular speed measurement output sequence
Figure FDA0002388131010000054
Or
Figure FDA0002388131010000055
And spacecraft attitude control moment sequence
Figure FDA0002388131010000056
Or
Figure FDA0002388131010000057
Known, therein for arbitrary variables
Figure FDA0002388131010000058
Sequence of
Figure FDA0002388131010000059
To represent
Figure FDA00023881310100000510
The optimization problem (15) can be solved using a multiple targeting method to obtain a unique existing optimal solution
Figure FDA00023881310100000511
4) At time T, the specific steps for solving the optimization problem (15) by using the multiple targeting method are as follows:
firstly, the angular momentum pi of the spacecraft is used herek＝JωkReplacing spacecraft angular velocity omegakThe optimization problem (15) can be equivalently transformed into the form of the following optimization problem (16):
Figure FDA00023881310100000512
Figure FDA00023881310100000513
Figure FDA00023881310100000514
Figure FDA00023881310100000515
wherein, Yk＝JykFrom the spacecraft angular velocity measurement y for time kkCalculating the angular momentum of the spacecraft;
Figure FDA00023881310100000516
3, estimating values of angular momentum of each main shaft of the spacecraft at the moment k;
Figure FDA00023881310100000517
Figure FDA00023881310100000518
respectively the minimum and maximum angular momentum of each main shaft of the spacecraft; equation (16c) is solving for the auxiliary variables
Figure FDA00023881310100000519
The implicit equation of (a) is,
considering system dynamics equality constraint (16b) and angular momentum time domain inequality constraint (16d), defining an augmented objective function:
Figure FDA0002388131010000061
wherein the content of the first and second substances,
Figure FDA0002388131010000062
is a Lagrangian multiplier,. mu. > 0 and
Figure FDA0002388131010000063
in order to be a penalty parameter,
Figure FDA0002388131010000064
the penalty function for the corresponding inequality constraint (16d) is expressed in the following way:
Figure FDA0002388131010000065
when T is larger than N, the broadening cost function is subjected to variation by using a variation method, wherein the variation expression is as follows:
Figure FDA0002388131010000066
wherein the content of the first and second substances,
Figure FDA0002388131010000067
as a penalty function
Figure FDA0002388131010000068
To pair
Figure FDA0002388131010000069
First order variational sum, variational coefficient
Figure FDA00023881310100000610
As follows:
Figure FDA00023881310100000611
wherein e isiI is 1,2 and 3 are unit vectors of main axes of inertia of the spacecraft respectively;
let δ ΦTaThe first order requirement for the 0 available optimization problem (16) is as follows:
Figure FDA0002388131010000071
wherein the content of the first and second substances,
Figure FDA0002388131010000072
the following first order requirements can be met by modifying equation (21):
Figure FDA0002388131010000073
wherein the content of the first and second substances,
Figure FDA0002388131010000074
for the vectors to be solved of the first order requirement, by lagrange multiplier sequence
Figure FDA0002388131010000075
And spacecraft angular momentum estimation sequence
Figure FDA0002388131010000076
Composition is carried out;
at this time, solving the optimization problem (16) can be equivalently converted into solving the root of the first-order requirement (22), namely, the root of the first-order requirement and the optimal solution of the optimization problem are the same, and the numerical iteration solving is carried out on the formula (22) by adopting a multiple targeting method, namely, a group of specific sequences s is required to be found*So that Fs(s*)＝06N+3If the first order requirement (22) is rooted and unique, then the root is the optimization problem (16)) The optimal solution of (2);
the root-finding iterative relationship in the multiple targeting method is as follows:
si+1＝si+κΔ(s)i(23)
wherein s isi+1For the i +1 th iterative solution, siFor the ith iterative solution,. kappa.. epsilon. (0, 1)]In order to iterate the step size,
Figure FDA0002388131010000081
for iterative descending direction, Δ(s)iIn (1)
Figure FDA0002388131010000082
A Jacobian matrix of formula (22);
a variational expression according to equation (22):
Figure FDA0002388131010000083
Figure FDA0002388131010000084
Figure FDA0002388131010000085
Figure FDA0002388131010000086
Figure FDA0002388131010000087
wherein, the coefficient matrix corresponding to each formula is as follows:
Figure FDA0002388131010000088
Figure FDA0002388131010000089
Figure FDA00023881310100000810
Figure FDA00023881310100000811
Figure FDA00023881310100000812
Figure FDA00023881310100000813
Figure FDA00023881310100000814
Figure FDA00023881310100000815
Figure FDA00023881310100000816
Figure FDA00023881310100000817
Figure FDA00023881310100000818
therefore, the Jacobian matrix of formula (22) can be obtained from formula (24) to formula (28)
Figure FDA0002388131010000091
The following were used:
Figure FDA0002388131010000092
wherein the content of the first and second substances,
Figure FDA0002388131010000093
the iteration termination condition in the algorithm flow is as follows:
||Fs(s)||≤σ (30)
wherein σ is the acceptable iteration precision error;
to this end, the root of the first order requirement (22) can be solved by the above steps
Figure FDA0002388131010000094
Then passing through the relation lambdak＝Q1wkAn optimal solution of the optimization problem (16) is obtained
Figure FDA0002388131010000095
Then passing through relation IIk＝JωkAn optimal solution of the optimization problem (15) can be obtained
Figure FDA0002388131010000096
When T is less than or equal to N, namely when the full information rolling time domain estimation is considered, the root of the first-order requirement (22) can be solved by only making all T-N in the algorithm steps equal to 0
Figure FDA0002388131010000097
Further by the relation lambdak＝Q1wkAn optimal solution of the optimization problem (16) is obtained
Figure FDA0002388131010000098
Then passing through relation IIk＝JωkAn optimal solution of the optimization problem (15) can be obtained
Figure FDA0002388131010000099
5) At time T, canOptimal solution of step 4)
Figure FDA00023881310100000910
Or
Figure FDA00023881310100000911
And attitude control moment sequence
Figure FDA00023881310100000912
Or
Figure FDA00023881310100000913
Substituting the attitude dynamics and the kinematic equation of the discrete spacecraft into iteration and recursion to obtain the optimal estimated value of the angular velocity of the spacecraft at the current moment
Figure FDA0002388131010000101
6) Calculating the next time estimation error covariance P according to the formula (8)T-N+1；
7) At the next time T +1, the data sequence is updated with the latest measured value and the latest control torque, and T: ═ T +1 is defined, and the process returns to step 2.
CN201811193664.9A 2018-10-15 2018-10-15 Spacecraft angular velocity estimation method based on rolling time domain estimation Active CN109343550B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811193664.9A CN109343550B (en) 2018-10-15 2018-10-15 Spacecraft angular velocity estimation method based on rolling time domain estimation

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811193664.9A CN109343550B (en) 2018-10-15 2018-10-15 Spacecraft angular velocity estimation method based on rolling time domain estimation

Publications (2)

Publication Number Publication Date
CN109343550A CN109343550A (en) 2019-02-15
CN109343550B true CN109343550B (en) 2020-04-21

Family

ID=65309998

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811193664.9A Active CN109343550B (en) 2018-10-15 2018-10-15 Spacecraft angular velocity estimation method based on rolling time domain estimation

Country Status (1)

Country Link
CN (1) CN109343550B (en)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109813312A (en) * 2019-03-05 2019-05-28 南京理工大学 The measuring device and method at a kind of card form measurement gestures of object angle
CN110006425A (en) * 2019-04-11 2019-07-12 南京航空航天大学 High dynamic Attitude rate estimator method based on carrier kinetic model auxiliary
CN110119153B (en) * 2019-05-10 2020-12-15 北京航空航天大学 Under-actuated spacecraft attitude control method under active assistance of light pressure moment
US11124196B2 (en) * 2019-07-02 2021-09-21 Mitsubishi Electric Research Laboratories, Inc. Receding horizon state estimator
CN110736468A (en) * 2019-09-23 2020-01-31 中国矿业大学 Spacecraft attitude estimation method assisted by self-adaptive kinematics model
CN111578931B (en) * 2020-05-21 2022-03-04 中国人民解放军海军航空大学 High-dynamic aircraft autonomous attitude estimation method based on online rolling time domain estimation
CN114167493B (en) * 2021-11-23 2023-08-04 武汉大学 Seismic rotation measurement system and method of GNSS double-antenna auxiliary gyroscope
CN114756047B (en) * 2022-06-14 2022-09-02 东方空间技术(山东)有限公司 Method and device for controlling movement of spacecraft

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP3555003B2 (en) * 1997-08-29 2004-08-18 株式会社日立製作所 Mobile device position measurement device
EP2908155A1 (en) * 2014-02-14 2015-08-19 Thales Method for anti-jamming processing of a radio signal, related module and computer program

Family Cites Families (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2001080597A (en) * 1999-09-13 2001-03-27 Mitsubishi Electric Corp Attitude control device for three-axis stability satellite
US6254036B1 (en) * 2000-02-16 2001-07-03 Hughes Electronics Corporation Method and apparatus for spacecraft wheel desaturation
CN101706512B (en) * 2009-11-25 2011-06-15 哈尔滨工业大学 Method for estimating pseudo rate of spacecraft based on attitude measurement information of star sensors and angular momentum measurement information of flywheels
CN101739655A (en) * 2009-12-17 2010-06-16 浙江工业大学 Method for scheduling public slow system dynamically based on rolling horizon scheduling algorithm
CN104063749B (en) * 2014-06-28 2017-03-29 中国人民解放军国防科学技术大学 A kind of autonomous mission planning algorithm of the imaging satellite based on roll stablized loop
US9646361B2 (en) * 2014-12-10 2017-05-09 Siemens Healthcare Gmbh Initialization independent approaches towards registration of 3D models with 2D projections
CN104406598B (en) * 2014-12-11 2017-06-30 南京航空航天大学 A kind of non-cooperative Spacecraft Attitude estimation method based on virtual sliding formwork control
CN105629732B (en) * 2016-01-29 2018-03-09 北京航空航天大学 A kind of spacecraft attitude output Tracking Feedback Control method for considering Control constraints
CN107483541A (en) * 2017-07-17 2017-12-15 广东工业大学 A kind of online task immigration method based on rolling time horizon
CN108153322B (en) * 2017-12-06 2019-03-29 北京航空航天大学 A kind of spacecraft attitude tracking adaptive fault tolerant control method for the rotary inertia considering time-varying
CN107942090B (en) * 2017-12-28 2019-10-29 北京航空航天大学 A kind of spacecraft Attitude rate estimator method for extracting Optic flow information based on fuzzy star chart

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP3555003B2 (en) * 1997-08-29 2004-08-18 株式会社日立製作所 Mobile device position measurement device
EP2908155A1 (en) * 2014-02-14 2015-08-19 Thales Method for anti-jamming processing of a radio signal, related module and computer program

Also Published As

Publication number Publication date
CN109343550A (en) 2019-02-15

Similar Documents

Publication Publication Date Title
CN109343550B (en) Spacecraft angular velocity estimation method based on rolling time domain estimation
CN110030994B (en) Monocular-based robust visual inertia tight coupling positioning method
CN111156987B (en) Inertia/astronomy combined navigation method based on residual compensation multi-rate CKF
CN108663936B (en) Model does not know spacecraft without unwinding Attitude Tracking finite-time control method
CN111189442B (en) CEPF-based unmanned aerial vehicle multi-source navigation information state prediction method
Bar-Itzhack et al. Recursive attitude determination from vector observations Euler angle estimation
Christian et al. Sequential optimal attitude recursion filter
CN110146092B (en) Double-body asteroid detection track optimization method based on navigation information evaluation
CN115046545A (en) Positioning method combining deep network and filtering
CN111578931B (en) High-dynamic aircraft autonomous attitude estimation method based on online rolling time domain estimation
CN113029173A (en) Vehicle navigation method and device
CN105136150A (en) Attitude determination method based on multiple star-sensor measure information fusion
Guo et al. Analysis and design of an attitude calculation algorithm based on elman neural network for SINS
Zhe et al. Adaptive complementary filtering algorithm for imu based on mems
Liu et al. A novel hybrid attitude fusion method based on LSTM neural network for unmanned aerial vehicle
CN113377006B (en) Global fast terminal sliding mode control method based on invariant flow observer
CN113587926B (en) Spacecraft space autonomous rendezvous and docking relative navigation method
CN105180946A (en) Wideband measurement-based satellite high-precision attitude determination method and system thereof
CN113359444B (en) Flexible spacecraft rigid-flexible coupling characteristic intelligent identification method based on neural network
CN110209190B (en) Satellite nominal orbit unbiased flight control method
Chen et al. High precision attitude estimation algorithm using three star trackers
Lu et al. Analysis and Evaluation of MEMS-IMU Attitude Estimation Algorithm
Li et al. Event-Triggered Moving Horizon Pose Estimation for Spacecraft Systems
CN116817932B (en) Integrated method for track maintenance, track determination and real-time mapping of gravitational field
CN116358564B (en) Unmanned aerial vehicle bee colony centroid motion state tracking method, system, equipment and medium

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