CN114889849A - Estimation method for multi-constraint flying-over optimal speed of spacecraft - Google Patents

Estimation method for multi-constraint flying-over optimal speed of spacecraft Download PDF

Info

Publication number
CN114889849A
CN114889849A CN202210737971.9A CN202210737971A CN114889849A CN 114889849 A CN114889849 A CN 114889849A CN 202210737971 A CN202210737971 A CN 202210737971A CN 114889849 A CN114889849 A CN 114889849A
Authority
CN
China
Prior art keywords
flying
spacecraft
speed
constraint
velocity
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.)
Granted
Application number
CN202210737971.9A
Other languages
Chinese (zh)
Other versions
CN114889849B (en
Inventor
罗亚中
严冰
张进
王柯茂
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
National University of Defense Technology
Original Assignee
National University of Defense Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by National University of Defense Technology filed Critical National University of Defense Technology
Priority to CN202210737971.9A priority Critical patent/CN114889849B/en
Publication of CN114889849A publication Critical patent/CN114889849A/en
Application granted granted Critical
Publication of CN114889849B publication Critical patent/CN114889849B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • B—PERFORMING OPERATIONS; TRANSPORTING
    • B64—AIRCRAFT; AVIATION; COSMONAUTICS
    • B64G—COSMONAUTICS; VEHICLES OR EQUIPMENT THEREFOR
    • B64G1/00—Cosmonautic vehicles
    • B64G1/22—Parts of, or equipment specially adapted for fitting in or to, cosmonautic vehicles
    • B64G1/24—Guiding or controlling apparatus, e.g. for attitude control
    • B64G1/242—Orbits and trajectories
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Chemical & Material Sciences (AREA)
  • Combustion & Propulsion (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Navigation (AREA)
  • Control Of Position, Course, Altitude, Or Attitude Of Moving Bodies (AREA)

Abstract

The invention discloses a method for estimating the optimal speed of multi-constraint flying of a spacecraft, which comprises the following steps: s1, establishing a terminal inequality constraint expression of the fly-by patrol; s2, calculating the optimal flying time and the plane orientation of the flying orbit according to the orbital out-of-plane degree of the spacecraft and the target; and S3, giving an initial value of the flying speed according to the flying orbital plane orientation, and establishing a multi-constraint flying optimal speed estimation in fixed time based on the Cohen-Tack condition. The method sets the magnitude and three direction components of the tail end flying speed as 4 unknowns, designs an objective function of an optimization problem as a speed increment under a single pulse, introduces 3 inequality constraints including a flying inspection sunlight angle, a flying relative speed and a track near-location height during flying, lists an 8-dimensional equation set by using a Cohen-Tack condition, obtains a flying optimal speed estimation, changes the flying estimation into a rendezvous estimation, and estimates the transfer cost by using a known track root number difference.

Description

Estimation method for multi-constraint flying optimal speed of spacecraft
Technical Field
The invention relates to the technical field of estimation of a spacecraft flying state, in particular to a method for estimating the optimal speed of multi-constraint flying of a spacecraft.
Background
With the construction of low-orbit giant constellations such as satellite chain satellites and the like, the on-orbit targets in the near-earth space are increased sharply, and the satellite group shielding the sun presents a serious challenge to the perception capability of the future space situation. Due to the development of space-based optical imaging technology, a spatial situation perception system based on a satellite provided with an optical sensor is further emphasized. The multi-spacecraft is utilized to perform out-of-plane flying, the flying can be approached to and a plurality of target satellites can be shot in a short time at a low fuel cost, so that image data of constellation satellites, formation satellites and other cluster satellites can be quickly obtained, and support is provided for subsequent task decision.
In studying one-to-many, many-to-many fly-by patrol sequence problems, efficient single-pair single-fly-by cost estimation is needed to support. In the intersection estimation problem, since the position speeds of the beginning and the end are known, the transfer cost at the corresponding time can be estimated according to the track root number difference. However, in the fly-by mission, since the terminal speed is unknown, the optimal fly-by speed satisfying the terminal speed constraint needs to be estimated, so that the whole fly-by transition cost is minimum, and therefore, estimating the fly-by mission cost under the condition that the fly-by terminal speed cannot be determined in advance is a technical problem to be solved at present.
Disclosure of Invention
The invention aims to provide a method for estimating the optimal speed of multi-constraint flying of a spacecraft, so as to overcome the defects in the prior art.
In order to achieve the purpose, the technical scheme adopted by the invention is as follows:
a method for estimating the optimal speed of multi-constraint flying of a spacecraft comprises the following steps:
s1, establishing a terminal inequality constraint expression of the fly-by patrol;
s2, calculating the optimal flying time and the plane orientation of the flying orbit according to the orbital out-of-plane degree of the spacecraft and the target;
and S3, giving an initial value of the flying speed according to the flying orbital plane orientation, and establishing a multi-constraint flying optimal speed estimation in fixed time based on the Cohen-Tack condition.
Further, the step S1 is specifically:
setting the starting time t of the spacecraft under the geocentric inertial coordinate system 0 And a position velocity at the time of flight t is [ r ] C0 ,v C0 ]And [ r C ,v C ]And the target is at the starting time t 0 And a position velocity at the time of flight t is [ r ] T0 ,v T0 ]And [ r T ,v T ]Setting the flying task as non-deviation flying with position coincidence and speed deviation; the motion of the spacecraft and the target in the set time period before the flying is set to be uniform and linear, so that the relative position vector delta r of the spacecraft and the target can be deduced to be r T –r C With relative velocity vector δ v ═ v T –v C The cross product of (a) is 0, and the sight line direction before the position coincidence is replaced by delta v; setting a sight line direction and a sun position vector r S The lower limit of the included angle of the direction is alpha smin The constraint expression on the viewing sun angle is obtained as follows:
arccos((r S /r S )·(δv/δv))-(π-α smin )≤0 (36)
setting the relative speed of the spacecraft in flying not to exceed the upper limit delta v max ：
δv-δv max ≤0 (37)
Setting the height of the near point of the spacecraft in flying to be not lower than the lower limit r pmin ：
r pmin -a C (1-e C )≤0 (38)
In the formula, a C And e C The semi-major axis and eccentricity of the spacecraft at the moment of flight t.
Further, the step S2 specifically includes the following steps:
s20, setting the spacecraft at the starting time t 0 And the number of tracks at the time of flight t is [ a ] C0 ,e C0 ,i C0 ,Ω C0 ,ω C0 ,M C0 ]And [ a ] C ,e C ,i C ,Ω C ,ω C ,M C ]And the target is at the starting time t 0 And the number of tracks at the time of flight t is [ a ] T0 ,e T0 ,i T0 ,Ω T0 ,ω T0 ,M T0 ]And [ a ] T ,e T ,i T ,Ω T ,ω T ,M T ]A, e, i, omega, M are respectively a semi-major axis, eccentricity, inclination, elevation intersection declination, near arch point latitude argument and near point angle, subscripts C0 and C are the space vehicle at t 0 And the number of tracks at time T, subscripts T0 and T targeting at T 0 And the number of tracks at time t, the track transfer duration Deltat being t-t 0 Calculating out the heterozygosity phi of the two initial orbital planes and the maximum heterozygosity phi allowed by the relative speed constraint according to a spherical triangle formula max Comprises the following steps:
Figure BDA0003709352020000021
in the formula, mu is the universal gravitation constant of the central celestial body;
s21, when phi is less than or equal to phi max The spacecraft allows flying through the target by coplanar maneuvers only, where the optimal flying time Δ t opt The expression of (a) is:
Figure BDA0003709352020000022
in the formula u T Phase, u, of the initial orbital plane intersection of the spacecraft and the target on the target orbit T0 An initial phase of the target, k being a non-negative integer, u T The expression of (a) is:
Figure BDA0003709352020000023
in the formula u T There are two solutions, only the solution on the side of the shadow,
Figure BDA0003709352020000024
and
Figure BDA0003709352020000025
aiming at considering only J 2 The orbit root precession rate under the long-term perturbation comprises the expression of the orbit root precession rate including a spacecraft:
Figure BDA0003709352020000031
in the formula, J 2 Is a perturbation term constant, R E Is the earth's average equatorial radius;
s22, according to phi and phi max And Δ t opt Determining a unit normal vector h of the flying-over orbital plane e The threshold value of the time deviation is delta t threshold When phi is less than or equal to phi max And | Δ t- Δ t opt |≤Δt threshold Considering that the spacecraft can perform coplanar maneuver, the angular edge (i) in the spherical triangle is used C0 -i T0 -u T ) Solving the elevation difference delta omega of the elevation intersection point as follows:
Figure BDA0003709352020000032
when phi is less than or equal to phi max And | Δ t- Δ t opt |>Δt threshold Or phi>φ max In time, the spacecraft needs to be maneuvered to orbit in different planes according to the corner angle (i) in the spherical triangle T0 -u T -φ max ) To solve for the inclination i C And the rise-crossing point declination difference delta omega is:
Figure BDA0003709352020000033
according to Δ Ω ═ Ω T0 –Ω C Thereby obtaining a unit normal vector h e Is composed of
Figure BDA0003709352020000034
Further, the step S3 specifically includes the following steps:
s30, establishing O, r with the target as the origin T The direction being the x-axis, h e Setting the flying speed with respect to a coordinate system of LVLH with the direction of the z-axisMagnitude and three directional components of v C ,γ x ,γ y ,γ z Velocity vector v C And h e Orthogonal, magnitude of flying velocity v C In the xOy plane, the magnitude v of the flying speed is set C The included angle between the positive direction of the y axis is theta epsilon-theta max ,θ max ]，θ max For boundary values of the angle, the three directional components of the velocity can be described by θ by the coordinate system change as:
Figure BDA0003709352020000035
in the formula (I), the compound is shown in the specification,
Figure BDA0003709352020000036
a rotation matrix from an LVLH coordinate system to a geocentric inertia coordinate system;
the magnitude of the velocity is deduced according to the displacement formula of the position velocity and the number of the tracks C Expression (c):
Figure BDA0003709352020000041
in the formula, e Cmax Is the allowable upper limit of the eccentricity of the fly-over orbit;
by discrete v C And the value of theta in the respective upper and lower bounds, and the guess value is:
Figure BDA0003709352020000042
in the formula, N 1 And N 2 For the number of discrete points, α and β are both integers, and θ is guess Carrying in formula (11) to obtain a guess initial value of the direction;
s31, if the set flying trace is estimated by a single pulse, the number of orbit deviations caused by the single pulse are:
Figure BDA0003709352020000043
the magnitude of the velocity increment Δ v of a single pulse is:
Figure BDA0003709352020000044
in the formula (I), the compound is shown in the specification,
Figure BDA0003709352020000045
the average orbit speed of the spacecraft at the initial moment is taken as the average orbit speed;
taking the square of the speed increment as an objective function of the optimal problem, the expression is as follows:
Figure BDA0003709352020000046
since the modulus values in the three unit directions of velocity are 1, the equality constraint is introduced as:
Figure BDA0003709352020000047
establishing a Lagrange equation, wherein the expression is as follows:
L＝f+λ T g+κ T σ (53)
in the formula, λ is an equality constraint multiplier, κ is an inequality constraint multiplier, and according to the kun-tak condition, the partial derivatives of the four unknowns are solved, so that an 8-dimensional equality equation set is obtained:
Figure BDA0003709352020000051
in the formula, relate to a C ,e C ,i C ,Ω C ,ω C For v C ,γ x ,γ y ,γ z The partial derivatives of (2) are expressed by the number of the tracks in terms of position and velocity, and are recorded in terms of h C ＝[h Cx ,h Cy ,h Cz ] T ，e C ＝[e Cx ,e Cy ,e Cz ] T The expression is:
Figure BDA0003709352020000052
wherein h is C For the angular momentum of the fly-over orbit, let γ ═ γ x ,γ y ,γ z ] T The partial derivatives can be derived as:
Figure BDA0003709352020000061
Figure BDA0003709352020000062
Figure BDA0003709352020000063
Figure BDA0003709352020000064
Figure BDA0003709352020000065
Figure BDA0003709352020000066
Figure BDA0003709352020000071
note i x ＝[1,0,0] T ,i y ＝[0,1,0] T ,i z ＝[0,0,1] T Angular momentum vector h C And its modulus value h C About v C The partial derivatives of γ are:
Figure BDA0003709352020000072
Figure BDA0003709352020000073
relative velocity vector δ v and its modulus value δ v with respect to v C The partial derivatives of γ are:
Figure BDA0003709352020000074
eccentricity vector e C About v C The partial derivatives of γ are:
Figure BDA0003709352020000075
based on the initial value [ v ] of the flying speed given in step S3 C ,γ x ,γ y ,γ z ] 0 And the initial values of the equality constraint and inequality constraint multipliers [ lambda, kappa ] 1 ,κ 2 ,κ 3 ] 0 When the flying speed is 0, a Powell-based mixed numerical value iterator iterates to obtain a final flying speed solution;
s32, obtaining the final flying speed solution, generating a phase difference
Figure BDA0003709352020000076
Wherein u is Cf Extrapolating the phase of the spacecraft to the moment of flight for the spacecraft,
Figure BDA0003709352020000077
and
Figure BDA0003709352020000078
for precession of spacecraft with respect to no maneuveringRate deviation, a first order approximation of the relative average angular velocity, can be obtained:
Figure BDA0003709352020000081
the phase correction is carried out by using the semimajor axis deviation, so that the delta u + delta n C Δ t is 0, and the modified semimajor axis expression is:
Figure BDA0003709352020000082
wherein N ∈ [ -1,1] represents the number of turns, and the magnitude of the flying speed after correction is as follows:
Figure BDA0003709352020000083
the resulting flying speed is
Figure BDA0003709352020000084
S33, converting the flying cost estimation into a rendezvous cost estimation, and estimating according to a rendezvous estimation algorithm;
s34, repeating the steps S30-S33, counting the speed increment delta v generated by each iteration, and selecting the minimum value to be recorded as delta v best The corresponding fly-over velocity is solved as [ v C ,γ x ,γ y ,γ z ] best 。
Compared with the prior art, the invention has the advantages that: the invention provides a method for estimating the optimal speed of multi-constraint flying of a spacecraft, which is characterized in that the size and three direction components of the tail end flying speed are set as 4 unknowns, the objective function of the design optimization problem is the speed increment under a single pulse, 3 inequality constraints including the sun angle of flying inspection, the relative speed of flying and the altitude of the orbit near place during flying are introduced, an 8-dimensional equation set is listed by utilizing the Kuen-Tack condition, the optimal speed estimation of flying is obtained by solving the equation set, the flying estimation is changed into intersection estimation, and therefore the transfer cost is estimated by utilizing the known orbit root number difference.
Drawings
In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the embodiments or the prior art will be briefly described below, it is obvious that the drawings in the following description are only some embodiments of the present invention, and for those skilled in the art, other drawings can be obtained according to the drawings without creative efforts.
FIG. 1 is a line-of-sight solar angle constraint schematic diagram of the spacecraft in the flying time of the invention.
FIG. 2 is a flow chart of the method for estimating the multi-constraint flying optimal speed of the spacecraft.
Detailed Description
The preferred embodiments of the present invention will be described in detail below with reference to the accompanying drawings so that the advantages and features of the present invention can be more easily understood by those skilled in the art, and the scope of the present invention will be more clearly and clearly defined.
Referring to fig. 1 and fig. 2, the embodiment discloses a method for estimating a multi-constraint flying optimal speed of a spacecraft, which includes the following steps:
and step S1, establishing a terminal inequality constraint expression for the fly-by patrol.
Specifically, the starting time t of the spacecraft under the geocentric inertial coordinate system is set 0 And a position velocity at the time of flight t is [ r ] C0 ,v C0 ]And [ r C ,v C ]And the target is at the starting time t 0 And a position velocity at the time of flight t is [ r ] T0 ,v T0 ]And [ r T ,v T ]Setting the flying task as non-deviation flying with position coincidence and speed deviation; the motion of the spacecraft and the target in the set time period before the flying is set to be uniform and linear, so that the relative position vector delta r of the spacecraft and the target can be deduced to be r T –r C With relative velocity vector δ v ═ v T –v C Cross product of (d) is 0 and δ v is used to replace the view before the position coincidenceA line direction; setting a sight line direction and sun position vector r S The lower limit of the included angle of the direction is alpha smin The constraint expression on the viewing sun angle is obtained as follows:
arccos((r S /r S )·(δv/δv))-(π-α smin )≤0 (71)
setting the relative speed of the spacecraft in flying not to exceed the upper limit delta v max ：
δv-δv max ≤0 (72)
Setting the height of the near point of the spacecraft in flying to be not lower than the lower limit r pmin ：
r pmin -a C (1-e C )≤0 (73)
In the formula, a C And e C The semi-major axis and eccentricity of the spacecraft at the moment of flight t.
And step S2, calculating the optimal flying time and the plane orientation of the flying orbit according to the orbital out-of-plane degree of the spacecraft and the target.
Specifically, step S2 includes the following steps:
step S20, setting the spacecraft at the starting time t 0 And the number of tracks at the time of flight t is [ a ] C0 ,e C0 ,i C0 ,Ω C0 ,ω C0 ,M C0 ]And [ a ] C ,e C ,i C ,Ω C ,ω C ,M C ]And the target is at the starting time t 0 And the number of tracks at the time of flight t is [ a ] T0 ,e T0 ,i T0 ,Ω T0 ,ω T0 ,M T0 ]And [ a ] T ,e T ,i T ,Ω T ,ω T ,M T ]A, e, i, omega, M are respectively a semi-major axis, eccentricity, inclination, elevation intersection declination, near arch point latitude argument and near point angle, subscripts C0 and C are the space vehicle at t 0 And the number of tracks at time T, subscripts T0 and T targeting at T 0 And the number of tracks at time t, the track transfer duration Deltat being t-t 0 Calculating out the heterozygosity phi of the two initial orbital planes and the maximum heterozygosity phi allowed by the relative speed constraint according to a spherical triangle formula max Comprises the following steps:
Figure BDA0003709352020000091
wherein mu is the gravitational constant of the central celestial body.
Step S21, when phi is less than or equal to phi max The spacecraft allows flying through the target by coplanar maneuvers only, where the optimal flying time Δ t opt The expression of (a) is:
Figure BDA0003709352020000101
in the formula u T Phase, u, of the initial orbital plane intersection of the spacecraft and the target on the target orbit T0 An initial phase of the target, k being a non-negative integer, u T The expression of (a) is:
Figure BDA0003709352020000102
in the formula u T There are two solutions, only the solution on the side of the shadow,
Figure BDA0003709352020000107
and
Figure BDA0003709352020000108
aiming at considering only J 2 The orbit root precession rate under the long-term perturbation comprises the expression of the orbit root precession rate including a spacecraft:
Figure BDA0003709352020000103
in the formula, J 2 Is a perturbation term constant, R E Is the average equatorial radius of the earth.
Step S22, according to phi and phi max And Δ t opt Determining a unit normal vector h of the flying-over orbital plane e The threshold value of the time deviation is delta t threshold When phi is less than or equal to phi max And | Δ t- Δ t opt |≤Δt threshold Considering that the spacecraft can perform coplanar maneuver, the angular edge (i) in the spherical triangle is used C0 -i T0 -u T ) Solving the elevation difference delta omega of the elevation intersection point as follows:
Figure BDA0003709352020000104
when phi is less than or equal to phi max And | Δ t- Δ t opt |>Δt threshold Or phi>φ max In time, the spacecraft needs to be maneuvered to orbit in different planes according to the corner angle (i) in the spherical triangle T0 -u T -φ max ) To solve for the inclination i C And the rise-crossing point declination difference delta omega is:
Figure BDA0003709352020000105
according to Δ Ω ═ Ω T0 –Ω C Thereby obtaining a unit normal vector h e Is composed of
Figure BDA0003709352020000106
And step S3, an initial value of the flying speed is given according to the flying orbital plane orientation, and a multi-constraint flying optimal speed estimation under fixed time is established based on the Cohen-Tack condition.
Specifically, step S3 specifically includes the following steps:
step S30, establishing O, r with the target as the origin T The direction being the x-axis, h e Setting the magnitude of the flying speed and three directional components as v in a relative coordinate system of the LVLH with the direction as a z-axis C ,γ x ,γ y ,γ z Velocity vector v C And h e Orthogonal, magnitude of flying velocity v C In the xOy plane, the magnitude v of the flying speed is set C The included angle theta and the positive direction of the y axis is epsilon[-θ max ,θ max ]，θ max For boundary values of the angle, the three directional components of the velocity can be described by θ by the coordinate system change as:
Figure BDA0003709352020000111
in the formula (I), the compound is shown in the specification,
Figure BDA0003709352020000112
a rotation matrix from an LVLH coordinate system to a geocentric inertia coordinate system;
the magnitude of the velocity is deduced according to the displacement formula of the position velocity and the number of the tracks C Expression (c):
Figure BDA0003709352020000113
in the formula, e Cmax Is the allowable upper limit of the eccentricity of the fly-over orbit;
by discrete v C And the value of theta in the respective upper and lower bounds, and the guess value is:
Figure BDA0003709352020000114
in the formula, N 1 And N 2 For the number of discrete points, α and β are both integers, and θ is guess Carrying in formula (11) to obtain a guess initial value of the direction;
step S31, if the estimated flying trace is set by a single pulse, the number of orbit deviations caused by the single pulse are:
Figure BDA0003709352020000115
the magnitude of the velocity increment Δ v of a single pulse is:
Figure BDA0003709352020000116
in the formula,. DELTA.v a =0.5v C0 Δa/a C0 ，Δv i =v C0 Δi，Δv Ω =v C0 ΔΩsini C0 ，
Figure BDA0003709352020000117
v C0 The average orbit speed of the spacecraft at the initial moment is taken as the average orbit speed;
taking the square of the speed increment as an objective function of the optimal problem, the expression is as follows:
Figure BDA0003709352020000121
since the modulus values in the three unit directions of velocity are 1, the equality constraint is introduced as:
Figure BDA0003709352020000122
establishing a Lagrange equation, wherein the expression is as follows:
L＝f+λ T g+κ T σ (88)
in the formula, λ is an equality constraint multiplier, κ is an inequality constraint multiplier, and according to the kun-tak condition, the partial derivatives of the four unknowns are solved, so that an 8-dimensional equality equation set is obtained:
Figure BDA0003709352020000123
in the formula, relate to a C ,e C ,i C ,Ω C ,ω C For v C ,γ x ,γ y ,γ z The partial derivatives of (2) are expressed by the number of the tracks in terms of position and velocity, and are recorded in terms of h C ＝[h Cx ,h Cy ,h Cz ] T ，e C ＝[e Cx ,e Cy ,e Cz ] T The expression is:
Figure BDA0003709352020000131
wherein h is C For the angular momentum of the fly-over orbit, let γ ═ γ x ,γ y ,γ z ] T The partial derivatives can be derived as:
Figure BDA0003709352020000132
Figure BDA0003709352020000133
Figure BDA0003709352020000134
Figure BDA0003709352020000135
Figure BDA0003709352020000141
Figure BDA0003709352020000142
Figure BDA0003709352020000143
note i x ＝[1,0,0] T ,i y ＝[0,1,0] T ,i z ＝[0,0,1] T Angular momentum vector h C And its modulus value h C About v C The partial derivatives of γ are:
Figure BDA0003709352020000144
Figure BDA0003709352020000145
relative velocity vector δ v and its modulus value δ v with respect to v C And the partial derivatives of γ are:
Figure BDA0003709352020000146
eccentricity vector e C About v C The partial derivatives of γ are:
Figure BDA0003709352020000151
based on the initial value [ v ] of the flying speed given in step S3 C ,γ x ,γ y ,γ z ] 0 And the initial values of the equality constraint and inequality constraint multipliers [ lambda, kappa ] 1 ,κ 2 ,κ 3 ] 0 A final solution to the speed of flight can be obtained by iterating a Powell-based mixed-value iterator Minpack-1 [ More j.j., garhow b.s., and hillboom k.e.user guide for Minpack-1.Argonne National lab.rept.anl-80-74, Illinois,1980 ].
After the final fly-over velocity solution is obtained in step S32, a phase difference is generated
Figure BDA0003709352020000152
Wherein u is Cf Extrapolating the phase of the spacecraft to the moment of flight for the spacecraft,
Figure BDA0003709352020000153
and
Figure BDA0003709352020000154
is a spacecraft havingThe deviation of precession rate when maneuvering is relative to that when no maneuvering, the relative average angular velocity is approximated to the first order, and the following can be obtained:
Figure BDA0003709352020000155
the phase correction is carried out by using the semimajor axis deviation, so that the delta u + delta n C Δ t is 0, and the modified semimajor axis expression is:
Figure BDA0003709352020000156
wherein N ∈ [ -1,1] represents the number of turns, and the magnitude of the flying speed after correction is as follows:
Figure BDA0003709352020000157
the resulting flying speed is
Figure BDA0003709352020000158
Step S33, when the terminal speed is obtained, the fly-over cost estimation is converted into the rendezvous cost estimation, and the estimation is carried out according to the rendezvous estimation algorithm, the rendezvous estimation algorithm can refer to the existing mature perturbation rendezvous estimation algorithm (Luoya, Huangsho Ying, Lihengzheng, Wushenggai, Zhang-in, and Populan) 2 Method for quickly estimating optimal speed increment of long-time rail crossing under perturbation [ P ]]The province of Hunan province: CN110789739B,2020-12-11.
Step S34, repeating steps S30-S33, counting the speed increment delta v generated in each iteration, and selecting the minimum value to be recorded as delta v best The corresponding fly-over velocity is solved as [ v C ,γ x ,γ y ,γ z ] best . As shown in the flow chart of fig. 2.
The invention is further described below by way of examples.
1) Establishing an inequality constraint expression of a terminal for fly-by patrol;
given alpha smin Given as 30deg, δ v max 2km/s, given r pmin ＝6600km。
2) And calculating the optimal flying time and the plane azimuth of the flying orbit according to the orbital out-of-plane degree of the spacecraft and the target.
2.1) recording the initial epoch as t 0 59274MJD shows the number of orbits of the spacecraft and target at the initial moment in table 1:
TABLE 1 spacecraft and target at t 0 59274MJD orbit number
Figure BDA0003709352020000161
According to the spherical triangle formula, the heterosurface degree phi of the two initial orbital planes is calculated to be 2.349deg, and the maximum heterosurface degree phi is calculated max ＝15.298deg。
2.2) setting the number of flying turns k to 11, and the optimal flying time delta t opt 68731.86s, the fly-by transition time Δ t is given 68731.86 s.
2.3) according to phi and phi max And Δ t opt Determining a unit normal vector h of the flying-over orbital plane e ＝[0.12908,-0.52305,0.84247] T 。
3) And giving an initial value of the flying speed according to the flying orbital plane azimuth, and establishing multi-constraint flying optimal speed estimation in fixed time based on the Cohen-Tack condition.
3.1) given θ max ＝1.5deg，e Cmax 0.02,. multidot.v Cmax ＝7.57km/s，v Cmin 7.42 km/s. Setting the value step length of theta to be 0.5deg, then N 2 7; v. the C The step length of the value of (1) is 0.01km/s, then N 1 15. Alpha and beta are both 0, and the initial value of the first flying speed is [ v C ,γ x ,γ y ,γ z ] 0 ＝[7.425,0.464,0.783,0.415]。
3.2) iterating the initial velocity value by using a Minpack-1 solver according to the established 8-dimensional equation set to obtain an iteration value [ v ] of the flying velocity C ,γ x ,γ y ,γ z ]＝[7.500,0.436,0.877,0.204]。
3.3) the phase difference produced at this speed is 96.085 DEG, and a new flying speed is obtained by trimming with a magnitude of v C ＝7.446km/s。
3.4) obtaining the flying speed, estimating the transition cost according to the orbit root difference between the new tail end flying state and the initial state of the spacecraft, and obtaining the corresponding speed increment delta v which is 215.6 m/s.
3.5) circularly taking values for 105 times, sequentially executing the steps 3.1) to 3.4), and taking the solution of the minimum delta v as the optimal flying speed [ v C ,γ x ,γ y ,γ z ] best And an optimum velocity increment Δ v best . In order to compare the accuracy and the efficiency of the algorithm, a differential evolution algorithm is adopted to carry out 3-pulse fly-by trajectory optimization on the same example, and the final optimization and estimation results and relative errors are shown in a table 2:
TABLE 2 fly-by estimation versus optimization of differential evolution algorithm
Figure BDA0003709352020000171
In the table, the relative error of 4 variables of the optimized flying speed is small, the maximum is not more than 5%, and the corresponding speed increment error is not more than 8%, so the precision of the estimation algorithm is considerable.
In optimization time, the 3-pulse flying track has 9 design variables, the population of a differential evolution algorithm is set as 108, the evolution algebra is 200, the number of parallel kernels is 6, and the average time for 20-time operation is 15 s; the algorithm provided by the invention only needs about 0.001s for each cycle, and needs about 0.1s for 105 times of operation, which is two orders of magnitude less than that of the differential algorithm.
The method sets the magnitude and three direction components of the tail end flying speed as 4 unknowns, designs an objective function of an optimization problem as a speed increment under a single pulse, introduces 3 inequality constraints including a flying inspection sunlight angle, a flying relative speed and a track near-location height during flying, lists an 8-dimensional equation set by using a Cohn-Tack condition, obtains a flying optimal speed estimation by solving the equation set, changes the flying estimation into a rendezvous estimation, and estimates the transfer cost by using a known track root number difference.
Although the embodiments of the present invention have been described with reference to the accompanying drawings, various changes or modifications may be made by the patentees within the scope of the appended claims, and within the scope of the invention, as long as they do not exceed the scope of the invention described in the claims.

