An energy-optimized high-precision magnetron and thruster combined attitude control method

Abstract

The application discloses a high-precision magnetic control and thruster combined attitude control method with optimal energy consumption, which comprises the following steps: firstly, calculating an expected control moment according to an attitude control amount and an attitude angular velocity control amount by using a PD control law, and combining a gravity gradient moment estimation to perform moment compensation; simultaneously, filtering and estimating a satellite residual magnetic moment according to angular velocity measurement information, combining a gravity gradient moment estimation and a previous period theoretical magnetic control moment estimation; then, distributing satellite three-axis control according to a satellite magnetic control axis / thruster control axis setting, performing jet phase plane protection calculation on the direction of thruster control, and outputting jet pulse width when necessary; distributing and calculating a small disturbance magnetic moment on the direction of magnetic control according to geomagnetic field information and the expected control moment, and combining the residual magnetic moment estimation result to perform magnetic moment compensation and three-axis magnetic moment output; finally, performing theoretical magnetic control moment estimation according to the three-axis magnetic moment calculation result, which is used for calculation in the next period.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of spacecraft attitude control, in particular to a high-precision magnetic control and thruster combined attitude control method with optimal energy consumption, which is suitable for all spacecrafts that realize high-precision three-axis attitude stabilization by combining magnetic control and thrust control, such as gravity field measurement satellites and high-quietness scientific exploration satellites, etc. The method can also be applied to all other spacecrafts that need magnetic control or specific working modes of spacecrafts. BACKGROUND

[0002] The satellite for precise gravity field measurement task has a load configuration including an extremely high-precision electrostatic suspension accelerometer. In order to reduce the influence of satellite disturbance on the measurement of the accelerometer, the attitude control system excludes rotating and shaking components such as momentum wheels, liquid fuel, etc. from the design, uses cold gas thrusters and continuous voltage type magnetic torque devices as the actuating mechanisms for attitude control, and uses a control method that optimizes the consumption of working medium to ensure that the satellite can normally work to meet the design life requirement. During the development process, on the basis of meeting the technical indicators in the early stage, the researchers made technical innovations and achieved the effect of greatly improving the control precision while further reducing the cold gas working medium consumption level. SUMMARY

[0003] The technical problem solved by the present application is to overcome the shortcomings of the prior art and provide a high-precision magnetic control and thruster combined attitude control method with optimal energy consumption. Unlike the previous feedback control design idea, the method combines feedback control law with direct compensation of satellite residual magnetic moment through filtering estimation, which is a novel design method combining feedback control with output feedforward compensation. The residual magnetic moment is converted into part of the control magnetic moment, avoiding the inefficiency of feedback control and the indirect loss of torque feedforward method, and effectively improving the satellite attitude control precision.

[0004] The technical solution of the present application is a high-precision magnetic control and thruster combined attitude control method with optimal energy consumption, which controls in each cycle according to the following steps:

[0005] First, the expected control moment is calculated using the PD control law according to the attitude control quantity and the attitude angular velocity control quantity, and the moment is compensated in combination with the gravity gradient moment estimation. At the same time, the satellite residual magnetic moment is estimated by filtering according to the angular velocity measurement information in combination with the gravity gradient moment estimation and the theoretical magnetic control moment estimation of the previous cycle;

[0006] Then, the three-axis control of the satellite is distributed according to the satellite magnetic control axis / thrust control axis setting, the direction of the thrust control is calculated for jet phase plane protection, and the jet pulse width is output when the phase plane threshold is triggered. The direction of the magnetic control is calculated for magnetic moment distribution according to the geomagnetic field information and the expected control moment, and the magnetic moment is compensated and limited in amplitude in combination with the above residual magnetic moment filtering estimation result, and the three-axis magnetic moment output is performed.

[0007] Finally, the theoretical magnetic control torque is estimated based on the output of the three-axis magnetic moment, and used for the calculation of the next cycle.

[0008] Preferably, the satellite magnetocontrol axis / thrust control axis is fixed throughout the entire control process.

[0009] Preferably, the allocation of satellite three-axis control includes:

[0010] Two of the three axes of the satellite are magnetically controlled, and the other axis is a thrust control axis. The axis with the weaker overall magnetic control effect is selected as the thrust control axis. The weaker overall magnetic control effect means that the satellite's inertia along this axis is relatively large compared to the other two axes, or that the direction of the geomagnetic field is closer to this axis than the other two axes for most of the time during the entire orbital period.

