Optimization Method and System for Electric Propulsion Orbit Transfer Strategy in Low Earth Orbit Satellite Constellation Networking

Abstract

The present invention provides a method and system for optimizing the electric propulsion orbit transfer strategy for low-Earth orbit satellite constellation networking, including: Step 1: Based on the engineering constraints of only two states of on and off according to the thrust magnitude of the electric propulsion, establish an optimization problem of the time-optimal constellation networking electric propulsion transfer orbit; Step 2: Divide the networking process into stages, and perform phased parameterization on the thrust magnitude and direction; Step 3: Determine the variables to be solved, and convert the optimization problem of the constellation networking electric propulsion transfer orbit into a shooting problem with constrained control variables; Step 4: Guess the initial value of the shooting problem, and use the Newton iteration algorithm to solve it, and finally obtain the optimal solution of the optimization problem of the constellation networking electric propulsion transfer orbit. The present invention can ensure that the transfer orbit meets the preset networking constraints, has low requirements for the attitude control system, has a simple structure, high solution efficiency, is conducive to on-board implementation, and can effectively realize the joint control of orbital elements.

Description

Technical Field

[0001] The present invention relates to the technical field of spacecraft orbit design and optimization, and in particular, to a method and system for optimizing the electric propulsion orbit transfer strategy for low-earth orbit satellite constellation networking. Background Art

[0002] In recent years, low-earth orbit constellations with high coverage, short revisit time, and even all-weather continuous coverage have been successively proposed and become the current research hotspot. The scale of such low-earth orbit constellations can reach dozens to tens of thousands of satellites, occupying a large amount of space orbital resources, and usually configured with an electric propulsion system to complete orbit control tasks such as constellation networking maneuver, configuration maintenance, and deorbiting in a fully electric propulsion mode. The constellation networking maneuver of low-earth orbit satellites refers to the transfer process in which the satellite changes its orbit and climbs from the initial orbit to the target orbit after the separation of the satellite and the rocket. Compared with general single-satellite missions, constellation networking imposes design constraints on all six orbital elements, and there are also significant differences in the design of its transfer orbits. Therefore, the design and optimization of the networking transfer orbit in the fully electric propulsion mode have become one of the key technologies for the development of low-earth orbit constellations.

[0003] Electric propulsion orbit optimization is a classic aerospace dynamics design problem, and many scholars have studied it, such as the design of asteroid exploration transfer orbits, the design of artificial sun-synchronous orbits, GTO-GEO multi-loop low-thrust transfers, and electric propulsion control for earth coverage of elliptical orbit satellites. However, there are still relatively few publicly available literatures on the design of low-earth orbit constellation networking transfer orbits. For low-earth orbit satellites, in order to reduce the probability of collision between the working satellites in the constellation, the working orbits of the satellites are often selected as near-circular orbits with small eccentricity frozen. Therefore, low-earth orbit constellation networking needs to solve the problem of joint control of multiple orbital elements to ensure that the phase difference, semi-major axis, eccentricity, and argument of perigee all meet the requirements after the orbit change and meet certain performance indicators. To address this problem, this patent considers the engineering constraints that the magnitude of the electric propulsion thrust is fixed and has only two states of on and off, constructs a time-optimal constellation networking electric propulsion transfer orbit optimization problem, proposes a thrust parameterization method, transforms the original problem into a shooting problem with constrained control variables, and solves it using a nonlinear programming algorithm. This method can achieve the joint control of the phase difference, semi-major axis, eccentricity, and argument of perigee, and is applicable to low-earth orbit constellation missions with small eccentricity frozen target orbits. Finally, the effectiveness of this method is verified using the 21st-order gravitational field model of the earth.

[0004] Currently, some research has been conducted on the spacecraft orbit transfer strategy. After retrieval, the main ones related to the design of orbit transfer strategies are as follows:

[0005] Patent "An Optimization Method for the Orbit Transfer Strategy of Geostationary Orbit Satellites" (Publication No.: CN102424116A) overcomes the deficiencies of the prior art and provides an optimization method for the orbit transfer strategy of geostationary orbit satellites, reasonably determining various constraint conditions for the design of the orbit transfer strategy to reduce manual intervention, calculation time, and computational workload during the design process of the orbit transfer strategy. The orbit transfer strategy of the geostationary orbit satellites concerned by this patent has significant differences in models and optimization methods from those of low-earth orbit satellite constellations, and its related design methods cannot be borrowed.

[0006] Patent "Calculation Method, System, and Medium for the Orbit Transfer Strategy of GEO Satellites Based on Particle Swarm Optimization Algorithm" (Publication No.: CN108216687A) also focuses on the orbit transfer of geostationary orbit satellites. Taking the ignition time and ignition direction as optimization variables and the propellant consumption as the objective function, by setting an initial particle population and performing evolutionary calculations according to the algorithm, the desired solution can be obtained more quickly, improving the calculation efficiency. The orbit transfer strategy of the geostationary orbit satellites concerned by this patent also has significant differences in models and optimization methods from those of low-earth orbit satellite constellations, and its related design methods cannot be borrowed.

