CN116841310A - In-plane small control quantity numerical iteration evaluation method based on orbit phase evolution - Google Patents
In-plane small control quantity numerical iteration evaluation method based on orbit phase evolution Download PDFInfo
- Publication number
- CN116841310A CN116841310A CN202310639091.2A CN202310639091A CN116841310A CN 116841310 A CN116841310 A CN 116841310A CN 202310639091 A CN202310639091 A CN 202310639091A CN 116841310 A CN116841310 A CN 116841310A
- Authority
- CN
- China
- Prior art keywords
- orbit
- satellite
- angle
- calculating
- temporary
- 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.)
- Pending
Links
- 238000011156 evaluation Methods 0.000 title claims abstract description 22
- 230000000694 effects Effects 0.000 claims abstract description 18
- 238000004364 calculation method Methods 0.000 claims abstract description 17
- 238000000034 method Methods 0.000 claims abstract description 16
- 239000013598 vector Substances 0.000 claims description 19
- 230000001174 ascending effect Effects 0.000 claims description 18
- 238000013459 approach Methods 0.000 claims description 9
- 238000012937 correction Methods 0.000 claims description 5
- 239000004243 E-number Substances 0.000 claims description 3
- 238000013213 extrapolation Methods 0.000 claims description 3
- 239000000446 fuel Substances 0.000 abstract description 3
- 238000004458 analytical method Methods 0.000 abstract description 2
- 238000005516 engineering process Methods 0.000 description 2
- 239000000126 substance Substances 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
Landscapes
- Control Of Position, Course, Altitude, Or Attitude Of Moving Bodies (AREA)
Abstract
According to the in-plane small control quantity numerical iteration evaluation method based on the orbit phase evolution, through the orbit quantity before and after satellite control and the high-precision orbit dynamics model, the time integral effect of the in-plane small control quantity is utilized, the phase difference of the controlled satellite orbit is taken as a control effect evaluation criterion, the limited calculation precision of an analysis model is considered, and the effectiveness of the calibration method is ensured by adopting numerical iteration in calibration. The method can improve the control effect evaluation precision of the spacecraft during the small control amount orbit in-plane adjustment, thereby improving the orbit control precision of the satellite, effectively saving the fuel consumption of the satellite, having high reliability, strong operability, easy popularization and use, having certain economic benefit for the on-orbit operation of the spacecraft and having important guiding significance for the task implementation.
Description
Technical Field
The invention belongs to the technical field of spacecraft measurement and control methods, and particularly relates to an in-plane small control quantity numerical iteration evaluation method based on orbit phase evolution.
Background
The on-orbit calibration of the orbit control thruster coefficient is to comprehensively calibrate various factor errors which mainly affect the orbit control effect, such as the thrust error of the orbit control thruster, the satellite attitude error and the like, according to the deviation of the actually measured orbit number and the target orbit number. In general, since the track-controlled thrusters are affected by factors such as on-track running conditions and environment, there is a certain deviation in the thrust amount of the thrusters, and this deviation is often a major factor affecting the track-controlled effect. Therefore, the thrust coefficient of the thruster needs to be calibrated after each track control, so that the calibrated coefficient is used in the subsequent track control, and the track control precision is improved.
With the continuous development of the aerospace technology, advanced propulsion technologies such as electric propulsion, laser propulsion and ion propulsion have the characteristics of higher than impulse and high efficiency, and complex and heavier propulsion related equipment in a chemical propulsion system can be omitted from a satellite, so that the satellite can be provided with less fuel to meet task requirements, and can carry more effective loads, and the satellite is widely favored in some engineering applications in recent years. However, since the thrust of the above type of thruster is usually tens milli-to-hundreds milli-newtons, compared with the thrust of the traditional chemical thruster, the track control amount generated in the same time is not large, and is limited by the track precision, and the traditional calibration of the small control amount by using the average control amount or the speed increment of the track number can result in low control effect evaluation precision and difficult accurate calibration of the performance of the thruster. When the satellite performs in-plane small control amount orbit adjustment, the control effect can be indirectly reflected through the time integral effect, so that the control effect under the in-plane small control amount can be indirectly estimated through orbit phase evolution, and the estimation precision of the satellite during orbit in-plane small control amount adjustment is improved through the time integral effect.
Disclosure of Invention
The invention aims to provide an in-plane small control amount numerical iteration evaluation method based on orbit phase evolution, which can improve the control effect evaluation precision when a spacecraft performs in-plane adjustment of a small control amount orbit.
The technical scheme adopted by the invention is as follows: the in-plane small control quantity numerical iteration evaluation method based on the orbit phase evolution utilizes the time integral effect of the in-plane small control quantity through the orbit quantity before and after satellite control and a high-precision orbit dynamics model, takes the phase difference of the controlled satellite orbit as a control effect evaluation criterion, and adopts numerical iteration in calibration to ensure the effectiveness of the calibration method.
The invention is also characterized in that:
the in-plane small control quantity numerical iteration evaluation method based on the orbit phase evolution specifically comprises the following steps: determining a satellite orbit at the current moment, controlling a pre-control satellite orbit at the middle moment, setting an iteration calibration initial value, calculating a temporary orbit latitude amplitude angle and a post-control orbit latitude amplitude angle, calculating a phase evolution quantity, calculating a control speed increment correction quantity according to the orbit phase evolution quantity, calculating an actual control speed increment, calculating a position and a speed vector of the pre-control satellite at the middle moment, calculating a position vector and a speed vector of the temporary orbit satellite, calculating a temporary orbit root number of the satellite, calculating a temporary orbit latitude amplitude angle based on a high-precision dynamics model, calculating a phase difference between the temporary orbit and a post-control satellite orbit, judging the phase difference to meet precision requirements, and calibrating a thruster coefficient.
The in-plane small control quantity numerical iteration evaluation method based on the orbit phase evolution specifically comprises the following steps:
step 1, determining the current T e Time of day satellite orbit parameters, including the half length of the satellite orbitShaft a e Eccentricity e e Inclination angle i e The ascending intersection point is right through omega e Near-site amplitude angle omega e Angle of closest point M e The method comprises the steps of carrying out a first treatment on the surface of the Determining a control intermediate time T s Is controlled by the satellite orbit parameter, including the semi-long axis a of the satellite orbit s Eccentricity e s Inclination angle i s The ascending intersection point is right through omega s Near-site amplitude angle omega s Angle of closest point M s ；
Step 2, setting an iteration calibration initial value including the time T of the temporary track m ＝T s Half major axis a m ＝a s Eccentricity e m ＝e s Inclination angle i m ＝i s The ascending intersection point is right through omega m ＝Ω s Near-site amplitude angle omega m ＝ω s Angle of closest point M m ＝M s Actual control speed increment Δv=0;
step 3, calculating a temporary orbit latitude amplitude angle u m And controlling the latitude amplitude angle u of the rear track e ；
Wherein arctan2 is an arctangent calculation function, E m Is the temporary track approach point angle E e The track is controlled to be close to the point angle;
step 4, calculating a phase evolution quantity delta u;
Δu＝u e -u m
step 5, calculating a control speed increment correction quantity delta v according to the orbit phase evolution quantity delta u;
wherein ,fm Is the true near point angle, calculated by the following formula:
step 6, calculating an actual control speed increment Deltav;
Δv＝Δv+δv
step 7, calculating the satellite position vector before controlling the middle moment under the J2000.0 coordinate systemAnd velocity vector
wherein ,F1 (T s ,a s ,e s ,i s ,Ω s ,ω s ,M s ) Based on satellite orbit time T s Semi-major axis a s Eccentricity e s Inclination angle i s The ascending intersection point is right through the meridian omega s Amplitude angle omega of near-spot s And a mean angle of approach M s Calculating satellite position vectorsAnd velocity vector->Is a function of (2);
step 8, calculating the temporary orbit satellite position vector under the J2000.0 coordinate systemAnd velocity vector->
Wherein, |x| is the modulo calculation of the vector;
step 9, according to the time T of the temporary orbit of the satellite m Position and locationSpeed->Updating and calculating satellite temporary orbit root number semi-long axis a m Eccentricity e m Inclination angle i m The ascending intersection point is right through omega m Near-site amplitude angle omega m Angle of closest point M m ；
wherein ,based on satellite orbit time T m Position->Speed->Calculating a function of the number of satellite orbits;
step 10, extrapolating the temporary orbit root number based on the high-precision dynamics model, and calculating T e Number of temporary orbits at a moment, the parameters including the semimajor axis a of the satellite orbit me Eccentricity e me Inclination angle i me The ascending intersection point is right through omega me Near-site amplitude angle omega me Angle of closest point M me ；
[T e ,a me ,e me ,i me ,Ω me ,ω me ,M me ]＝F 3 (T m ,a m ,e m ,i m ,Ω m ,ω m ,M m )
wherein ,F3 (T m ,a m ,e m ,i m ,Ω m ,ω m ,M m ) The method is a function for extrapolation of the satellite orbit number based on a high-precision dynamics model;
step 11, calculating T e Moment temporary orbit latitude amplitude angle u me ；
wherein ,Eme Is the temporary track near point angle;
step 12, calculating the phase difference delta u between the temporary orbit and the controlled satellite orbit;
δu＝|u me -u e |
step 13, judging whether the phase difference delta u meets the precision requirement, and if the phase difference delta u meets the requirement, namely delta u < delta, wherein delta is an iterative calculation threshold value, performing step 14; otherwise, turning to step 3;
step 14, calibrating the thruster coefficient;
wherein ,Δvs Is the theoretical velocity increment, k s Is the thruster coefficient, deltav s and ks Obtained from the derailment control parameters.
In step 3, the temporary track is close to the point angle E m And controlling the track approaching point angle E e From the equationAnd (5) performing iterative calculation to obtain the product.
Temporary orbit approach Point Angle E in step 11 me From equation E me ＝M me +e me sinE me And (5) performing iterative calculation to obtain the product.
The beneficial effects of the invention are as follows: according to the in-plane small control quantity numerical iteration evaluation method based on the orbit phase evolution, through the orbit quantity before and after satellite control and the high-precision orbit dynamics model, the time integral effect of the in-plane small control quantity is utilized, the phase difference of the controlled satellite orbit is taken as a control effect evaluation criterion, the limited calculation precision of an analysis model is considered, and the effectiveness of the calibration method is ensured by adopting numerical iteration in calibration. The method can improve the control effect evaluation precision of the spacecraft during the small control amount orbit in-plane adjustment, thereby improving the orbit control precision of the satellite, effectively saving the fuel consumption of the satellite, having high reliability, strong operability, easy popularization and use, having certain economic benefit for the on-orbit operation of the spacecraft and having important guiding significance for the task implementation.
Drawings
Fig. 1 is a flow chart of an in-plane small control amount numerical iterative evaluation method based on track phase evolution.
Detailed Description
The invention will be described in detail with reference to the accompanying drawings and detailed description.
The invention provides an in-plane small control quantity numerical iteration evaluation method based on orbital phase evolution, wherein the adopted spacecraft is a near-earth satellite, an electric propulsion system is adopted by the spacecraft, the thrust is 40 millinewtons, semi-long axis control is carried out in the orbit plane, and the control quantity is 5 meters. As shown in fig. 1, the specific steps are as follows:
(1) Determining the current T e The satellite orbit at the moment, the parameter comprises a semi-long axis a of the satellite orbit e Eccentricity e e Inclination angle i e The ascending intersection point is right through omega e Near-site amplitude angle omega e Angle of closest point M e The method comprises the steps of carrying out a first treatment on the surface of the Determining a control intermediate time T s Is controlled by a satellite orbit, the parameter comprises a semi-long axis a of the satellite orbit s Eccentricity e s Inclination angle i s The ascending intersection point is right through omega s Near-site amplitude angle omega s Angle of closest point M s 。
(2) Setting an iteration calibration initial value including the time T of a temporary track m ＝T s Half major axis a m ＝a s Eccentricity e m ＝e s Inclination angle i m ＝i s The ascending intersection point is right through omega m ＝Ω s Near-site amplitude angle omega m ＝ω s Angle of closest point M m ＝M s The actual control speed increment Δv=0.
(3) Calculating a temporary orbit latitude amplitude angle u m And controlling the latitude amplitude angle u of the rear track e ；
Wherein arctan2 is an arctangent calculation function, temporary orbit approach point angle E m And controlling the track approaching point angle E e From the equationAnd (5) performing iterative calculation to obtain the product.
(4) The phase evolution quantity deltau is calculated.
Δu＝u e -u m
(5) And calculating a control speed increment correction quantity delta v according to the track phase evolution quantity.
wherein ,fm Is the true near point angle, calculated by the following formula:
(6) The actual control speed increment deltav is calculated.
Δv＝Δv+δv
(7) Calculating the satellite position vector before controlling the middle moment under the J2000.0 coordinate systemAnd velocity vector->
wherein ,F1 (T s ,a s ,e s ,i s ,Ω s ,ω s ,M s ) Based on satellite orbit time T s Semi-major axis a s Eccentricity e s Inclination angle i s The ascending intersection point is right through the meridian omega s Amplitude angle omega of near-spot s And a mean angle of approach M s Calculating satellite position vectorsAnd velocity vector->Is a function of (2).
(8) Calculating a temporary orbital satellite location vector in a J2000.0 coordinate systemAnd velocity vector->
Where, |x| is the modulo calculation of the vector.
(9) According to the time T of the temporary orbit of the satellite m Position and locationSpeed->Updating and calculating satellite temporary orbit root number semi-long axis a m Eccentricity e m Inclination angle i m The ascending intersection point is right through omega m Near-site amplitude angle omega m Angle of closest point M m ；
wherein ,based on satellite orbit time T m Position->Speed->And calculating a function of the number of satellite orbits.
(10) Extrapolation of temporary orbit root number based on high-precision dynamics model to calculate T e Number of temporary orbits at a moment, the parameters including the semimajor axis a of the satellite orbit me Eccentricity e me Inclination angle i me The ascending intersection point is right through omega me Near-site amplitude angle omega me Angle of closest point M me ；
[T e ,a me ,e me ,i me ,Ω me ,ω me ,M me ]＝F 3 (T m ,a m ,e m ,i m ,Ω m ,ω m ,M m )
wherein ,F3 (T m ,a m ,e m ,i m ,Ω m ,ω m ,M m ) Is a function of extrapolating the number of satellite orbits based on a high-precision dynamics model.
(11) Calculate T e Moment temporary orbit latitude amplitude angle u me ；
Wherein the temporary track is close to the point angle E me From equation E me ＝M me +e me sinE me And (5) performing iterative calculation to obtain the product.
(12) And calculating the phase difference delta u of the temporary orbit and the controlled satellite orbit.
(13) Judging whether the phase difference delta u meets the precision requirement, if so, namely delta u < delta, wherein delta is an iterative calculation threshold value which can be manually selected according to the requirement, performing step 14, otherwise, turning to step 3.
(14) Calibrating the coefficient of the thruster;
wherein ,Δvs Is the theoretical velocity increment, k s Is the thruster coefficient, deltav, used in this control s and ks Obtained from the derailment control parameters.
Claims (5)
1. The in-plane small control quantity numerical iteration evaluation method based on the orbit phase evolution is characterized in that the time integral effect of the in-plane small control quantity is utilized through the orbit quantity before and after satellite control and a high-precision orbit dynamics model, the phase difference of the controlled satellite orbit is used as a control effect evaluation criterion, and numerical iteration is adopted in calibration to ensure the effectiveness of the calibration method.
2. The method for iteratively evaluating the in-plane small control quantity value based on the orbit phase evolution according to claim 1, which is characterized by comprising the following steps: determining a satellite orbit at the current moment, controlling a pre-control satellite orbit at the middle moment, setting an iteration calibration initial value, calculating a temporary orbit latitude amplitude angle and a post-control orbit latitude amplitude angle, calculating a phase evolution quantity, calculating a control speed increment correction quantity according to the orbit phase evolution quantity, calculating an actual control speed increment, calculating a position and a speed vector of the pre-control satellite at the middle moment, calculating a position vector and a speed vector of the temporary orbit satellite, calculating a temporary orbit root number of the satellite, calculating a temporary orbit latitude amplitude angle based on a high-precision dynamics model, calculating a phase difference between the temporary orbit and a post-control satellite orbit, judging the phase difference to meet precision requirements, and calibrating a thruster coefficient.
3. The method for iteratively evaluating the in-plane small control quantity value based on the orbit phase evolution according to claim 1, comprising the steps of:
step 1, determining the current T e The satellite orbit parameters of the moment include the semi-long axis a of the satellite orbit e Eccentricity e e Inclination angle i e The ascending intersection point is right through omega e Near-site amplitude angle omega e Angle of closest point M e The method comprises the steps of carrying out a first treatment on the surface of the Determining a control intermediate time T s Is controlled by the satellite orbit parameter, including the semi-long axis a of the satellite orbit s Eccentricity e s Inclination angle i s The ascending intersection point is right through omega s Near-site amplitude angle omega s Angle of closest point M s ；
Step 2, setting an iteration calibration initial value including the time T of the temporary track m ＝T s Half major axis a m ＝a s Eccentricity e m ＝e s Inclination angle i m ＝i s The ascending intersection point is right through omega m ＝Ω s Near-site amplitude angle omega m ＝ω s Angle of closest point M m ＝M s Actual control speed increment Δv=0;
step 3, calculating a temporary orbit latitude amplitude angle u m And controlling the latitude amplitude angle u of the rear track e ；
Wherein arctan2 is an arctangent calculation function, E m Is the temporary track approach point angle E e The track is controlled to be close to the point angle;
step 4, calculating a phase evolution quantity delta u;
Δu＝u e -u m
step 5, calculating a control speed increment correction quantity delta v according to the orbit phase evolution quantity delta u;
wherein ,fm Is the true near point angle, calculated by the following formula:
step 6, calculating an actual control speed increment Deltav;
Δv＝Δv+δv
step 7, calculating the satellite position vector before controlling the middle moment under the J2000.0 coordinate systemAnd velocity vector->
wherein ,F1 (T s ,a s ,e s ,i s ,Ω s ,ω s ,M s ) Based on satellite orbit time T s Semi-major axis a s Eccentricity e s Inclination angle i s The ascending intersection point is right through the meridian omega s Amplitude angle omega of near-spot s And a mean angle of approach M s Calculating satellite position vectorsAnd velocity vector->Is a function of (2);
step 8, calculating the temporary orbit satellite position vector under the J2000.0 coordinate systemAnd velocity vector->
Wherein, |x| is the modulo calculation of the vector;
step 9, according to the time T of the temporary orbit of the satellite m Position and locationSpeed->Updating and calculating satellite temporary orbit root number semi-long axis a m Eccentricity e m Inclination angle i m The ascending intersection point is right through omega m Near-site amplitude angle omega m Angle of closest point M m ；
wherein ,based on satellite orbit time T m Position->Speed->Calculating a function of the number of satellite orbits;
step 10, extrapolating the temporary orbit root number based on the high-precision dynamics model, and calculating T e Number of temporary orbits at a moment, the parameters including the semimajor axis a of the satellite orbit me Eccentricity e me Inclination angle i me The ascending intersection point is right through omega me Near-site amplitude angle omega me Angle of closest point M me ；
[T e ,a me ,e me ,i me ,Ω me ,ω me ,M me ]＝F 3 (T m ,a m ,e m ,i m ,Ω m ,ω m ,M m )
wherein ,F3 (T m ,a m ,e m ,i m ,Ω m ,ω m ,M m ) The method is a function for extrapolation of the satellite orbit number based on a high-precision dynamics model;
step 11, calculating T e Moment temporary orbit latitude amplitude angle u me ；
wherein ,Eme Is the temporary track near point angle;
step 12, calculating the phase difference delta u between the temporary orbit and the controlled satellite orbit;
δu＝|u me -u e |
step 13, judging whether the phase difference delta u meets the precision requirement, and if the phase difference delta u meets the requirement, namely delta u < delta, wherein delta is an iterative calculation threshold value, performing step 14; otherwise, turning to step 3;
step 14, calibrating the thruster coefficient;
wherein ,Δvs Is the theoretical velocity increment, k s Is the thruster coefficient, deltav s and ks Obtained from the derailment control parameters.
4. The method for iterative evaluation of in-plane small control quantity values based on orbital phase evolution according to claim 3, wherein the temporary orbital approach point angle E in step 3 m And controlling the track approaching point angle E e From the equationAnd (5) performing iterative calculation to obtain the product.
5. The method for iterative evaluation of in-plane small control quantity values based on orbital phase evolution according to claim 3, wherein the temporary orbital approach point angle E in step 11 me From equation E me ＝M me +e me sinE me And (5) performing iterative calculation to obtain the product.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310639091.2A CN116841310A (en) | 2023-05-31 | 2023-05-31 | In-plane small control quantity numerical iteration evaluation method based on orbit phase evolution |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310639091.2A CN116841310A (en) | 2023-05-31 | 2023-05-31 | In-plane small control quantity numerical iteration evaluation method based on orbit phase evolution |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116841310A true CN116841310A (en) | 2023-10-03 |
Family
ID=88173382
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310639091.2A Pending CN116841310A (en) | 2023-05-31 | 2023-05-31 | In-plane small control quantity numerical iteration evaluation method based on orbit phase evolution |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116841310A (en) |
-
2023
- 2023-05-31 CN CN202310639091.2A patent/CN116841310A/en active Pending
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106292287B (en) | A kind of UUV path following method based on adaptive sliding-mode observer | |
| CN110794863B (en) | Heavy carrier rocket attitude control method capable of customizing control performance indexes | |
| CN106697333B (en) | A kind of robust analysis method of spacecraft orbit control strategy | |
| Pesch | A practical guide to the solution of real-life optimal control problems | |
| CN105573337B (en) | A kind of braking Closed Loop Guidance method that leaves the right or normal track for meeting reentry angle and voyage constraint | |
| CN113665849B (en) | Autonomous phase control method combining EKF filtering algorithm and neural network | |
| CN110262241B (en) | Spacecraft orbit control method based on Gaussian process prediction control | |
| CN108959734B (en) | Real-time recursion-based solar light pressure moment identification method and system | |
| CN113900448B (en) | Aircraft prediction correction composite guidance method based on sliding mode interference observer | |
| CN112560343B (en) | J2 perturbation Lambert problem solving method based on deep neural network and targeting algorithm | |
| CN112306075A (en) | High-precision off-orbit reverse iterative guidance method | |
| CN112393835B (en) | Small satellite on-orbit thrust calibration method based on extended Kalman filtering | |
| CN111605736B (en) | Earth-moon L2 point transfer orbit optimal error correction point selection method | |
| CN110597274B (en) | SGCMG dynamic frame angular velocity determination method adaptive to attitude redirection | |
| CN116841310A (en) | In-plane small control quantity numerical iteration evaluation method based on orbit phase evolution | |
| CN112363538A (en) | AUV (autonomous underwater vehicle) area tracking control method under incomplete speed information | |
| CN115655285B (en) | Unscented quaternion attitude estimation method for correcting weight and reference quaternion | |
| CN115993777A (en) | Track perturbation model inversion-based diameter-cut joint control decoupling iteration calibration method | |
| CN110955255A (en) | High-precision orbit control attitude maintaining method, system and medium based on CMG | |
| CN112904719B (en) | Annular area tracking control method suitable for underwater robot position | |
| Huang et al. | Pseudospectral method for optimal propellantless rendezvous using geomagnetic Lorentz force | |
| CN115422496A (en) | Combined correction identification method for carrier rocket mass and thrust parameters under thrust fault | |
| CN115113638A (en) | Fuel optimal active drift three-dimensional imaging track control method | |
| CN112455725A (en) | Method for transferring and converting pulse orbit transfer direction to limited thrust orbit | |
| Shen et al. | High-accuracy optimal finite-thrust trajectories for Moon escape |
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 |