[0011] Preferably, the perturbation magnetic moment distribution is calculated as follows:

[0012] First, determine which two of the pitch, yaw, and roll axes are magnetic control axes, and which third axis is the thrust control axis. Based on the component values ​​of the desired control torque after torque compensation in the two magnetic control axes, calculate two candidate magnetic moment values. Determine the magnetic moment of the thrust control axis based on whether the signs of the two candidate magnetic moment values ​​are the same.

[0013] Calculate the magnetic moments of the two magnetic control shafts based on the magnetic moments of the thrust control shaft mentioned above.

[0014] Preferably, two candidate magnetic moment values ​​M are calculated. i1 =-T cj / B k M i2 =T ck / B j j and k represent two magnetic control axes, and i represents the thrust control axis; where T cj T ck These represent the calculated values ​​of the desired control torque after torque compensation in the j and k components of the magnetic control axis, respectively. j B k Let M be the components of the geomagnetic induction intensity along the magnetostriction axes j and k in the body coordinate system. i1 With M i2 If the signs are the same, then the magnetic moment M of the thrust control shaft i is... i The calculation formula is: M i =sgn(M i1 min(M) i1 M i2 Otherwise, let M i =0;

[0015] Then calculate the magnetic moment M of the magnetically controlled axis k. k =(T cj +M i Bk ) / B i , the magnetic control axis j the magnetic moment M j = (M i B j -T ck ) / B i , wherein B i is the component value of the geomagnetic induction intensity in the body coordinate system of the thrust control axis.

[0016] Preferably, the expected control moment T excomp after the torque compensation is ex = T g + T ex ;

[0017] wherein T g is the expected control moment, the gravity gradient torque estimation T 2 = 3ω0 BO (E x JE), ω0is the satellite orbit angular velocity value, J is the 3*3 dimension satellite inertia matrix value, E is the 3rd column vector of the satellite body system relative to the orbit system conversion matrix A rem , and x represents the vector cross product calculation.

[0018] Preferably, the filtering estimation of the satellite residual magnetic moment includes:

[0019] Calculating the residual magnetic moment estimation T lst = (Jω - Jω g ) / Δt - T c -T lst , wherein ω and ω c are the three-axis angular velocity of the satellite body relative to the inertial system in the current period and the previous period respectively, Δt is the period length, T v is the theoretical magnetic control moment estimation calculated in the previous period;

[0020] Calculating the magnetic field direction vector B x = B / B|, wherein B = [B y B z B T ] rem is the three-axis component representation of the geomagnetic induction intensity vector in the satellite body coordinate system, and |·| represents the vector length calculation of the vector;

[0021] If the previous period does not perform thrust control, the filtering update value of the satellite residual magnetic moment is calculated as:

[0022] M remlst = M Mrem + K remlst ((M v x B v ) x B rem - (T v) / B),

[0023] where M remlst is the upper cycle calculated value of the satellite residual magnetic moment, K Mrem is the filter coefficient value;

[0024] If the upper cycle has been thrust controlled, the satellite residual magnetic moment maintains the upper cycle estimated value:

[0025] M rem = M remlst ;

[0026] The residual magnetic moment value M remlst = M rem for the next cycle calculation is updated.

[0027] Preferably, the magnetic moment compensation and clipping are as follows:

[0028] The three-axis magnetic moment calculated by the small perturbation magnetic moment distribution is written as a vector M c = [M x M y M z ] T , and the three-axis magnetic moment after the magnetic moment compensation is calculated as M out = M c - M rem ; M rem is the residual magnetic moment filter estimation result.

[0029] Then, the three-axis magnetic moment after the compensation M out is subjected to clipping processing and outputted as M outlim , that is, if the absolute value of the magnetic moment of a certain axis does not exceed the maximum magnetic moment value allowed for the axis, the magnetic moment value is outputted, otherwise, the magnetic moment is outputted as the maximum magnetic moment value allowed for the axis while keeping the sign of the magnetic moment.

[0030] Preferably, the theoretical magnetic control moment estimation is as follows:

[0031] T c = M outlim × B;

[0032] M outlim is the value of the three-axis magnetic moment after the clipping processing, and B is the representation of the geomagnetic induction intensity vector in the satellite body coordinate system.

[0033] A computer readable storage medium, the computer readable storage medium stores a computer program, the computer program is executed by a processor to implement the steps of the method of any one of claims 1-8.

[0034] A high-precision magnetic control and thruster combined attitude control method with optimal energy consumption, which is applied to gravity field measurement satellites, high-quietness scientific exploration satellites, all spacecrafts needing magnetic control or specific working modes of spacecrafts.

[0035] Compared with the prior art, the application has the following advantages:

[0036] (1) The application converts the residual magnetic moment into a part of the control moment by filtering and estimating the satellite residual magnetic moment and direct compensation, avoids the inefficiency of feedback control and the indirect loss of torque feedforward method, effectively improves the satellite attitude control precision, and reduces the consumption of propellant;

[0037] (2) The application fixes the magnetic control axis and the thrust control axis, ensures the continuity of magnetic control, avoids the disturbance of magnetic control / jet switching to the satellite attitude, further reduces the consumption of propellant, and further improves the satellite attitude control precision;

[0038] (3) The application performs low-propellant consumption and high-precision satellite attitude control based on magnetic control and thruster, and can still realize high-precision satellite attitude control when the momentum wheel or CMG and other actuators fail to be used, or other actuators are not used or not configured due to low-power consumption requirements or low-cost requirements, greatly improves the satellite operation reliability, and prolongs the service life of the satellite. BRIEF DESCRIPTION OF DRAWINGS

[0039] Figure 1 A flowchart of the application, a high-precision magnetic control and thruster combined attitude control method with optimal energy consumption; DETAILED DESCRIPTION

[0040] A high-precision magnetic control and thruster combined attitude control method with optimal energy consumption, comprising the following steps:

[0041] (1) The expected control moment is calculated according to the attitude control quantity and the attitude angular velocity control quantity using the PD control law, and the moment compensation is performed in combination with the gravity gradient moment estimation;

[0042] Calculate the expected control moment Wherein K p and K d are 3*3-dimensional control coefficients, and are three-axis Euler attitude angles and three-axis Euler attitude angular velocities, respectively.

[0043] Calculate the gravity gradient moment estimation T g = 3ω0 2 (E x JE), wherein ω0 is the satellite orbit angular velocity value, J is the 3*3-dimensional satellite inertia matrix value, and E is the conversion matrix of the satellite body system relative to the orbit system A BOThe third column vector, × indicates the vector cross product calculation;

[0044] Calculate the desired control torque T after torque compensation excomp =T ex +T g .

[0045] (2) Based on the angular velocity measurement information, combined with the gravity gradient torque estimation and the theoretical magnetic control torque estimation of the previous period, the satellite's remanent magnetic moment is filtered and estimated.

[0046] Calculate the estimated remanent torque T rem =(Jω-Jω) lst ) / Δt-T g -T c , where ω and ω lst These are the current and previous period values ​​(ω) of the satellite's three-axis angular velocities relative to the inertial frame. lst The initial value is ω (initial value), Δt is the cycle duration of this algorithm call (this cycle refers to the satellite time when the algorithm is called this time - the satellite time when the algorithm was called last time; for satellites, this cycle is usually called the control cycle and is a constant value), T c This is an estimate of the theoretical magnetic control torque calculated in the previous cycle;