[0007] Patent "Satellite Constellation Reconfiguration Method, Device, and Storage Device Based on Inclined-Plane Orbit Transfer Strategy" (Publication No.: CN107885917A) provides a satellite constellation reconfiguration method, device, and storage device based on the inclined-plane orbit transfer strategy. By means of satellite maneuvering orbit transfer, the orbital parameters of the satellites are optimized, the spatial configuration of the satellite constellation is changed, and the satellite constellation is reconfigured to meet the performance requirements of emergency missions for satellite network constellation. This patent focuses on the constellation reconfiguration method, while this patent focuses on the constellation maneuvering during the initial stage of satellite launch into orbit, and the design objects are different.

[0008] Patent "An Autonomous Orbit Transfer Method for Satellites" (Publication No.: CN101219713) overcomes the deficiencies of the prior art and proposes an autonomous orbit transfer method for on-board attitude determination, attitude maneuver, and autonomous orbit control, solving the orbit transfer problem of deep space exploration satellites and ensuring the accurate, reliable, and automatic completion of orbit transfer. This patent focuses on the autonomy of orbit transfer of deep space exploration spacecraft, which is different from the research object of this patent.

[0009] The paper "Research on the Constellation Orbit Maneuver Method under Continuous Low-Thrust Conditions" (Chinese Space Science and Technology, 2020, Issue 01: 1-16) studies the constellation orbit maneuver problem of the multi-satellite constellation mission. The combined adjustment of the phase difference and semi-major axis is achieved by using time-sharing ascending orbit maneuvers, but the combined adjustment of the semi-major axis, eccentricity, and argument of perigee is not involved, so it cannot be directly applied to the constellation mission with a target orbit of a nearly circular orbit with a small eccentricity freeze. And this patent mainly focuses on the combined adjustment of the semi-major axis, eccentricity, and argument of perigee.

[0010] Patent document CN105799954A (application number: CN201410845134.3) discloses a modular aircraft for space-based dispersed deployment of micro-nano payloads and its orbit transfer guidance method. The aircraft includes: an integrated bearing and propulsion unit, an integrated measurement and control unit, and a lightweight inter-stage adaptation and separation device; it can achieve remote autonomous and rapid orbit transfer to the target orbit. However, this patent only studies the orbit transfer of the aircraft, which is different from the research object of this patent. Summary of the Invention

[0011] Aiming at the defects in the prior art, the purpose of the present invention is to provide a method and system for optimizing the electric propulsion orbit transfer strategy for low-earth orbit satellite constellation networking.

[0012] According to the method for optimizing the electric propulsion orbit transfer strategy for low-earth orbit satellite constellation networking provided by the present invention, it includes:

[0013] Step 1: Based on the engineering constraints of only two states of on and off, establish an optimal-time constellation networking electric propulsion transfer orbit optimization problem according to the thrust magnitude of the electric propulsion;

[0014] Step 2: Divide the networking process into stages, and perform phased parameterization processing on the thrust magnitude and direction;

[0015] Step 3: Determine the variables to be solved, and convert the constellation networking electric propulsion transfer orbit optimization problem into a shooting problem with constrained control variables;

[0016] Step 4: Guess the initial value of the shooting problem, and use the Newton iteration algorithm to solve it, and finally obtain the optimal solution of the constellation networking electric propulsion transfer orbit optimization problem.

[0017] Preferably, the step 1 includes:

[0018] Step 1.1: Establish the orbital dynamics equation of the satellite in the geocentric inertial reference coordinate system:

[0019]

[0020]

[0021]

[0022] Among them, r and v are the position and velocity vectors respectively; f P is the perturbation acceleration including J2 and atmospheric drag; μ is the gravitational constant of the earth; r is the radius of the earth; g e is the gravitational acceleration at the earth's sea level; m is the mass of the spacecraft; I sp is the specific impulse of the thruster; α is the unit vector of the thrust direction; T maxis the maximum thrust; u is the ratio of the actual thrust to the maximum thrust;

[0023] Based on the engineering constraint of only two states of on and off, the value of u is 0 or 1;

[0024] The unit vector in the thrust direction is further expressed in the inertial coordinate system as:

[0025]

[0026] where is the thrust azimuth angle, and δ is the thrust elevation angle;

[0027] Let the instantaneous orbital elements of the satellite be:

[0028] k osc = [a e i Ω ω M]

[0029] where a is the semi-major axis, e is the eccentricity, i is the orbital inclination, Ω is the right ascension of the ascending node, ω is the argument of perigee, and M is the mean anomaly;

[0030] Let represent the mean orbital elements of the satellite, x = [r v] be the rectangular coordinate parameters, then k osc 、k ave and x can be converted to each other;

[0031] Step 1.2: Let the initial time be t0 and the end time be t f , then the networking constraint in the satellite plane at the initial and end times is expressed as:

[0032] r(t0) = r0, v(t0) = v0, m(t0) = m0

[0033]

[0034]

[0035] In the formula, and are constants, representing the expected semi-major axis, eccentricity, and argument of perigee at the end time, represents the along-track angle, defined as:

[0036]

[0037] At the end time, the expected is determined according to the following formula:

[0038]

[0039] In the formula, represents the initial at time t0, Indicates the rate of change of the along-track angle;

[0040] Step 1.3: Set the control variables and state variables:

[0041] u = [u α] T

[0042] y = [r v m] T

[0043] Then the orbital dynamics equation is expressed as:

[0044]

[0045] The corresponding boundary constraints are expressed as:

[0046] b ini [y(t0)] = y(t0) - y0 = 0

[0047]

[0048] where y0 = [r0 v0 m0] T