Claims (4)

1.A method for estimating the optimal speed of multi-constraint flying of a spacecraft is characterized by comprising the following steps:
s1, establishing a terminal inequality constraint expression of the fly-by patrol;
s2, calculating the optimal flying time and the plane orientation of the flying orbit according to the orbital out-of-plane degree of the spacecraft and the target;
and S3, giving an initial value of the flying speed according to the flying orbital plane orientation, and establishing a multi-constraint flying optimal speed estimation in fixed time based on the Cohen-Tack condition.
2. The method for estimating the multi-constraint flying optimal speed of the spacecraft of claim 1, wherein the step S1 specifically comprises:
setting the starting time t of the spacecraft under the geocentric inertial coordinate system 0 And a position velocity at the time of flight t is [ r ] C0 ,v C0 ]And [ r C ,v C ]And the target is at the starting time t 0 And a position velocity at the time of flight t is [ r ] T0 ,v T0 ]And [ r T ,v T ]Setting the flying task as non-deviation flying with position coincidence and speed deviation; the motion of the spacecraft and the target in the set time period before the flying is set to be uniform and linear, so that the relative position vector delta r of the spacecraft and the target can be deduced to be r T –r C With relative velocity vector δ v ═ v T –v C The cross product of (a) is 0, and the sight line direction before the position coincidence is replaced by delta v; setting a sight line direction and sun position vector r S The lower limit of the included angle of the direction is alpha smin The constraint expression on the viewing sun angle is obtained as follows:
arccos((r S /r S )·(δv/δv))-(π-α smin )≤0 (1)
setting the relative speed of the spacecraft in flying not to exceed the upper limit delta v max ：
δv-δv max ≤0 (2)
Setting the height of the near point of the spacecraft in flying to be not lower than the lower limit r pmin ：
r pmin -a C (1-e C )≤0 (3)
In the formula, a C And e C The semi-major axis and eccentricity of the spacecraft at the moment of flight t.
3. The method for estimating multi-constraint flyover optimal speed of a spacecraft according to claim 1, wherein the step S2 specifically comprises the following steps:
s20, setting the spacecraft at the starting time t 0 And the number of tracks at the time of flight t is [ a ] C0 ,e C0 ,i C0 ,Ω C0 ,ω C0 ,M C0 ]And [ a ] C ,e C ,i C ,Ω C ,ω C ,M C ]And the target is at the starting time t 0 And the number of tracks at the time of flight t is [ a ] T0 ,e T0 ,i T0 ,Ω T0 ,ω T0 ,M T0 ]And [ a ] T ,e T ,i T ,Ω T ,ω T ,M T ]A, e, i, omega, M are respectively a semi-major axis, eccentricity, inclination, elevation intersection declination, near arch point latitude argument and near point angle, subscripts C0 and C are the space vehicle at t 0 And the number of tracks at time T, subscripts T0 and T targeting at T 0 And the number of tracks at time t, the track transfer duration Deltat being t-t 0 Calculating out the heterozygosity phi of the two initial orbital planes and the maximum heterozygosity phi allowed by the relative speed constraint according to a spherical triangle formula max Comprises the following steps:
Figure FDA0003709352010000011
in the formula, mu is the universal gravitation constant of the central celestial body;
s21, when phi is less than or equal to phi max The spacecraft allows flying over the target by coplanar maneuvers only, at which time the optimal flying time Δ t opt The expression of (a) is:
Figure FDA0003709352010000021
in the formula u T Phase, u, of the initial orbital plane intersection of the spacecraft and the target on the target orbit T0 An initial phase of the target, k being a non-negative integer, u T The expression of (a) is:
Figure FDA0003709352010000022
in the formula u T There are two solutions, only the solution on the side of the shadow,
Figure FDA0003709352010000023
and
Figure FDA0003709352010000024
aiming at considering only J 2 The orbit root precession rate under the long-term perturbation comprises the expression of the orbit root precession rate including a spacecraft:
Figure FDA0003709352010000025
in the formula, J 2 Is a perturbation term constant, R E Is the earth's average equatorial radius;
s22, according to phi and phi max And Δ t opt Determining a unit normal vector h of the flying-over orbital plane e Time biasThe threshold value of the difference is Δ t threshold When phi is less than or equal to phi max And | Δ t- Δ t opt |≤Δt threshold Considering that the spacecraft can perform coplanar maneuver, the angular edge (i) in the spherical triangle is used C0 -i T0 -u T ) Solving the elevation difference delta omega of the elevation intersection point as follows:
Figure FDA0003709352010000026
when phi is less than or equal to phi max And | Δ t- Δ t opt |>Δt threshold Or phi>φ max In time, the spacecraft needs to be maneuvered to orbit in different planes according to the corner angle (i) in the spherical triangle T0 -u T -φ max ) To solve for the inclination i C And the rise-crossing point declination difference delta omega is:
Figure FDA0003709352010000027
according to Δ Ω ═ Ω T0 –Ω C Thereby obtaining a unit normal vector h e Is composed of
Figure FDA0003709352010000028
4. The method for estimating multi-constraint flying optimal speed of a spacecraft according to claim 1, wherein the step S3 specifically comprises the following steps:
s30, establishing O, r with the target as the origin T The direction being the x-axis, h e Setting the magnitude of the flying speed and three directional components as v in a relative coordinate system of the LVLH with the direction as a z-axis C ,γ x ,γ y ,γ z Velocity vector v C And h e Orthogonal, magnitude of flying velocity v C In the xOy plane, the magnitude v of the flying speed is set C The included angle between the positive direction of the y axis is theta epsilon-theta max ,θ max ]，θ max For boundary values of the angle, the three directional components of the velocity can be described by θ by the coordinate system change as:
Figure FDA0003709352010000031
in the formula (I), the compound is shown in the specification,
Figure FDA0003709352010000032
a rotation matrix from an LVLH coordinate system to a geocentric inertia coordinate system;
the magnitude of the velocity is deduced according to the displacement formula of the position velocity and the number of the tracks C Expression (c):
Figure FDA0003709352010000033
in the formula, e Cmax Is the allowable upper limit of the eccentricity of the fly-over orbit;
by discrete v C And the value of theta in the respective upper and lower bounds, and the guess value is:
Figure FDA0003709352010000034
in the formula, N 1 And N 2 For the number of discrete points, α and β are both integers, and θ is guess Carrying in formula (11) to obtain a guess initial value of the direction;
s31, if the set flying trace is estimated by a single pulse, the number of orbit deviations caused by the single pulse are:
Figure FDA0003709352010000035
the magnitude of the velocity increment Δ v of a single pulse is:
Figure FDA0003709352010000036
in the formula,. DELTA.v a ＝0.5v C0 Δa/a C0 ,Δv i ＝v C0 Δi,Δv Ω ＝v C0 ΔΩsini C0 ,
Figure FDA0003709352010000037
v C0 The average orbit speed of the spacecraft at the initial moment is taken as the average orbit speed;
taking the square of the speed increment as an objective function of the optimal problem, the expression is as follows:
Figure FDA0003709352010000041
since the modulus values in the three unit directions of velocity are 1, the equality constraint is introduced as:
Figure FDA0003709352010000042
establishing a Lagrange equation, wherein the expression is as follows:
L＝f+λ T g+κ T σ (18)
in the formula, λ is an equality constraint multiplier, κ is an inequality constraint multiplier, and according to the kun-tak condition, the partial derivatives of the four unknowns are solved, so that an 8-dimensional equality equation set is obtained:
Figure FDA0003709352010000043
in the formula, relate to a C ,e C ,i C ,Ω C ,ω C For v C ,γ x ,γ y ,γ z The partial derivatives of (2) are expressed by the number of the tracks in terms of position and velocity, and are recorded in terms of h C ＝[h Cx ,h Cy ,h Cz ] T ，e C ＝[e Cx ,e Cy ,e Cz ] T The expression is:
Figure FDA0003709352010000051
wherein h is C For the angular momentum of the fly-over orbit, let γ ═ γ x ,γ y ,γ z ] T The partial derivatives can be derived as:
Figure FDA0003709352010000052
Figure FDA0003709352010000053
Figure FDA0003709352010000054
Figure FDA0003709352010000055
Figure FDA0003709352010000061
Figure FDA0003709352010000062
Figure FDA0003709352010000063
note i x ＝[1,0,0] T ,i y ＝[0,1,0] T ,i z ＝[0,0,1] T Angular momentum vector h C And its modulus value h C About v C The partial derivatives of γ are:
Figure FDA0003709352010000064
Figure FDA0003709352010000065
relative velocity vector δ v and its modulus value δ v with respect to v C The partial derivatives of γ are:
Figure FDA0003709352010000066
eccentricity vector e C About v C The partial derivatives of γ are:
Figure FDA0003709352010000071
based on the initial value [ v ] of the flying speed given in step S3 C ,γ x ,γ y ,γ z ] 0 And the initial values of the equality constraint and inequality constraint multipliers [ lambda, kappa ] 1 ,κ 2 ,κ 3 ] 0 When the flying speed is 0, a Powell-based mixed numerical value iterator iterates to obtain a final flying speed solution;
s32, obtaining the final flying speed solution, generating a phase difference
Figure FDA0003709352010000072
Wherein u is Cf Extrapolating the phase of the spacecraft to the moment of flight for the spacecraft,
Figure FDA0003709352010000073
and
Figure FDA0003709352010000074
for the precession rate deviation of the spacecraft when the spacecraft has movement relative to no movement, the relative average angular velocity is approximated to the first order, and the following can be obtained:
Figure FDA0003709352010000075
the phase correction is carried out by using the semimajor axis deviation, so that the delta u + delta n C Δ t is 0, and the modified semimajor axis expression is:
Figure FDA0003709352010000076
wherein N ∈ [ -1,1] represents the number of turns, and the magnitude of the flying speed after correction is as follows:
Figure FDA0003709352010000077
the resulting flying speed is
Figure FDA0003709352010000078
S33, converting the flying cost estimation into a rendezvous cost estimation, and estimating according to a rendezvous estimation algorithm;
s34, repeating the steps S30-S33, counting the speed increment delta v generated by each iteration, and selecting the minimum value to be recorded as delta v best The corresponding fly-over velocity is solved as [ v C ,γ x ,γ y ,γ z ] best 。
CN202210737971.9A 2022-06-23 2022-06-23 Estimation method for multi-constraint fly-through optimal speed of spacecraft Active CN114889849B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210737971.9A CN114889849B (en) 2022-06-23 2022-06-23 Estimation method for multi-constraint fly-through optimal speed of spacecraft

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210737971.9A CN114889849B (en) 2022-06-23 2022-06-23 Estimation method for multi-constraint fly-through optimal speed of spacecraft