[0047] Calculate the magnetic field direction vector B v =B / |B|, where B = [B x B y B z ] T The three-axis component representation of the geomagnetic induction intensity vector in the satellite body coordinate system, |·| represents the vector length calculated for the vector;

[0048] If thrust control was not performed in the previous cycle, calculate the filtered update value of the satellite's remanent magnetic moment:

[0049] M rem =M remlst +K Mrem ((M remlst ×B v )×B v -(T rem ×B v ) / |B|),

[0050] Where M remlst This is the calculated value of the satellite's remanent magnetic moment over the previous period, K. Mrem These are the filter coefficient values;

[0051] If thrust control was used in the previous cycle, the satellite's remanent magnetic moment will remain at the estimated value from the previous cycle.

[0052] M rem =M remlst ;

[0053] Update the residual magnetic moment value M for the next cycle calculation remlst = M rem .

[0054] (3) According to the satellite magnetic control axis / thrust control axis setting, the satellite three-axis control is distributed, the direction of thrust control is calculated in the following step (4), and the direction of magnetic control is calculated in the following step (5);

[0055] In step (3), the satellite magnetic control axis / thrust control axis is fixedly set and the three-axis control is distributed as follows:

[0056] Set two axes of the satellite three-axis as the magnetic control axis, perform the magnetic control calculation according to step (5), and set the other axis as the thrust control axis, perform the thrust control calculation according to step (4); when setting the magnetic control axis / thrust control axis, generally, the axis with weaker overall magnetic control effect should be set as the thrust control axis, for example, the inertia of the satellite along the axis is relatively large, or the direction of the geomagnetic field is more close to the axis most of the time during the entire orbit period; the set magnetic control axis / thrust control axis is fixed and unchanged during the entire control process, so as to ensure the continuity of the magnetic control.