[0049] The problem P of the electric propulsion transfer orbit in the networking process is to find the optimal control variables:

[0050]

[0051] while satisfying:

[0052]

[0053] b ini [y(t0)] = 0, b ter [y(t f )] = 0

[0054] u = 0 or u = 1, α T α = 1.

[0055] Preferably, the said step 2 includes:

[0056] Step 2.1: According to the performance index requirements of the fastest orbit injection and the motion characteristics of the networking orbit, divide u and α in the time domain at times t1, t2 and t3. On the basis of full throttle, divide the transfer orbit into a cruise stage, a rising orbit stage and an eccentricity vector adjustment stage in chronological order;

[0057] Step 2.2: Set the thrust magnitude and direction in each stage as:

[0058]

[0059]

[0060] In the formula, α a and α eω are respectively the unit vectors of the thrust direction in the ascending orbit stage and the eccentricity vector adjustment stage;

[0061] Step 2.3: To meet the performance index requirements of the fastest orbit injection, in the ascending orbit stage and the eccentricity vector adjustment stage, the thrust directions are respectively set along the tangential direction and the direction perpendicular to the apsidal line in the orbital plane, which is expressed as:

[0062]

[0063] α eω = (cosσ)i p + (sinσ)j p

[0064] where σ represents the angle between the thrust direction and the apsidal line, and:

[0065]

[0066]

[0067] where, i p and j p are calculated based on the average orbital parameters to ensure a smooth change in the thrust direction.

[0068] Preferably, the said step 3 includes:

[0069] Step 3.1: For the given initial state t0 and y0, the variables to be solved are defined as:

[0070] z = [Δt1 Δt2 Δt3 σ] T

[0071] where, Δt1 = t1 - t0, Δt2 = t2 - t1, Δt3 = t3 - t2;

[0072] Step 3.2: Convert the problem P into the following shooting problem:

[0073]

[0074] Preferably, the said step 4 includes:

[0075] Step 4.1: Make guesses for Δt1, Δt2, Δt3 and σ to determine the initial values of the iterative variables;

[0076] Step 4.2: Use the Runge - Kutta method to calculate the cost function F(z), and use the Newton iteration algorithm for iterative calculation until ||F(z)|| converges to less than the specified threshold.

[0077] The low-Earth orbit satellite constellation networking electric propulsion orbit transfer strategy optimization system provided by the present invention includes:

[0078] Module M1: Based on the engineering constraints of only two states of on and off for the thrust of the electric propulsion, establish an optimization problem of the time-optimal constellation networking electric propulsion transfer orbit.

[0079] Module M2: Divide the networking process into stages, and parameterize the thrust magnitude and direction in stages.

[0080] Module M3: Determine the variables to be solved, and convert the optimization problem of the constellation networking electric propulsion transfer orbit into a shooting problem with constrained control variables.

[0081] Module M4: Guess the initial value of the shooting problem, and use the Newton iteration algorithm to solve it, and finally obtain the optimal solution of the optimization problem of the constellation networking electric propulsion transfer orbit.

[0082] Preferably, the Module M1 includes:

[0083] Module M1.1: Establish the orbital dynamics equation of the satellite in the geocentric inertial reference coordinate system:

[0084]

[0085]

[0086]

[0087] where r and v are the position and velocity vectors respectively; f P is the perturbation acceleration including J2 and atmospheric drag; μ is the Earth's gravitational constant; r is the Earth's radius; g e is the Earth's sea-level gravitational acceleration; m is the mass of the spacecraft; I sp is the specific impulse of the thruster; α is the unit vector in the thrust direction; T max is the maximum thrust; u is the ratio of the actual thrust to the maximum thrust;

[0088] Based on the engineering constraints of only two states of on and off, the value of u is 0 or 1;

[0089] The unit vector in the thrust direction is further expressed in the inertial coordinate system as:

[0090]

[0091] where, is the thrust azimuth angle, and δ is the thrust elevation angle;

[0092] Let the instantaneous orbital elements of the satellite be:

[0093] k osc = [a e i Ω ω M]

[0094] where a is the semi-major axis, e is the eccentricity, i is the inclination, Ω is the right ascension of the ascending node, ω is the argument of perigee, and M is the mean anomaly;

[0095] Let represent the mean orbital elements of the satellite, x = [r v] be the rectangular coordinate parameters, then k osc and k ave and x can be converted to each other;

[0096] Module M1.2: Let the initial time be t0 and the end time be t f , then the formation flying constraints in the plane of the satellite at the initial and end times are expressed as:

[0097] r(t0) = r0, v(t0) = v0, m(t0) = m0

[0098]

[0099]

[0100] In the formula, and are constants, representing the desired semi-major axis, eccentricity, and argument of perigee at the end time, represents the along-track angle, defined as:

[0101]

[0102] At the end time, the desired is determined according to the following formula:

[0103]

[0104] In the formula, represents the initial value at t0, represents the rate of change of the along-track angle;

[0105] Module M1.3: Let the control variables and state variables:

[0106] u = [u α] T

[0107] y = [r v m] T

[0108] Then the orbital dynamics equation is expressed as:

[0109]

[0110] The corresponding boundary constraints are expressed as:

[0111] b ini [y(t0)] = y(t0) - y0 = 0

[0112]

[0113] where y0 = [r0 v0 m0] T