Publications (2)

Publication Number Publication Date
CN114889849A true CN114889849A (en) 2022-08-12
CN114889849B CN114889849B (en) 2023-10-20

Family

ID=82730118

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210737971.9A Active CN114889849B (en) 2022-06-23 2022-06-23 Estimation method for multi-constraint fly-through optimal speed of spacecraft

Country Status (1)

Country Link
CN (1) CN114889849B (en)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115790607A (en) * 2023-01-31 2023-03-14 南京航空航天大学 Short arc historical data-based non-cooperative target maneuvering characteristic detection method
CN116495196A (en) * 2023-06-29 2023-07-28 南京航空航天大学 Determination method for deep space exploration spacecraft star capturing orbit

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6302354B1 (en) * 2000-06-07 2001-10-16 The Aerospace Corporation Space vehicular fly-by guidance method
CN107992682A (en) * 2017-12-05 2018-05-04 北京理工大学 A kind of optimal multiple-pulse transfer method of interplanetary multi-body system asteroid detection
CN112231943A (en) * 2020-12-17 2021-01-15 中国人民解放军国防科技大学 Multi-star fly-over sequence searching method and system containing 'one stone and multiple birds' fly-over segments
CN113343561A (en) * 2021-05-25 2021-09-03 中国科学院国家空间科学中心 Method and system for solving optimal moon fly-by transfer orbit
CN114019995A (en) * 2021-11-09 2022-02-08 北京航空航天大学 Different-plane multi-target rapid fly-by trajectory planning method with optimal fuel

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6302354B1 (en) * 2000-06-07 2001-10-16 The Aerospace Corporation Space vehicular fly-by guidance method
CN107992682A (en) * 2017-12-05 2018-05-04 北京理工大学 A kind of optimal multiple-pulse transfer method of interplanetary multi-body system asteroid detection
CN112231943A (en) * 2020-12-17 2021-01-15 中国人民解放军国防科技大学 Multi-star fly-over sequence searching method and system containing 'one stone and multiple birds' fly-over segments
CN113343561A (en) * 2021-05-25 2021-09-03 中国科学院国家空间科学中心 Method and system for solving optimal moon fly-by transfer orbit
CN114019995A (en) * 2021-11-09 2022-02-08 北京航空航天大学 Different-plane multi-target rapid fly-by trajectory planning method with optimal fuel

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
朱彦伟，杨乐平: "航天器在轨服务接近策略研究", 中国空间科学技术, no. 01, pages 14 - 20 *
潘迅，泮斌峰，唐硕: "求解中途飞越燃料最优转移轨道的同伦方法", 宇航学报, vol. 38, no. 04, pages 393 - 400 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115790607A (en) * 2023-01-31 2023-03-14 南京航空航天大学 Short arc historical data-based non-cooperative target maneuvering characteristic detection method
CN115790607B (en) * 2023-01-31 2023-05-12 南京航空航天大学 Non-cooperative target maneuvering characteristic detection method based on short arc historical data
CN116495196A (en) * 2023-06-29 2023-07-28 南京航空航天大学 Determination method for deep space exploration spacecraft star capturing orbit
CN116495196B (en) * 2023-06-29 2023-09-29 南京航空航天大学 Determination method for deep space exploration spacecraft star capturing orbit