[0057] (4) The direction of thrust control is calculated for jet phase plane protection, and jet pulse width output is performed when the phase plane threshold is triggered; the algorithm can refer to the book “Satellite Attitude Dynamics and Control” edited by Academician Tu Shancheng (Aerospace Press, 2001) page 442.

[0058] (5) The direction of magnetic control is calculated according to the geomagnetic field information and the expected control torque, and the residual magnetic moment estimation result in the above step (2) is combined to perform magnetic moment compensation and limiting, and three-axis magnetic moment output is performed;

[0059] The small perturbation magnetic moment distribution calculation in step (5) is as follows:

[0060] a. If the pitch axis and the yaw axis are the magnetic control axes, then:

[0061] First, calculate two alternative magnetic moment values M x1 = -T cy / B z , M x2 = T cz / B y , wherein T cy , T cz are the pitch axis and yaw axis component values of the expected control torque after the torque compensation calculated in step (1), B y , B z are the pitch axis and yaw axis component values of the geomagnetic induction intensity in the body coordinate system, and if M x1 and M x2If the same sign, then the roll axis magnetic moment M x The calculation formula is: M x x1 )min(M x1 , M x2 ); otherwise, let M x = 0;

[0062] Then calculate the yaw axis magnetic moment M z = (T cy + M x B z ) / B x , and the pitch axis magnetic moment M y = (M x B y -T cz ) / B x , where B x is the component value of the geomagnetic induction in the roll axis of the body coordinate system.

[0063] b. If the roll axis and the pitch axis are the magnetic control axes, then:

[0064] First, calculate two alternative magnetic moment values M z1 = -T cx / B y , M z2 = T cy / B x , where T cx , T cy are the component values of the expected control moment after torque compensation in the roll axis and the pitch axis respectively, B x , B y are the component values of the geomagnetic induction in the roll axis and the pitch axis of the body coordinate system respectively, if M z1 and M z2 are of the same sign, then the yaw axis magnetic moment M z The calculation formula is M z = sgn(M z1 )min(M z1 , M z2 ); otherwise, let M z = 0;

[0065] Then calculate the pitch axis magnetic moment M y = (T cx + M z B y ) / B z , and the roll axis magnetic moment M x = (M z B x -T cy ) / B z , where B z ​Bz is the component value of geomagnetic induction intensity in the yaw axis of the body coordinate system.

[0066] c. If the roll axis and the yaw axis are the magnetic control axes, then:

[0067] First, calculate two alternative magnetic moment values M y1 = T cx / B z , M y2 = -T cz / B x , where T cx , T cz are the component values of the moment-compensated expected control moment in the roll axis and the yaw axis calculated in step (1), and B x , B z are the component values of the geomagnetic induction intensity in the roll axis and the yaw axis of the body coordinate system, respectively. If M y1 and M y2 are of the same sign, then the pitch axis magnetic moment M y is calculated as M y = sgn(M y1 )min(M y1 , M y2 ); otherwise, let M y = 0.

[0068] Then, calculate the roll axis magnetic moment M x = (T cz + M y B x ) / B y and the yaw axis magnetic moment M z = (M y B z - T cx ) / B y , where B y is the component value of the geomagnetic induction intensity in the pitch axis of the body coordinate system.

[0069] The magnetic moment compensation and amplitude limiting output in step (5) are as follows:

[0070] Write the three-axis magnetic moment calculated in step (4) in vector form M c = [M x M y M z ] T , and calculate the three-axis magnetic moment after compensation M out = M c - M rem .