[0114] In the process of network formation, for the electric propulsion transfer orbit problem P, the optimal control variables are to be found:

[0115]

[0116] while satisfying:

[0117]

[0118] b ini [y(t0)] = 0, b ter [y(t f )] = 0

[0119] u = 0 or u = 1, α T α = 1.

[0120] Preferably, the module M2 includes:

[0121] Module M2.1: According to the performance index requirements of the fastest orbit injection and the motion characteristics of the network formation orbit, u and α are divided in the time domain at times t1, t2, and t3. On the basis of full throttle, the transfer orbit is divided into a cruise stage, a lift stage, and an eccentricity vector adjustment stage in chronological order;

[0122] Module M2.2: Set the magnitude and direction of the thrust in each stage as:

[0123]

[0124]

[0125] where α a and α eω are the unit vectors of the thrust direction in the lift stage and the eccentricity vector adjustment stage respectively;

[0126] Module M2.3: To meet the performance index requirements of the fastest orbit injection, in the lift stage and the eccentricity vector adjustment stage, set the thrust direction along the tangential direction and the direction perpendicular to the apsidal line in the orbital plane respectively, expressed as:

[0127]

[0128] αeω = (cosσ)i p + (sinσ)j p

[0129] where σ represents the angle between the thrust direction and the arch line, and:

[0130]

[0131]

[0132] where i p and j p are calculated based on the average orbital parameters to ensure a smooth change in the thrust direction.

[0133] Preferably, the module M3 includes:

[0134] Module M3.1: For a given initial state t0 and y0, define the variables to be solved as:

[0135] z = [Δt1 Δt2 Δt3 σ] T

[0136] where Δt1 = t1 - t0, Δt2 = t2 - t1, Δt3 = t3 - t2;

[0137] Module M3.2: Convert the problem P into the following shooting problem:

[0138]

[0139] Preferably, the module M4 includes:

[0140] Module M4.1: Make guesses for Δt1, Δt2, Δt3, and σ to determine the initial values of the iterative variables;

[0141] Module M4.2: Use the Runge - Kutta method to calculate the cost function F(z), and use the Newton iteration algorithm for iterative calculation until ||F(z)|| converges to less than the specified threshold.

[0142] Compared with the prior art, the present invention has the following beneficial effects:

[0143] The method in the present invention considers that the magnitude of the electric propulsion thrust is fixed, with only two states of on and off. According to the requirement of the fastest orbit injection, it makes full use of the motion characteristics of the networking orbit, and proposes a thrust parameterization method, which can ensure that the transfer orbit meets the pre - set networking constraints, has low requirements for the attitude control system, and the thrust parameters are only 4, with a simple structure, high solution efficiency, is conducive to on - satellite implementation, and can effectively achieve the combined control of orbital elements. BRIEF DESCRIPTION OF THE DRAWINGS

[0144] Other features, objects, and advantages of the present invention will become more apparent from the following detailed description of non - limiting embodiments read in conjunction with the accompanying drawings:

[0145] Figure 1 This is the flow chart of the method of the present invention;

[0146] Figure 2 This is the thrust profile in the specific embodiment;

[0147] Figure 3 This is the semi - major axis and along - track angle deviation in the specific embodiment;

[0148] Figure 4 This is the eccentricity vector and thrust magnitude in the specific embodiment;

[0149] Figure 5 This is the thrust direction angle and fuel consumption in the inertial coordinate system in the specific embodiment. Specific Embodiment

[0150] The present invention will be described in detail below in conjunction with specific embodiments. The following embodiments will help those skilled in the art to further understand the present invention, but do not limit the present invention in any form. It should be noted that for those of ordinary skill in the art, without departing from the concept of the present invention, several changes and improvements can still be made. These all belong to the protection scope of the present invention.

[0151] Embodiment:

[0152] The present invention provides an optimization algorithm for the electric - propulsion orbit - transfer strategy of a low - Earth - orbit satellite constellation networking. Considering that the magnitude of the electric - propulsion thrust is fixed and there are only two states: on and off, according to the requirement of the fastest orbit injection, taking full advantage of the motion characteristics of the networking orbit, a thrust parameterization method is proposed, which can ensure that the transfer orbit meets the pre - set networking constraints, has low requirements for the attitude - control system, and the thrust parameters are only 4, with a simple structure and high solution efficiency, which is conducive to on - satellite implementation and can effectively achieve the combined control of orbital elements.

[0153] As Figure 1 , the present invention provides an optimization algorithm for the electric - propulsion orbit - transfer strategy of a low - Earth - orbit satellite constellation networking, including the following steps:

[0154] Step 1: Considering the engineering constraints that the magnitude of the electric - propulsion thrust is fixed and there are only two states: on and off, establish an optimization problem for the time - optimal constellation - networking electric - propulsion transfer orbit.

[0155] Among them, the specific process includes:

[0156] Step 101: In the geocentric inertial reference coordinate system, establish the orbital dynamics equation of the satellite as:

[0157]

[0158]

[0159] where r and v are the position and velocity vectors respectively, and f P is the perturbation acceleration including J2, atmospheric drag, etc., μ is the gravitational constant of the Earth, r is the radius of the Earth, and g e is the gravitational acceleration at the sea level of the Earth, m is the mass of the spacecraft, and I sp is the specific impulse of the thruster, α is the unit vector in the thrust direction, and T max is the maximum thrust, and u is the ratio of the actual thrust to the maximum thrust. In this paper, it is considered that the engine thrust cannot be adjusted continuously and there are only two states of on and off. Therefore, the value of u is 0 or 1. The unit vector in the thrust direction can be further expressed in the inertial coordinate system as:

[0160]

[0161] where is the thrust azimuth angle and δ is the thrust elevation angle.

[0162] Let the instantaneous orbital elements of the satellite be:

[0163] k osc =[a e i Ω ω M]

[0164] where a is the semi-major axis, e is the eccentricity, i is the orbital inclination, Ω is the right ascension of the ascending node, ω is the argument of perigee, and M is the mean anomaly. Let represent the mean orbital elements of the satellite, x = [r v] be the rectangular coordinate parameters, then k osc , k ave and x can be converted to each other.

[0165] Step 102: Let the initial time be t0 and the end time be t f , then the formation flying constraints in the plane of the satellite at the initial and end times can be expressed as:

[0166] r(t0) = r0, v(t0) = v0, m(t0) = m0

[0167]

[0168]

[0169] In the formula, and are constants, representing the desired semi-major axis, eccentricity and argument of perigee at the end time, represents the along-track angle, which is defined as:

[0170]

[0171] At the terminal moment, the desired can be determined according to the following formula:

[0172]

[0173] In the formula, represents the initial value at time t0, represents the rate of change of the track angle.

[0174] Step 103: For the convenience of description, let the control variables and state variables be:

[0175] u = [u α] T

[0176] y = [r v m] T

[0177] Then the orbital dynamics equation can be expressed as:

[0178]

[0179] The corresponding boundary constraints are expressed as:

[0180] b ini [y(t0)] = y(t0) - y0 = 0

[0181]

[0182] where y0 = [r0 v0 m0] T .

[0183] The problem P of the electric propulsion transfer orbit during the networking process can be described as follows: find the optimal control variables;

[0184]

[0185] while satisfying:

[0186]

[0187] b ini [y(t0)] = 0, b ter [y(t f )] = 0

[0188] u = 0 or 1, α T α = 1

[0189] Step 2: Divide the networking process into stages, and parameterize the thrust magnitude and direction in stages.

[0190] Among them, Step 2 specifically includes:

[0191] Step 201: According to the performance index requirements of the fastest orbit injection, fully considering the motion characteristics of the networking orbit, divide u and α in the time domain at times t1, t2, and t3. On the basis that the thrust is basically fully open, divide the transfer orbit into a cruise phase, a raising orbit phase, and an eccentricity vector adjustment phase in chronological order, as Figure 2 .

[0192] Step 202: Set the magnitude and direction of the thrust in each phase as:

[0193]

[0194]

[0195] where α a and α eω are the unit vectors of the thrust direction in the raising orbit phase and the eccentricity vector adjustment phase, respectively.

[0196] Step 203: To meet the performance index requirements of the fastest orbit injection, in the raising orbit phase and the eccentricity vector adjustment phase, set the thrust direction along the tangential direction and the direction perpendicular to the apsidal line in the orbital plane, respectively, specifically expressed as:

[0197]

[0198] α eω =(cosσ)i p +(sinσ)j p

[0199] where σ represents the angle between the thrust direction and the apsidal line, and:

[0200]

[0201]

[0202] It should be noted that to ensure the smooth change of the thrust direction, the calculations of i p and j p are based on the average orbital parameters.

[0203] Step 3: Determine the variables to be solved and establish a shooting problem with control variables constrained.

[0204] The specific process of this step is as follows:

[0205] Step 301: After the transformation in Step 2, for the given initial state t0 and y0, define the variables to be solved as

[0206] z = [Δt1 Δt2 Δt3 σ] T

[0207] where, Δt1 = t1 - t0, Δt2 = t2 - t1, Δt3 = t3 - t2.

[0208] Step 302: Convert problem P into the following shooting problem:

[0209]

[0210] It can be seen that in the above equations, the number of unknown variables is equal to the number of constraint equations, so the solution of the equations is completely determined.

[0211] Step 4: Guess the initial value of the shooting problem and solve it using the Newton iteration algorithm.

[0212] The specific process of this step is as follows:

[0213] Step 401: Guess Δt1, Δt2, Δt3 and σ to determine the initial value of the iteration variable;

[0214] Step 402: Calculate the cost function F(z) using the Runge-Kutta numerical calculation algorithm and perform iterative calculation using the Newton iteration algorithm until ||F(z)|| converges to less than the specified threshold.

[0215] The following is a numerical simulation verification of an optimization algorithm for the electric propulsion orbit transfer strategy of a low-Earth orbit satellite constellation networking.

[0216] Assume that the satellite is launched into orbit at a relatively low altitude. After the satellites are separated from each other, the initial orbital altitude of the satellite is 1000 km, and the target orbital altitude is 1200 km. Assume that the initial mass of the satellite is 170 kg, the thrust magnitude is 20 mN, and the specific impulse is 1600 s. In the table, the argument of perigee and eccentricity of the target orbit are selected as 90 deg and 0.001 to achieve the freezing of the target orbit. Considering that the inclination of the target orbit is close to 90 deg and its orbit freezing is greatly affected by the Earth's high-order gravitational field, in this simulation, a 21st-order gravitational field model of the Earth will be considered. Since the initial orbital eccentricity is close to 0 and the argument of perigee cannot be defined, two components of the eccentricity vector are used to describe it in the subsequent analysis, that is:

[0217] e x = ecosω, e y = esinω

[0218] In the target orbit, the target values of e x and e y are 0 and 0.001. In the simulation calculation results and analysis, unless otherwise specified, all orbital elements refer to the mean orbital elements.

[0219] According to the above calculation conditions, numerical simulation calculations are carried out, and the results are as shown in Figure 3 、 Figure 4, Figure 5 as shown

[0220] The low-Earth orbit satellite constellation networking electric propulsion orbit transfer strategy optimization system provided by the present invention includes: Module M1: Based on the engineering constraints of only two states of on and off for the thrust of the electric propulsion, establish an optimization problem of the time-optimal constellation networking electric propulsion transfer orbit; Module M2: Divide the networking process into stages, and perform staged parameterization on the thrust magnitude and direction; Module M3: Determine the variables to be solved, and convert the optimization problem of the constellation networking electric propulsion transfer orbit into a shooting problem with constrained control variables; Module M4: Guess the initial value of the shooting problem, and use the Newton iteration algorithm to solve it, and finally obtain the optimal solution of the optimization problem of the constellation networking electric propulsion transfer orbit.

[0221] The said Module M1 includes: Module M1.1: Establish the orbital dynamics equation of the satellite in the geocentric inertial reference coordinate system:

[0222]

[0223]

[0224]

[0225] where r and v are the position and velocity vectors respectively; f P is the perturbation acceleration including J2 and atmospheric drag; μ is the Earth's gravitational constant; r is the Earth's radius; g e is the Earth's sea-level gravitational acceleration; m is the mass of the spacecraft; I sp is the specific impulse of the thruster; α is the unit vector of the thrust direction; T max is the maximum thrust; u is the ratio of the actual thrust to the maximum thrust; based on the engineering constraints of only two states of on and off, the value of u is 0 or 1; the unit vector of the thrust direction is further expressed in the inertial coordinate system as: where is the thrust azimuth angle, δ is the thrust elevation angle; let the instantaneous orbital elements of the satellite be: k osc =[a e i Ω ω M], where a is the semi-major axis, e is the eccentricity, i is the orbital inclination, Ω is the right ascension of the ascending node, ω is the argument of perigee, M is the mean anomaly; let represent the mean orbital elements of the satellite, x = [r v] is the rectangular coordinate parameter, then k osc , k ave and x can be converted to each other;

[0226] Module M1.2: Let the initial time be t0 and the end time be t f, the in-plane networking constraints of the satellite at the initial and final times are expressed as: r(t0) = r0, v(t0) = v0, m(t0) = m0

[0227]

[0228]

[0229] where, and are constants, representing the desired semi-major axis, eccentricity, and argument of perigee at the end time, represents the along-track angle, defined as:

[0230] At the end time, the desired is determined according to the following formula:

[0231] where, represents the initial value at time t0, represents the rate of change of the along-track angle;

[0232] Module M1.3: Let the control variables and state variables be: u = [u α] T , y = [r v m] T

[0233] Then the orbital dynamics equation is expressed as:

[0234] The corresponding boundary constraints are expressed as: b ini [y(t0)] = y(t0) - y0 = 0

[0235]

[0236] where, y0 = [r0 v0 m0] T

[0237] The problem P of the electric propulsion transfer orbit during the networking process is to find the optimal control variables:

[0238] while satisfying: b ini [y(t0)] = 0, b ter [y(t f )] = 0, u = 0 or u = 1, α T α = 1.

[0239] The module M2 includes: Module M2.1: According to the performance index requirements of the fastest orbit injection and the motion characteristics of the networking orbit, u and α are divided in the time domain at times t1, t2, and t3. On the basis of full throttle, the transfer orbit is divided into a cruise phase, an orbit-raising phase, and an eccentricity vector adjustment phase in chronological order; Module M2.2: Set the thrust magnitude and direction in each phase as: where α a and α eω are the unit vectors of the thrust direction in the orbit-raising phase and the eccentricity vector adjustment phase, respectively; Module M2.3: To meet the performance index requirements of the fastest orbit injection, in the orbit-raising phase and the eccentricity vector adjustment phase, set the thrust directions along the tangential direction and the direction perpendicular to the apsidal line in the orbital plane, respectively, expressed as: α eω =(cosσ)i p +(sinσ)j p , where σ represents the angle between the thrust direction and the apsidal line, and: where i p and j p are calculated based on the average orbital parameters to ensure a smooth change in the thrust direction.

[0240] The module M3 includes: Module M3.1: For the given initial states t0 and y0, define the variables to be solved as: z = [Δt1 Δt2 Δt3 σ] T , where Δt1 = t1 - t0, Δt2 = t2 - t1, Δt3 = t3 - t2;

[0241] Module M3.2: Convert the problem P into the following shooting problem:

[0242]

[0243] The module M4 includes: Module M4.1: Make guesses for Δt1, Δt2, Δt3, and σ to determine the initial values of the iteration variables; Module M4.2: Use the Runge-Kutta method to calculate the cost function F(z) and use the Newton iteration algorithm for iterative calculation until ||F(z)|| converges to less than the specified threshold.