Also Published As

Publication number Publication date
CN114889849B (en) 2023-10-20

Similar Documents

Publication Publication Date Title
CN109823571B (en) Multi-stage attitude control method for remote sensing micro-nano satellite
CN114889849A (en) Estimation method for multi-constraint flying-over optimal speed of spacecraft
CN109911249B (en) Interstellar transfer limited thrust orbit-entering iterative guidance method for low thrust-weight ratio aircraft
CN109269504B (en) Attitude maneuver path planning method with terminal constraint
CN112572835B (en) Satellite in-orbit angular momentum management and control method with attitude switching function
CN110632935B (en) Autonomous control method for formation satellite flying around
CN110244767B (en) Formation configuration reconstruction optimization method adopting finite element method
CN113602532A (en) Solid carrier rocket in-orbit correction method
CN111605737A (en) Spacecraft three-phase control multi-level collaborative planning and agile maneuvering method
CN111338367A (en) Method for determining middle track under double-pulse control of same track for freezing eccentricity ratio
CN112713922A (en) Visibility rapid forecasting algorithm of multi-beam communication satellite
CN115373423A (en) Formation capture method for commercial satellite
CN113447043A (en) GNSS-based satellite astronomical navigation system error autonomous calibration method and system
McAdams et al. MESSENGER mission design and navigation
CN110955255A (en) High-precision orbit control attitude maintaining method, system and medium based on CMG
CN114063645B (en) Eccentricity inclination angle vector-based inclination flying-around holding control effect evaluation method
CN114684389A (en) Moon-to-earth transfer window considering reentry constraint and accurate transfer orbit determination method
CN114070403B (en) Feedforward tracking control method and system for inter-satellite laser communication system
CN111891402B (en) Mars detection ground antenna pointing recovery method based on autonomous maneuvering
Lincoln et al. Components of a vision assisted constrained autonomous satellite formation flying control system
Somov et al. A space robot control at approaching and inspecting a geostationary satellite state
Han et al. Engineering Optimization Method of Orbit Transfer Strategy for All-electric Propulsion Satellites
CN112455725A (en) Method for transferring and converting pulse orbit transfer direction to limited thrust orbit
Xie et al. Autonomous guidance, navigation, and control of spacecraft
CN113071712B (en) Rapid calculation method for monthly transition injection-to-orbit strategy

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