[0071] Then, perform amplitude limiting processing on the compensated three-axis magnetic moment M out and output M outlimThat is, if the absolute value of the magnetic moment of an axis does not exceed the maximum magnetic moment value allowed for the axis, the magnetic moment value is output, otherwise the maximum magnetic moment value allowed for the axis is output based on the sign of the magnetic moment.

[0072] (6) Theoretical magnetic control moment estimation is performed according to the result of three-axis magnetic moment calculation in step (5), which is used for calculation in the next period. The theoretical magnetic control moment estimation is as follows:

[0073] T c = M outlim × B.

[0074] Example 1:

[0075] Taking a low-orbit satellite with an orbit inclination near 90 degrees as an example, as shown in FIG. 1, the specific steps of the present application are as follows: Figure 1

[0076] (1) The expected control moment is calculated using a PD control law according to the attitude control quantity and the attitude angular velocity control quantity, and the moment compensation is performed in combination with the gravity gradient moment estimation. First, the expected control moment is calculated using a PD control law where K p and K d are 3*3-dimensional control coefficients, which can be set according to the general PD control design method, and Φ and are the three-axis Euler attitude angle and the three-axis Euler attitude angular velocity, respectively, which are provided by the satellite attitude determination algorithm; then the gravity gradient moment estimation T g = 3ω0 2 (E×JE) is calculated, where ω0 is the satellite orbit angular velocity value, which is given by satellite orbit calculation, J is the 3*3-dimensional satellite inertia matrix value, E is the third column vector of the satellite body system relative to the orbit system conversion matrix A BO , which is provided by the satellite attitude determination algorithm, and × represents vector cross multiplication calculation; subsequently, the expected control moment after moment compensation T excomp = T ex + T g is calculated.

[0077] (2) The satellite residual magnetic moment is estimated by filtering according to the angular velocity measurement information in combination with the gravity gradient moment estimation and the theoretical magnetic control moment estimation in the last period; first, the residual magnetic moment estimation T rem = (Jω-Jω lst ) / Δt-T g -T c is calculated, where ω and ω lst are the three-axis angular velocity of the satellite body relative to the inertial system in the current period and in the last period, respectively, which are provided by the satellite attitude determination algorithm, Δt is the period length called by the algorithm, and T c is the theoretical magnetic control moment estimation calculated in the last period; then the magnetic field direction vector B​v = B / |B|, where B = [B x B y B z ] T is the three-axis component of geomagnetic induction intensity vector in satellite body coordinate system, given by satellite magnetometer measurement or geomagnetic field calculation, | · | represents the vector length calculation of vector; then if no thrust control is performed in the last period, the filtered update value of satellite residual magnetic moment is calculated:

[0078] M rem = M remlst + K Mrem ((M remlst x B v ) x B v - (T rem x B v ) / |B|),

[0079] where M remlst is the calculated value of satellite residual magnetic moment in the last period, K Mrem is the filter coefficient value, generally not more than 0.01, to ensure filter stability;

[0080] If thrust control is performed in the last period, the satellite residual magnetic moment maintains the estimated value M rem = M remlst in the last period; finally, the residual magnetic moment value M remlst = M rem for the next period calculation is updated.

[0081] (3) According to the setting of satellite magnetic control axis / thrust control axis, the satellite three-axis control is distributed, the direction of thrust control is calculated according to step (4), and the direction of magnetic control is calculated according to step (5); two axes of the satellite three-axis are set as magnetic control axes, magnetic control calculation is performed according to step (5), and the other axis is set as thrust control axis, thrust control calculation is performed according to step (4); when setting the magnetic control axis / thrust control axis, the axis with generally weaker magnetic control effect should be set as the thrust control axis, for example, the inertia of the satellite along the axis is relatively large, or the direction of the geomagnetic field is more close to the axis most of the time during the entire orbit period, etc. Since the orbit inclination of the satellite is near 90 degrees, the direction of the geomagnetic field is mainly located near the plane spanned by the roll axis and the yaw axis, and the satellite configuration is long strip type, the inertia of the roll axis is obviously smaller than that of the other two axes, so the yaw axis is set as the thrust control axis, and the roll axis and the pitch axis are set as the magnetic control axes; the set magnetic control axis / thrust control axis is fixed and unchanged during the entire control process, to ensure the continuity of magnetic control.