[0244] Those skilled in the art know that, in addition to implementing the systems, devices, and their respective modules provided by the present invention in the form of pure computer-readable program code, it is entirely possible to logically program the method steps to enable the systems, devices, and their respective modules provided by the present invention to be implemented in the form of logic gates, switches, application-specific integrated circuits, programmable logic controllers, and embedded microcontrollers, etc. to achieve the same program. Therefore, the systems, devices, and their respective modules provided by the present invention can be considered as a kind of hardware component, and the modules included therein for implementing various programs can also be regarded as the structures within the hardware component; the modules for implementing various functions can also be regarded as either software programs for implementing the method or the structures within the hardware component.

[0245] The specific embodiments of the present invention have been described above. It should be understood that the present invention is not limited to the above specific embodiments, and those skilled in the art can make various changes or modifications within the scope of the claims, which does not affect the essence of the present invention. Without conflict, the embodiments of the present application and the features in the embodiments can be combined with each other arbitrarily.

Claims

1. An optimization method for the electric propulsion orbit transfer strategy of a low-earth orbit satellite constellation networking, characterized in that, including: Step 1: Based on the engineering constraints of only two states of on and off, establish an optimization problem of the time-optimal constellation networking electric propulsion transfer orbit according to the thrust magnitude of the electric propulsion; Step 2: Divide the networking process into stages, and parameterize the thrust magnitude and direction in stages; Step 3: Determine the variables to be solved, and convert the optimization problem of the constellation networking electric propulsion transfer orbit into a shooting problem with constrained control variables; Step 4: Guess the initial values of the shooting problem, and use the Newton iteration algorithm to solve it. Finally, obtain the optimal solution of the optimization problem of the constellation networking electric propulsion transfer orbit; The said Step 1 includes: Step 1.1: Establish the orbital dynamics equation of the satellite in the geocentric inertial reference coordinate system: where r and v are the position and velocity vectors respectively; f P is the perturbation acceleration including J2 and atmospheric drag; μ is the Earth's gravitational constant; r is the Earth's radius; g e is the Earth's sea-level gravitational acceleration; m is the mass of the spacecraft; I sp is the specific impulse of the thruster; α is the unit vector in the thrust direction; T max is the maximum thrust; u is the ratio of the actual thrust to the maximum thrust; Based on the engineering constraints of only two states of on and off, the value of u is 0 or 1; The unit vector of the thrust direction is further expressed in the inertial coordinate system as: Among them, is the thrust azimuth angle, and δ is the thrust elevation angle; Let the instantaneous orbital elements of the satellite be: k osc = [a e i Ω ω M] where, a is the semi-major axis, e is the eccentricity, i is the orbital inclination, Ω is the right ascension of the ascending node, ω is the argument of perigee, and M is the mean anomaly; Let represent the mean orbital elements of a satellite, and let \(x = [r\ v]\) be the rectangular coordinate parameters. Then \(k osc , \(k ave and \(x\) can be converted to each other; Step 1.2: Let the initial time be t0 and the end time be t f , then the networking constraints in the satellite plane at the initial and end times are expressed as: r(t0) = r0, v(t0) = v0, m(t0) = m0 wherein, and are constant values, representing the semi-major axis, eccentricity, and argument of perigee expected at the end time, represents the along-track angle, defined as: At the end moment, the desired is determined according to the following formula: In the formula, represents the initial value at time t0, represents the rate of change of the along-track angle; Step 1.3: Let the control variables and state variables: u = [uα] T y = [r v m] T Then the orbital dynamics equation is expressed as: The corresponding boundary constraints are expressed as: b ini [y(t0)] = y(t0) - y0 = 0 where y0 = [r0 v0 m0] T The electric propulsion transfer orbit problem P in the networking process is to find the optimal control variables: while satisfying: b ini [y(t0)] = 0, b ter [y(t f )] = 0 u = 0 or u = 1, α T α = 1.

2. The method for optimizing the electric propulsion orbit transfer strategy of the low-Earth orbit satellite constellation networking according to claim 1, wherein The said Step 2 includes: Step 2.1: According to the performance index requirements of the fastest orbit injection and the motion characteristics of the networking orbit, divide u and α in the time domain at times t1, t2, and t3. On the basis of full thrust, divide the transfer orbit into a cruise stage, a rising orbit stage, and an eccentricity vector adjustment stage in chronological order; Step 2.2: Set the thrust magnitude and direction of each stage as: where α a and α eω are the unit vectors of the thrust directions in the ascending orbit phase and the eccentricity vector adjustment phase, respectively; Step 2.3: To meet the performance index requirements of the fastest orbit injection, in the rising orbit stage and the eccentricity vector adjustment stage, set the thrust direction along the tangential direction and the direction perpendicular to the apsidal line in the orbital plane respectively, expressed as: α eω = (cosσ)i p + (sinσ)j p where, σ represents the angle between the thrust direction and the apsidal line, and: Among them, i p and j p are calculated based on the average orbital parameters to ensure a smooth change in the thrust direction.

3. The method for optimizing the electric propulsion orbit transfer strategy of the low-earth orbit satellite constellation networking according to claim 2, wherein, The said Step 3 includes: Step 3.1: For the given initial state t0 and y0, define the variables to be solved as: z = [Δt1 Δt2 Δt3 σ] T where, Δt1 = t1 - t0, Δt2 = t2 - t1, Δt3 = t3 - t2; Step 3.2: Convert the problem P into the following shooting problem:

4. The low-orbit satellite constellation networking electric propulsion orbit transfer strategy optimization method according to claim 3, characterized in that The said Step 4 includes: Step 4.1: Guess Δt1, Δt2, Δt3, and σ, and determine the initial values of the iterative variables; Step 4.2: Use the Runge-Kutta method to calculate the cost function F(z), and use the Newton iteration algorithm for iterative calculation until ||F(z)|| converges to less than the specified threshold.

5. A system for optimizing the electric propulsion orbit transfer strategy of a low-Earth orbit satellite constellation networking, characterized in that, including: Module M1: Based on the thrust magnitude of the electric propulsion, establish an optimization problem of the time-optimal constellation networking electric propulsion transfer orbit under the engineering constraints of only two states of on and off; Module M2: Divide the networking process into stages, and parameterize the thrust magnitude and direction in stages; Module M3: Determine the variables to be solved, and convert the optimization problem of the constellation networking electric propulsion transfer orbit into a shooting problem with constrained control variables; Module M4: Guess the initial values of the target shooting problem, solve it using the Newton iteration algorithm, and finally obtain the optimal solution to the constellation networking electric propulsion transfer orbit optimization problem; The said module M1 includes: Module M1.1: Establish the orbital dynamics equation of the satellite in the geocentric inertial reference coordinate system: where r and v are the position and velocity vectors, respectively; f P is the perturbation acceleration including J2 and atmospheric drag; μ is the gravitational constant of the Earth; r is the radius of the Earth; g e is the gravitational acceleration at the sea level of the Earth; m is the mass of the spacecraft; I sp is the specific impulse of the thruster; α is the unit vector in the thrust direction; T max is the maximum thrust; u is the ratio of the actual thrust to the maximum thrust; Based on the engineering constraints with only two states of on and off, the value of u is 0 or 1; The unit vector in the thrust direction is further expressed in the inertial coordinate system as: wherein, is the thrust azimuth angle, and δ is the thrust elevation angle; Let the instantaneous orbital elements of the satellite be: k osc = [a e i Ω ω M] where, a is the semi-major axis, e is the eccentricity, i is the orbital inclination, Ω is the right ascension of the ascending node, ω is the argument of perigee, and M is the mean anomaly; Let represent the mean orbital elements of the satellite, and let \(x = [r\ v]\) be the rectangular coordinate parameters. Then \(k\) osc , \(k\) ave and \(x\) can be converted to each other; Module M1.2: Let the initial time be t0 and the end time be t f , then the networking constraints in the satellite plane at the initial and end times are expressed as: r(t0) = r0, v(t0) = v0, m(t0) = m0 In the formula, and are constant values, representing the semi-major axis, eccentricity, and argument of perigee expected at the end time, represents the along-track angle, defined as: At the end moment, the desired is determined according to the following formula: wherein, represents the initial value at time t0, represents the rate of change of the along-track angle; Module M1.3: Let the control variables and state variables: u = [u α] T y = [r v m] T Then the orbital dynamics equation is expressed as: The corresponding boundary constraints are expressed as: b ini [y(t0)] = y(t0) - y0 = 0 where y0 = [r0 v0 m0] T The electric propulsion transfer orbit problem P during the networking process is to find the optimal control variables: While satisfying: b ini [y(t0)] = 0, b ter [y(t f )] = 0 u = 0 or u = 1, α T α = 1.

6. The low-earth orbit satellite constellation networking electric propulsion orbit transfer strategy optimization system according to claim 5, characterized in that The said module M2 includes: Module M2.1: According to the performance index requirements of the fastest orbit injection and the motion characteristics of the networking orbit, divide u and α in the time domain at times t1, t2, and t3. On the basis of full thrust, divide the transfer orbit into a cruise stage, an orbit-raising stage, and an eccentricity vector adjustment stage in chronological order; Module M2.2: Set the thrust magnitude and direction in each stage as: where α a and α eω are the unit vectors of the thrust directions in the ascending orbit phase and the eccentricity vector adjustment phase, respectively; Module M2.3: To meet the performance index requirements of the fastest orbit injection, in the orbit-raising stage and the eccentricity vector adjustment stage, set the thrust direction along the tangential direction and the direction perpendicular to the apsidal line in the orbital plane respectively, expressed as: α eω = (cosσ)i p + (sinσ)j p where, σ represents the angle between the thrust direction and the apsidal line, and: where i p and j p are calculated based on the average orbital parameters to ensure a smooth change in the thrust direction.

7. The low-earth orbit satellite constellation networking electric propulsion orbit transfer strategy optimization system according to claim 6, wherein The said module M3 includes: Module M3.1: For the given initial state t0 and y0, define the variables to be solved as: z = [Δt1 Δt2 Δt3 σ] T where, Δt1 = t1 - t0, Δt2 = t2 - t1, Δt3 = t3 - t2; Module M3.2: Convert the problem P into the following target shooting problem:

8. The low-earth orbit satellite constellation networking electric propulsion orbit transfer strategy optimization system according to claim 7, characterized in that, The said module M4 includes: Module M4.1: Guess Δt1, Δt2, Δt3, and σ to determine the initial values of the iteration variables; Module M4.2: Use the Runge-Kutta method to calculate the cost function F(z), and perform iterative calculations using the Newton iteration algorithm until ||F(z)|| converges to less than the specified threshold.

Images