[0082] (4) Perform jet phase plane protection calculation for the direction of thrust control, and output jet pulse width if necessary; the specific algorithm is referred to in the book Attitude Dynamics and Control of Satellite (Aerospace Publishing House, 2001) edited by Academician Shancheng Tu, page 442.

[0083] (5) Perform small perturbation magnetic moment distribution calculation for the direction of magnetic control according to geomagnetic field information and expected control moment, and perform magnetic moment compensation and limiting according to the residual magnetic moment estimation result of step (2) above, to output three-axis magnetic moment; since the satellite sets the roll axis and the pitch axis as the magnetic control axes, first calculate two candidate magnetic moment values M z1 = -T cx / B y , M z2 = T cy / B x , where T cx , T cy are the roll axis and pitch axis component values of the expected control moment after moment compensation calculated in step (1), B x , B y are the roll axis and pitch axis component values of the geomagnetic induction intensity in the body coordinate system, if M z1 and M z2 are of the same sign, then the yaw axis magnetic moment M z is calculated as M z = sgn(M z1 )min(M z1 , M z2 ), otherwise M z = 0; then calculate the pitch axis magnetic moment M y = (T cx + M z B y ) / B z , and the roll axis magnetic moment M x = (M z B x - T cy ) / B z , where B z is the yaw axis component value of the geomagnetic induction intensity in the body coordinate system; subsequently calculate the three-axis magnetic moment M out = M c - M rem after moment compensation, where M c = [M x M y M z ] T is the vector form of the three-axis magnetic moment calculated in step (4); finally, perform limiting processing on the compensated three-axis magnetic moment M out and output M outlimThat is, if the absolute value of the magnetic moment of a certain axis does not exceed the maximum allowable magnetic moment value of that axis, then the output is based on that magnetic moment value; otherwise, the output is based on the maximum allowable magnetic moment value of that axis while retaining the sign of the magnetic moment.

[0084] (6) Based on the results of the triaxial magnetic moment calculation in step (5), perform a theoretical magnetic control torque estimate for the next cycle calculation; the specific formula is T. c =M outlim ×B;

[0085] (7) Repeat the process throughout the entire operation of the satellite. Figure 1 The diagram shows a high-precision attitude control process that optimizes energy consumption by combining magnetic control and thruster, thereby achieving high-precision attitude control with low propellant consumption and ensuring that the satellite can safely and effectively complete its scientific exploration mission.

[0086] The above description is only the best specific embodiment of the present invention, but the protection scope of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the protection scope of the present invention.

[0087] The contents not described in detail in this specification are common knowledge to those skilled in the art.

Claims

1. A high-precision magnetic control and thruster combined attitude control method with optimal energy consumption, characterized in that The control is performed in each cycle according to the following steps: First, the expected control moment is calculated according to the attitude control amount and the attitude angular velocity control amount using a PD control law, and the moment compensation is combined with the gravity gradient moment estimation, and the satellite residual magnetic moment is estimated by filtering according to the angular velocity measurement information combined with the gravity gradient moment estimation and the theoretical magnetic control moment estimation of the last cycle; Then, the satellite three-axis control is distributed according to the satellite magnetic control axis / thrust control axis setting, the jet phase plane protection calculation is performed on the direction of the thrust control, and the jet pulse width output is performed when the phase plane threshold is triggered, the magnetic moment distribution calculation is performed on the direction of the magnetic control according to the geomagnetic field information and the expected control moment, and the magnetic moment compensation and limiting are combined with the above-mentioned residual magnetic moment filtering estimation result, and the three-axis magnetic moment output is performed; Finally, the theoretical magnetic control moment estimation is performed according to the output result of the three-axis magnetic moment, which is used for the calculation of the next cycle.

2. The united attitude control method according to claim 1, characterized by: The satellite magnetic control axis / thrust control axis is fixedly set in the whole control process.

3. The united attitude control method according to claim 2, characterized by: The distribution of the satellite three-axis control includes: Two axes of the satellite three-axis are set as magnetic control, and the other axis is set as thrust control axis, the axis with weak magnetic control effect is set as thrust control axis, and the weak magnetic control effect is that the inertia of the satellite along the axis is relatively large than the other two axes, or the direction of the geomagnetic field is more close to the axis than the other two axes in most of the whole orbit period.

4. The united attitude control method according to claim 2, characterized by: The magnetic moment distribution calculation is as follows: First, it is judged which two axes of the pitch, yaw and roll three axes are the magnetic control axes, and the third axis is the thrust control axis; according to the component values of the expected control moment after the moment compensation in the two magnetic control axes, two alternative magnetic moment values are calculated, and the magnetic moment of the thrust control axis is determined according to whether the signs of the two alternative magnetic moment values are the same; The magnetic moments of the two magnetic control axes are calculated according to the magnetic moment of the above-mentioned thrust control axis.

5. The joint attitude control method according to claim 4, characterized in that: Calculate two alternative magnetic moment values M i1 = -T cj / B k , M i2 = T ck / B j , j, k represent two magnetic control axes, and i represents a thrust control axis; wherein T cj , T ck are respectively the component values of the expected control moment after compensation of the calculated force moment on the magnetic control axes j, k, B j , B k are respectively the component values of the geomagnetic induction intensity on the magnetic control axes j, k of the body coordinate system, if M i1 and M i2 are of the same sign, then the magnetic moment M i of the thrust control axis i is calculated according to the formula: M i = sgn(M i1 )min(M i1 , M i2 ); otherwise, let M i = 0; The magnetic control axis k magnetic moment M k = (T cj + M i B k ) / B i , the magnetic control axis j magnetic moment M j = (M i B j -T ck ) / B i , wherein B i is the geomagnetic induction intensity in the body coordinate system thrust control axis component value.

6. The method of claim 1, wherein: moment compensated desired control moment T excomp = T ex + T g ; where T ex is the desired control torque, the gravity gradient torque estimate T g = 3ω0 2 (E x JE), ω0is the satellite orbit angular velocity value, J is the 3*3 dimensional satellite inertia matrix value, E is the 3rd column vector of the satellite body frame to orbit frame transformation matrix A BO , and x denotes the vector cross product computation.

7. The method of combined gesture control of claim 6, wherein: The filtering estimation of the satellite residual magnetic moment includes: Calculate the remanent torque estimate T rem = (Jω - Jω lst ) / Δt - T g - T c , where ω and ω lst are the three-axis angular velocity of the satellite body relative to the inertial system, the current period value and the previous period value, Δt is the period length, and T c is the theoretical magnetic control torque estimate calculated in the previous period. Computing the magnetic field direction vector B v = B / |B|, where B = [B x B y B z ] T is the three-axis component representation of the geomagnetic induction intensity vector in the satellite body coordinate system, | · | denotes the vector length calculation on the vector; If the thrust control is not performed in the last cycle, the filtering update value of the satellite residual magnetic moment is calculated: M rem = M remlst + K Mrem ((M remlst x B v ) x B v - (T rem x B v ) / | B | ), where M remlst is the upper periodic value of the satellite remanent moment, K Mrem is the filter coefficient value; If the thrust control is performed in the last cycle, the satellite residual magnetic moment maintains the estimation value in the last cycle: M rem = M remlst ; updating the residual magnetization moment value M for the next cycle calculation remlst = M rem .

8. The method of claim 1, wherein: The magnetic moment compensation and limiting is as follows: The three-axis magnetic moment calculated by the small perturbation magnetic moment distribution is written as a vector form M c = [M x M y M z ] T , the three-axis magnetic moment after compensation of the magnetic moment M out = M c - M rem ; M rem is the remanent magnetic moment filtering estimation result; The compensated three-axis magnetic moment M out is then subjected to limiting processing and output M outlim , i.e. if the absolute value of a certain axis magnetic moment does not exceed the maximum magnetic moment value allowed for that axis, the magnetic moment value is output, otherwise the maximum magnetic moment value allowed for that axis is output on the basis of retaining the positive or negative sign of the magnetic moment.

9. The method of claim 1, wherein: The theoretical magnetic control moment estimation is as follows: T c = M outlim x B; M outlim is the value after the limiting process for the three-axis magnetic moment, and B is the geomagnetic induction vector in the satellite body coordinate system.

10. The optimal energy consumption high-precision magnetic control and thruster combined attitude control method according to claim 1, characterized in that: It is applied to the gravity field measurement satellite, the high-quietness scientific exploration satellite and all spacecrafts that need magnetic control.

Images