LOW-ORBIT SATELLITE DEORBIT CONTROL METHOD AND SYSTEM BASED ON PARTICLE SWARM ALGORITHM

Abstract

The disclosure describes a low-orbit satellite deorbit control method and system based on a particle swarm algorithm. The method includes: constructing an objective function according to a position parameter and a perturbation acceleration of each of particles in a particle swarm as well as a preset Lagrangian multiplier and a penalty factor; determining an objective fitness of each particle and a population fitness of the particle swarm based on the objective function, updating the position parameter and the velocity of each particle, and obtaining a population optimal fitness at a maximum number of iterations; updating the Lagrangian multiplier and the penalty factor according to the population optimal fitness, comparing the objective deviation value with a preset convergence condition, and determining an objective population fitness and an objective position parameter according to a comparison result to obtain an objective deorbit mode.

Claims

1. A low-orbit satellite deorbit control method based on a particle swarm algorithm, comprising: (S1) constructing an objective function according to a position parameter and a perturbation acceleration of each of particles in a particle swarm as well as a preset Lagrangian multiplier and a preset penalty factor; (S2) calculating an objective function value of each particle in the particle swarm based on the objective function, determining an objective fitness of each particle and a population fitness of the particle swarm according to a calculation result, updating the position parameter and velocity of each particle, and obtaining a population optimal fitness at a maximum number of iterations and the position parameters of particles corresponding to the population optimal fitness; and (S3) updating the Lagrangian multiplier and the penalty factor according to the population optimal fitness, obtaining a current number of iterations and an objective deviation value based on the population optimal fitness, comparing the current number of iterations and the objective deviation value with a preset iteration condition and a preset convergence condition, and determining an objective population fitness and an objective position parameter according to a comparison result to obtain an objective deorbit mode.

2. The method according to claim 1, wherein the step of constructing an objective function according to a position parameter and a perturbation acceleration of each of particles in a particle swarm as well as a preset Lagrangian multiplier and a preset penalty factor further comprises: (S11) obtaining a rate of change of each orbital element over time under the thrust applied by a thruster, and determining a control rate containing current costate variables according to the rate of change; (S12) calculating a thrust vector, a fuel consumption rate, and a perturbation acceleration of the thrust applied by the thruster; (S13) determining a particle position vector according to the current costate variables, and using the particle position vector as the position parameter of each of the particles in the particle swarm; (S14) performing orbit integration of the control rate based on the position parameter and the perturbation acceleration, and constructing the objective function according to a integration result, the preset Lagrangian multiplier and the preset penalty factor.

3. The method according to claim 2, wherein the steps of comparing the current number of iterations and the objective deviation value with a preset iteration condition and a preset convergence condition respectively and determining an objective population fitness and an objective position parameter according to a comparison result to obtain an objective deorbit mode further comprises: if the current number of iterations meets the preset iteration condition and the objective deviation value meets the preset convergence condition, determining the population optimal fitness as the objective population fitness, and determining the position parameter corresponding to the objective population fitness as the objective position parameter; if the current number of iterations does not meet the preset iteration condition or the objective deviation value does not meet the preset convergence condition, updating the objective function based on the updated Lagrangian multiplier and the updated penalty factor, and continuing to perform step S14 and subsequent steps, until the current number of iterations meets the preset iteration condition and the objective deviation value meets the preset convergence condition.

4. The method according to claim 2, wherein in step S11, the step of determining a control rate containing current costate variables according to the rate of change further comprises: based on the rate of change and the current costate variables, obtaining a Hamiltonian function by equation (7): $\begin{array}{cc}H={\lambda }_a\frac{da}{dt}+{\lambda }_e\frac{de}{dt}+{\lambda }_i\frac{di}{dt};& \left(7\right)\end{array}$ where H is a Hamiltonian function value, λ.sub.a, λ.sub.e and λ.sub.i are the current costate variables, a, e and i are semi-major axis, eccentricity and orbital inclination of an operational orbit respectively, and t is time; deriving a preset angle corresponding to the thrust applied by the thruster according to the Hamilton function (7) to obtain the control rate containing the current costate variables.

5. The method according to claim 2, wherein in step S12, the thrust vector is obtained by equation (11), the fuel consumption rate is obtained by equation (12), and the perturbation acceleration is obtained by equation (14); $\begin{array}{cc}\left|{\stackrel{.fwdarw.}{F}}_{\mathrm{thrust}}\right|=\frac{2\eta P}{{\mathrm{mgI}}_{sp}};& \left(11\right)\\ \frac{dm}{dt}=\frac{-2\eta P}{{\left({\mathrm{gI}}_{sp}\right)}^2};& \left(12\right)\\ \left[\begin{array}{c}{a}_t\\ a_n\\ a_h\end{array}\right]=\left[\begin{array}{c}\frac{\mathrm{pv}}{h\left(1+e^2+2e\cos \theta \right)}\left(e\sin \theta *D+\left(1+e\cos \theta \right)*E\right)\\ -\frac{\mathrm{pv}}{h\left(1+e^2+2e\cos \theta \right)}\left(-\left(1+e\cos \theta \right)*D+e\sin \theta *E\right)\end{array}\right];& \left(14\right)\\ \mathrm{where}& \\ \{\begin{array}{c}D=\begin{array}{c}-\frac{3}{2}\frac{G{M}_e}{r^2}{\left(\frac{R_e}{r}\right)}^2\left(1-3{\sin }^{2}i{\sin }^2\left(\omega +\theta \right)\right)-\\ \frac{1}{2}{K}_D\rho v\sqrt{\frac{G{M}_e}{p}}e\sin \theta \end{array}\\ E=\begin{array}{c}-\frac{3}{2}\frac{G{M}_e}{r^2}{\left(\frac{R_e}{r}\right)}^2{\sin }^{2}i\sin \left(2\left(\omega +\theta \right)\right)-\\ \frac{1}{2}{K}_D\rho v\left(\sqrt{\frac{G{M}_e}{p}}\frac{p}{r}-\sigma r\cos i\right)\end{array}\\ F=\begin{array}{c}-\frac{3}{2}\frac{G{M}_e}{r^2}{\left(\frac{R_e}{r}\right)}^2\sin 2i\sin \left(\omega +\theta \right)-\\ \frac{1}{2}{K}_D\rho v\sigma r\sin i\cos \left(\theta +\omega \right)\end{array}\end{array};& \left(13\right)\end{array}$ where η is small thruster efficiency, P is small thruster power, m is spacecraft mass, g is acceleration of gravity, I.sub.sp is small thruster specific impulse, GM.sub.e is gravitational coefficient of the earth, having a value of 3.986×10.sup.14, R.sub.e is radius of the earth, having a size of 6378137 m, $K_D=C_D\frac{A}{m}\sqrt{k_r},C_D,\frac{A}{m}$ is the satellite drag coefficient and area-to-mass ratio, $k_r=1-\frac{2\sigma h\cos i}{v^2},\sigma =7.292\times 10^{-5}$ is rotation velocity of the earth, ρ is atmospheric density at the altitude of the orbit on which the satellite is located, v is tangential velocity, r is radial distance, θ is true anomaly, ω is argument of perigee, i is orbital inclination.

6. The method according to claim 2, wherein in step S14, the objective function J shown in equation (16) is constructed:
J=f(x.sub.i.sup.k)+Σ.sub.j=1.sup.6λ.sub.jh.sub.j(x.sub.i.sup.k)+Σ.sub.j=1.sup.6r.sub.j[h.sub.j(x.sub.i.sup.k)].sup.2  (16); where h.sub.j(x.sub.i.sup.k)=|σ.sub.j−σ.sub.j*|, j=1, 2, . . . , 6 (17); h.sub.j(x.sub.i.sup.k) is an absolute value of the difference between the current orbital elements and the objective orbital elements generated by the control rate of a particle i in the k-th iteration, μ.sub.j,r.sub.j are the Lagrangian multiplier and the penalty factor respectively, f(x.sub.i.sup.k) is deorbit duration generated by the control rate of the particle i in the k-th iteration, which is obtained by orbital integration of the control rate based on the position parameter and the perturbation acceleration.

7. The method according to claim 6, wherein in the process of orbit integration, if the orbit altitude increases and/or the eccentricity deviates from a threshold, the obtained integration result is deleted.

8. The method according to claim 2, wherein in step S2, the position parameter and the velocity of each particle are updated by equation (18): $\begin{array}{cc}\{\begin{array}{c}{x}_i^{k+1}=\omega v_i^k+p_g^k+{\rho }_k\left(1-2m_1\right)\\ v_i^{k+1}=-x_i^k+p_g^k+\omega v_i^k+{\rho }_k\left(1-2m_1\right)\end{array};& \left(18\right)\\ \mathrm{where}& \\ \{\begin{array}{cc}\omega ={\omega }_{\max }-\frac{{\omega }_{\max }-{\omega }_{\min }}{k_{\max }^2}{k}^2& \\ {\rho }_{k+1}=\{\begin{array}{cc}2{\rho }_k& n_s>s\\ 0.5{\rho }_k& n_f>f\\ 2{\rho }_k& \mathrm{other}\end{array}& \end{array};& \left(19\right)\end{array}$ wherein the position and the velocity of a certain particle in the k-th iteration are x.sub.i.sup.k and v.sub.i.sup.k respectively, ω is inertia weight, p.sub.i.sup.k is the best historical position of an individual in the past k iterations, p.sub.g.sup.k is the best population historical position of the particle swarm after k iterations, m.sub.1 is a random number between [0,1], ρ.sub.k is a random direction parameter, n.sub.s,n.sub.f represents the relationship between current best population historical position p.sub.g.sup.k and that of the previous iteration p.sub.g.sup.k-1, and if both are the same, n.sub.f is accumulated by one, otherwise n.sub.s is accumulated by one, s, f represent thresholds of n.sub.s, n.sub.f respectively, ω.sub.max is maximum value of inertia weight, ω.sub.min is minimum value of inertia weight, k.sub.max is maximum number of iterations; wherein in each previous iteration, an individual particle is compared with an individual historical value, and the particle swarm is compared with a population historical value, so that respective particle positions p.sub.i.sup.k and p.sub.g.sup.k corresponding to the minimum fitness are obtained.

9. The method according to claim 8, wherein in step S14, the Lagrangian multiplier and the penalty factor are updated by equation (20): $\begin{array}{cc}\{\begin{array}{c}{\lambda }_j^{n+1}={\lambda }_j^n+2r_j^n{h}_j\left(x_{\mathrm{best}}^n\right)\\ r_j^{n+1}=\{\begin{array}{cc}2r_j^n& {h_j}^2\left(x_{\mathrm{best}}^n\right)>{h_j}^2\left(x_{\mathrm{best}}^{n-1}\right).Math.{h_j}^2\left(x_{\mathrm{best}}^n\right)>{\mathrm{.Math.}}_f\\ 0.5r_j^n& {h_j}^2\left(x_{\mathrm{best}}^n\right)\le {\mathrm{.Math.}}_f\\ r_j^n& \mathrm{other}\end{array}\end{array}& \left(20\right)\end{array}$ wherein value of x.sub.best is equal to the population optimal fitness p.sub.g.sup.k obtained in step S2, ε.sub.f is a constraint threshold, λ.sub.j.sup.n,r.sub.j.sup.n are the Lagrangian multiplier and the penalty factor respectively, where n is sequence number, h.sub.j(x.sub.i.sup.k) is an absolute value of the difference between the current orbital elements and the objective orbital elements generated by the control rate of a particle i in the k-th iteration.

10. The method according to claim 2, wherein the orbital elements comprise semi-major axis a, eccentricity e, orbital inclination i, right ascension of the ascending node Ω, argument of perigee ω and true anomaly θ.

11. The method according to claim 10, wherein the objective fitness and the population fitness are fitnesses corresponding to a minimum value of the objective function.

12. A readable storage medium storing one or more readable programs adapted to cause a processor is executable by one or more processors to implement the steps of a low-orbit satellite deorbit control method based on a particle swarm algorithm according to claim 1.

13. A computer device comprising a memory configured to store computer programs thereon, and a processor configured to execute the computer programs stored on the memory to implement the steps of a low-orbit satellite deorbit control method based on a particle swarm algorithm according to claim 1.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0045] The accompanying drawings, which are incorporated in and constitute a part of this specification, illustrate exemplary embodiments of the present invention, and together with the description serve to explain the principles of the invention. Obviously, the following descriptions with reference to the drawings are merely descriptions of exemplary embodiments of the invention. For those skilled in the art, all other drawings may be obtained from those drawings without creative efforts.

[0046] FIG. 1 is a flowchart illustrating a low-orbit satellite deorbit control method based on a particle swarm algorithm according to an exemplary embodiment of the disclosure;

[0047] FIG. 2 is an outer block diagram illustrating an augmented Lagrangian particle swarm algorithm according to an exemplary embodiment of the disclosure;

[0048] FIG. 3 is a block diagram illustrating a sub-program of a particle swarm algorithm according to an exemplary embodiment of the disclosure;

[0049] FIG. 4 is a schematic diagram illustrating the thrust direction of a small thruster according to an exemplary embodiment of the disclosure;

[0050] FIG. 5 shows the result of variation of semi-major axis over time under natural conditions according to an exemplary embodiment of the disclosure;

[0051] FIG. 6 shows the iterative process of a first disposal orbit routine according to an exemplary embodiment of the disclosure;

[0052] FIG. 7 shows the iterative process of a second disposal orbit routine according to an exemplary embodiment of the disclosure;

[0053] FIG. 8 shows the change of the semi-major axis, eccentricity, and orbit inclination of the satellite under the optimal control rate according to an exemplary embodiment of the disclosure (the first disposal orbit); and

[0054] FIG. 9 shows the change of the semi-major axis, eccentricity, and orbital inclination of the satellite under the optimal control rate according to an exemplary embodiment of the disclosure (the second disposal orbit).

DETAILED DESCRIPTION OF THE EMBODIMENTS

[0055] The technical problems to be solved and technical solutions and beneficial effects of the present disclosure will become apparent from the following detailed description with reference to the accompanying drawings and embodiments. It should be noted that the exemplary embodiments described herein are merely intended to explain the present invention, rather than limit the scope of the present invention.

[0056] The exemplary embodiment of the disclosure provides a low-orbit satellite deorbit control method based on a particle swarm algorithm. Firstly, three dimensions of particles are used to represent three costate variables, and all particles in the population are randomly initialized within a feasible region and a Lagrangian multiplier and a penalty factor are initialized respectively. Subsequently, an objective function is used as a fitness function of the particle swarm algorithm, a fitness function value is calculated for each particle, an optimal fitness of each particle and a population optimal fitness are updated, and the position and velocity of each of the particles are updated until a maximum number of iterations is reached; Finally, augmented Lagrangian function parameters are update, a new objective function is generate, the previous step is repeat until the number of iterations and constraint threshold requirements are met, and an optimization result is updated. This process combines the Hamiltonian function and the augmented Lagrangian particle swarm algorithm. Under the augmented Lagrangian function, the particle swarm algorithm can obtain costate variables corresponding to the shortest deorbit duration after multiple iterations, thereby achieving the optimal control rate.

[0057] Illustratively, FIG. 1 shows a flowchart of a low-orbit satellite deorbit control method based on a particle swarm algorithm according to an exemplary embodiment of the disclosure; and FIG. 2 shows an outer block diagram of an augmented Lagrangian particle swarm algorithm according to an exemplary embodiment of the disclosure. The method of embodiments of the disclosure will be described below with reference to FIGS. 1 and 2.

[0058] At S1, an objective function is constructed according to a position parameter and a perturbation acceleration of each of particles in a particle swarm as well as a preset Lagrangian multiplier and a penalty factor.

[0059] The method further includes the following steps.

[0060] At step S11, a rate of change of each orbital element over time under the thrust applied by a thruster is obtained, and a control rate containing current costate variables is determined according to the rate of change.

[0061] In this embodiment, the orbital elements include semi-major axis a, eccentricity e, orbital inclination i, right ascension of the ascending node Ω, argument of perigee ω and true anomaly θ, and the control rate containing current costate variables is obtained by the following steps S111, S112 and S113.

[0062] At step S111, the rate of change of each orbital element over time under the thrust applied by the thruster is obtained. In this case, the influence of the thrust by the thruster on the orbital element is calculated by:

[00007]$\begin{array}{cc}\frac{da}{dt}=\frac{2a^{2}v}{\mu }{a}_t;& \left(1\right)\\ \frac{de}{dt}=\frac{1}{v}\left(2\left(e+\cos \theta \right)a_t+\frac{r}{a}{a}_n\sin \theta \right);& \left(2\right)\\ \frac{di}{dt}=r{h}^{-1}{a}_h\cos \left(\omega +\theta \right);& \left(3\right)\\ \frac{d\Omega }{dt}=\frac{r\sin \left(\omega +\theta \right)}{h\sin i}{a}_h;& \left(4\right)\\ \frac{d\omega }{dt}=-\frac{d\Omega }{dt}\cos i+\frac{1}{ev}\left(2a_t\sin \theta -\left(2e+\frac{r}{a}\cos \theta \right)a_n\right);& \left(5\right)\\ \frac{d\theta }{dt}=\frac{h}{r^2}-\frac{1}{ev}\left(2a_t\sin \theta -\left(2e+\frac{r}{a}\cos \theta \right)a_n\right);\mathrm{where}p=a\left(1-e^2\right),h=\sqrt{\mu p},v=\sqrt{\frac{2\mu }{r}-\frac{\mu }{a}},r=\frac{p}{\left(1+e\cos \theta \right)},\mu =3.986\times 10^{14},& \left(6\right)\end{array}$

where p, h, r, v, μ are indicated in these equations, v is satellite velocity, μ is gravitational coefficient of the earth, a.sub.n,t,h is perturbation acceleration of occlusive orbit in the radial direction (radially outward), in the velocity direction (along the direction of velocity) and in the angular momentum direction (the direction of angular momentum of the occlusive orbit may be positive).

[0063] At S112, based on the rate of change and the current costate variables, a Hamiltonian function is obtained by equation (7).

[0064] In this embodiment, the Hamiltonian function is calculated by:

[00008]$\begin{array}{cc}H={\lambda }_a\frac{da}{dt}+{\lambda }_e\frac{de}{dt}+{\lambda }_i\frac{di}{dt};& \left(7\right)\end{array}$

[0065] where H is a Hamiltonian function value, λ.sub.a, λ.sub.e and λ.sub.i are the current costate variables, a, e and i are semi-major axis, eccentricity and orbital inclination of the operational orbit respectively, and t is time.

[0066] At step S113, a preset angle corresponding to the thrust applied by the thruster is derived according to the Hamilton function (7) to obtain the control rate containing the current costate variables.

[0067] Herein, FIG. 4 shows a schematic diagram of the thrust direction of a small thruster according to an exemplary embodiment of the disclosure. The process of determining the control rate in the embodiment will be described below with reference to FIG. 4. During the process, a control rate containing current costate variables at the angles α and β as shown in FIG. 4 is calculated, where the partial derivatives of the Hamilton function may be taken with respect to α and β respectively, and a plurality of simultaneous equations may be set as follows:

[00009]$\begin{array}{cc}A={\lambda }_e\frac{r}{a}\sin \theta B={\lambda }_a\frac{a^2{v}^2}{\mu }+{\lambda }_e\left(e+\cos \theta \right)\sin {\alpha }^{*}=\frac{-A}{\sqrt{4B^2+A^2}}\cos {\alpha }^{*}=\frac{-2B}{\sqrt{4B^2+A^2}};& \left(8\right)\\ C={\lambda }_i\frac{rv}{h}\cos \left(\omega +\theta \right)D={\lambda }_a\frac{2a^2{v}^2}{\mu }\cos {\alpha }^{*}E={\lambda }_e\left(2\left(e+\cos \theta \right)\cos {\alpha }^{*}+\frac{r}{a}\sin \theta \sin {\alpha }^{*}\right)\sin {\beta }^{*}=\frac{-C}{\sqrt{C^2+D^2+E^2}}\cos {\beta }^{*}=\frac{-\left(D+E\right)}{\sqrt{C^2+D^2+E^2}}.& \left(9\right)\end{array}$

[0068] At step S12, a thrust vector, a fuel consumption rate, and a perturbation acceleration of the thrust applied by the thruster are calculated.

[0069] In this embodiment, the vector of the small thrust and the fuel consumption rate shown in FIG. 4 are calculated, and the non-spherical perturbation and the perturbation acceleration of the atmospheric drag perturbation are calculated. J.sub.2

[0070] At step S121, the thrust vector of the thruster is calculated by:

[00010]$\begin{array}{cc}{\begin{array}{ccc}[a_n& a_t& {a_h]}^T=.Math.\stackrel{.fwdarw.}{F}\end{array}}_{\mathrm{thrust}}.Math.[\begin{array}{ccc}\sin \alpha \cos \beta & \cos \alpha \cos \beta & {\sin \beta ]}^T;\end{array}& \left(10\right)\\ .Math.{\stackrel{.fwdarw.}{F}}_{\mathrm{thrust}}.Math.=\frac{2\eta P}{{\mathrm{mgI}}_{\mathrm{sp}}};& \left(11\right)\end{array}$

[0071] The fuel consumption rate is calculated by:

[00011]$\begin{array}{cc}\frac{dm}{dt}=\frac{-2\eta P}{{\left({\mathrm{gI}}_{sp}\right)}^2};& \left(12\right)\end{array}$

[0072] where η is small thruster efficiency, P is small thruster power, m is spacecraft mass, g is acceleration of gravity, and I.sub.sp is small thruster specific impulse.

[0073] At step S122, the nonspherical J.sub.2 perturbation and the perturbation acceleration of the atmospheric drag perturbation are calculated.

[0074] The perturbation acceleration is calculated by:

[00012]$\begin{array}{cc}\{\begin{array}{c}D=\begin{array}{c}-\frac{3}{2}\frac{G{M}_e}{r^2}{\left(\frac{R_e}{r}\right)}^2\left(1-3{\sin }^{2}i{\sin }^2\left(\omega +\theta \right)\right)-\\ \frac{1}{2}{K}_D\rho v\sqrt{\frac{G{M}_e}{p}}e\sin \theta \end{array}\\ E=\begin{array}{c}-\frac{3}{2}\frac{G{M}_e}{r^2}{\left(\frac{R_e}{r}\right)}^2{\sin }^{2}i\sin \left(2\left(\omega +\theta \right)\right)-\\ \frac{1}{2}{K}_D\rho v\left(\sqrt{\frac{G{M}_e}{p}}\frac{p}{r}-\sigma r\cos i\right)\end{array}\\ F=\begin{array}{c}-\frac{3}{2}\frac{G{M}_e}{r^2}{\left(\frac{R_e}{r}\right)}^2\sin 2i\sin \left(\omega +\theta \right)-\\ \frac{1}{2}{K}_D\mathrm{\rho v\sigma }r\sin i\cos \left(\theta +\omega \right)\end{array}\end{array};& \left(13\right)\\ \left[\begin{array}{c}{a}_t\\ a_n\\ a_h\end{array}\right]=\left[\begin{array}{c}\begin{array}{c}\frac{\mathrm{pv}}{h\left(1+e^2+2e\cos \theta \right)}(e\sin \theta *D+\\ \left(1+e\cos \theta \right)*E)\end{array}\\ \begin{array}{c}-\frac{\mathrm{pv}}{h\left(1+e^2+2e\cos \theta \right)}(-\\ \left(1+e\cos \theta \right)*D+e\sin \theta *E)\end{array}\\ F\end{array}\right];& \left(14\right)\end{array}$

[0075] where, GM.sub.e is gravitational coefficient of the earth, having a value of 3.986×10.sup.14, R.sub.e is radius of the earth, having a size of 6378137 m,

[00013]$K_D=C_D\frac{A}{m}\sqrt{k_r},C_D,\frac{A}{m}$

is the satellite drag coefficient and area-to-mass ratio,

[00014]$k_r=1-\frac{2\sigma h\cos i}{v^2},$

σ=7.292×10.sup.−5 is rotation velocity of the earth, ρ is atmospheric density at the altitude of the orbit on which the satellite is located, v is tangential velocity, r is radial distance, θ is true anomaly, ω is argument of perigee, and i is orbital inclination.

[0076] At step S13, a particle position vector is determined according to the current costate variables, and the particle position vector is used as the position parameter of each of the particles in the particle swarm.

[0077] In this embodiment, the costate variable is used as the particle position vector, and the particle position vector is calculated by:

λ.sub.i=[λ.sub.aλ.sub.eλ.sub.i]  (15).

[0078] At step S14, orbit integration of the control rate is performed based on the position parameter and the perturbation acceleration, and the objective function is constructed according to a integration result, the preset Lagrangian multiplier and the penalty factor;

[0079] In this embodiment, the objective function J shown in equation (16) is constructed:

J=f(x.sub.i.sup.k)+Σ.sub.j=1.sup.6λ.sub.jh.sub.j(x.sub.i.sup.k)+Σ.sub.j=1.sup.6r.sub.j[h.sub.j(x.sub.i.sup.k)].sup.2  (16);

[0080] where h.sub.j(x.sub.i.sup.k)=|σ.sub.j−σ.sub.j*|, j=1, 2, . . . , 6 (17); h.sub.j(x.sub.i.sup.k) is an absolute value of the difference between the current orbital elements and the objective orbital elements generated by the control rate of a particle i in the k-th iteration, λ.sub.j,r.sub.j are the Lagrangian multiplier and the penalty factor respectively, f(x.sub.i.sup.k) is deorbit duration generated by the control rate of the particle i in the k-th iteration, which is obtained by orbital integration of the control rate based on the position parameter and the perturbation acceleration (implemented by integration of equations (1) to (6)). λ.sub.j,r.sub.j remains unchanged during the sub-problem solving operation.

[0081] It should be noted that some constraints may be implemented in the orbital integration process through corresponding processing. Since the initial position value of each of the particles is random, some particles may cause abnormal conditions such as an increase in orbital altitude and abnormal eccentricity during orbital integration. To solve such a problem, an interruption condition can be set during the orbital integration process and the deorbit duration is returned as an infinite value. The fitness value of each of these particles also becomes infinite, which is eliminated in fitness screening.

[0082] At step S2, an objective function value of each particle in the particle swarm is calculated based on the objective function, an objective fitness of each particle and a population fitness of the particle swarm are determined according to a calculation result, the position parameter and the velocity of each particle are updated, and a population optimal fitness at a maximum number of iterations and the position parameters of particles corresponding to the population optimal fitness is obtained, wherein the population optimal fitness is a population fitness that is output when the preset maximum number of iterations is met.

[0083] In this embodiment, the objective fitness and the population fitness are fitnesses corresponding to the minimum value of the objective function.

[0084] The optimal fitness of each particle p.sub.i.sup.k and the population optimal fitness p.sub.g.sup.k are obtained by comparing respective objective function values of respective particles in the population (the smallest is the optimum), the position and the velocity of each particle are updated, and an optimization result of the subroutine is output when the preset number of iterations is met (see FIG. 3).

[0085] In this case, the position parameter and the velocity of each particle are updated by equation (18):

[00015]$\begin{array}{cc}\{\begin{array}{c}{x}_i^{k+1}=\omega v_i^k+p_g^k+{\rho }_k\left(1-2m_1\right)\\ v_i^{k+1}=-x_i^k+p_g^k+\omega v_i^k+{\rho }_k\left(1-2m_1\right)\end{array};& \left(18\right)\\ \mathrm{where}& \\ \{\begin{array}{c}\omega ={\omega }_{\max }-\frac{{\omega }_{\max }-{\omega }_{\min }}{k_{\max }^2}{k}^2\\ {\rho }_{k+1}=\{\begin{array}{cc}2{\rho }_k& n_s>s\\ 0.5{\rho }_k& n_f>f\\ 2{\rho }_k& \mathrm{other}\end{array}\end{array};& \left(19\right)\end{array}$

[0086] wherein the position and velocity of a certain particle in the k-th iteration are x.sub.i.sup.k and v.sub.i.sup.k respectively, ω is inertia weight, p.sub.i.sup.k is the best historical position of the individual in the past k iterations, p.sub.g.sup.k is the best population historical position of the particle swarm after k iterations, m.sub.1 is a random number between [0,1], ρ.sub.k is a random direction parameter, n.sub.s,n.sub.f represents the relationship between current best population historical position p.sub.g.sup.k and that of the previous iteration p.sub.g.sup.k-1, and if both are the same, n.sub.f is accumulated by one, otherwise n.sub.s is accumulated by one, s, f represent thresholds of n.sub.s,n.sub.f respectively, ω.sub.max is maximum value of inertia weight, ω.sub.min is minimum value of inertia weight, k.sub.max is maximum number of iterations;

[0087] wherein in each previous iteration, an individual particle is compared with an individual historical value, and the particle swarm is compared with a population historical value, so that respective particle positions p.sub.i.sup.k and p.sub.g.sup.k corresponding to the minimum fitness are obtained.

[0088] At step S3, the Lagrangian multiplier and the penalty factor are updated according to the population optimal fitness, a current number of iterations and an objective deviation value are obtained based on the population optimal fitness, the current number of iterations and the objective deviation value are compared with a preset iteration condition and a preset convergence condition, and an objective population fitness and an objective position parameter are determined according to a comparison result to obtain an objective deorbit mode. In this case, the objective deviation value based on the population optimal fitness is the deviation value between the semi-major axis or the perigee generated by the position parameter integration corresponding to the population optimal fitness obtained in step S2 and the target.

[0089] In this embodiment, if the current number of iterations meets the preset iteration condition and the objective deviation value meets the convergence condition, the population optimal fitness is determined as the objective population fitness, and the position parameter corresponding to the objective population fitness is determined as the objective position parameter.

[0090] If the current number of iterations does not meet the preset iteration condition or the objective deviation value does not meet the convergence condition, the objective function is updated based on the updated Lagrangian multiplier and the penalty factor, and step S14 and subsequent steps are continued to be performed, until the current number of iterations meets the preset iteration condition and the objective deviation value meets the convergence condition.

[0091] In this case, the Lagrange multiplier λ.sub.j.sup.n and the penalty factor r.sub.j.sup.n is updated by calculation of:

[00016]$\begin{array}{cc}\{\begin{array}{c}{\lambda }_j^{n+1}={\lambda }_j^n+2r_j^n{h}_j\left(x_{\mathrm{best}}^n\right)\\ r_j^{n+1}=\{\begin{array}{cc}2r_j^n& {h_j}^2\left(x_{\mathrm{best}}^n\right)>{h_j}^2\left(x_{\mathrm{best}}^{n-1}\right).Math.{h_j}^2\left(x_{\mathrm{best}}^n\right)>{\mathrm{.Math.}}_f\\ 0.5r_j^n& {h_j}^2\left(x_{\mathrm{best}}^n\right)\le {\mathrm{.Math.}}_f\\ r_j^n& \mathrm{other}\end{array}\end{array}& \left(1\right)\end{array}$

where the value of x.sub.best.sup.n is equal to the population optimal fitness p.sub.g.sup.k obtained in step S2 (an optimal solution to the sub-problem shown in FIG. 3 is obtained), ε.sub.f is a constraint threshold, λ.sub.j.sup.n,r.sub.j.sup.n are the Lagrangian multiplier and the penalty factor respectively, where n is sequence number, h.sub.j(x.sub.i.sup.k) is an absolute value of the difference between the current orbital elements and the objective orbital elements generated by the control rate of a particle i in the k-th iteration.

[0092] In this case, if the optimal solution x.sub.best.sup.n does not satisfy the constraint condition, the procedure returns to step S14, a new objective function is calculated by the updated Lagrangian multiplier and penalty factor, and then the procedure proceeds to subsequent steps until the constraint condition is satisfied ε.sub.f.

[0093] In the low-orbit satellite deorbit control method based on the particle swarm algorithm of the disclosure, firstly a control rate with costate variables is derived from a Hamiltonian function; subsequently, costate variable vectors are used as position parameters of particles in a particle swarm to operate a particle swarm algorithm subroutine with a current objective function (FIG. 3), an optimal fitness p.sub.i.sup.k of each particle, a population optimal fitness p.sub.g.sup.k and the position and velocity of each of the particles are updated in each iteration, and an optimization result is output when a maximum number of iterations is met; finally, according to the optimization result obtained by the subroutine, parameters in the augmented Lagrangian function are updated, and then a new objective function is obtained for the next subroutine iteration, and an optimization result is output when the number of iterations and constraint threshold requirements are met. The method can be utilized to obtain an optimal control rate that minimizes the low-thrust deorbit duration of low-orbit satellites, which introduces the advantages of insensitivity to initial values and fast convergence speed, and provides an effective solution for the engineering application of low-orbit satellite deorbit with small thrust.

Embodiment 1

[0094] The low-orbit satellite deorbit control method based on the particle swarm algorithm of the disclosure will be described below in conjunction with this embodiment. The embodiment may be achieved by the following steps: at the first step, parameters of a spacecraft itself and parameters of a small thruster are given; at the second step, semi-major axis, eccentricity and orbital inclination of an operational orbit and two disposal orbits of the spacecraft are given; at the third step, parameters of an augmented Lagrangian particle swarm algorithm are given; at the fourth step, deorbit duration without thrust is calculate; at the fifth step, an optimal control rate, deorbit duration and fuel consumption of the two disposal orbits are calculated.

[0095] Herein the parameters of the spacecraft itself and the parameters of the small thruster are shown in Table 1, the semi-major axis, eccentricity and orbit inclination on an operational orbit and two disposal orbits of the spacecraft are shown in Table 2, the parameters of the augmented Lagrangian particle swarm algorithm are shown in Table 3, and the optimal control rate, deorbit duration and fuel consumption of the two disposal orbits are shown in Table 4.

TABLE-US-00001 TABLE 1 Spacecraft Parameters Numerical Numerical parameter value parameter value Mass of Satellite 1000.0 kg Specific Impulse 4660 s Atmospheric Drag Coefficient 2 Thruster Power 115 W Area-to-Mass Ratio 0.02 m.sup.2/kg Engine Thrust 3 mN Thruster efficiency 24.5% Atmospheric SA76 model

TABLE-US-00002 TABLE 2 Operational Orbit and Disposal Orbit Parameters Second Operational First Disposal Disposal Orbit Elements Orbit Orbit Orbit Semi-Major Axis (km) 7200 6528 7200 Eccentricity 0.001 0.001 0.093 Orbital Inclination (degrees) 98.5 98.5 98.5

TABLE-US-00003 TABLE 3 Augmented Lagrangian Particle Swarm Algorithm Parameters Algorithm Parameter Numerical Value Number of Populations N 500 Maximum Number of Iterations k.sub.max 50 Maximum Inertia Weight ω.sub.max 0.9 Minimum Inertia Weight ω.sub.min 0.4 Threshold of n.sub.s s 15 Threshold of n.sub.f f 5 Initial Value of Random Direction Parameter ρ.sub.0 1 Initial Penalty Factor r.sub.0 1000 Constraint Threshold ε.sub.f 0.0001

TABLE-US-00004 TABLE 4 Optimal Control Rate, Deorbit Duration and Fuel Consumption on Two Disposal Orbits Deorbit Fuel Parameter λ.sub.a λ.sub.e λ.sub.i Duration Consumption Particle Swarm [−10/a, 10/a] [−10, 10] [−10, 10] Search Range (Type 1) Particle Swarm [−10/a, 10/a] [−1 × 10.sup.6, 1 × 10.sup.6] [−10, 10] Search Range (Type 2) Optimal −1.389 × 10.sup.−6 −4.074 6.563  857 days 2.00 kg Solution (First Disposal Orbit) Optimal −1.389 × 10.sup.−7 −2223.345 10 1607 days 3.75 kg Solution (Second Disposal Orbit)

[0096] FIG. 5 is a graph showing variation of semi-major axis over time without thrust according to Embodiment 1 of the disclosure. As seen from FIG. 5, the deorbit duration without thrust under natural conditions will be 140 years, far beyond the 25 years deorbit requirement. FIGS. 6 and 7 show the results of iterations in the two disposal orbits respectively. FIG. 6 illustrates a total of 60 iterations in the process, and the semi-major axis deviation convergence condition will be met at the 57th iteration. FIG. 7 illustrates a total of 25 iterations in the process, and the perigee deviation convergence condition will be met at the 12th iteration. FIGS. 8 and 9 show graphs of variation of semi-major axis, eccentricity and orbital inclination over time on the two disposal orbit respectively. As seen from FIGS. 8 and 9, the spacecraft can reach the corresponding objective orbit (disposal orbit) under the optimal control rate obtained by the disclosure. Based on this effect, the effectiveness of the optimal control rate obtained by the disclosed algorithm is proved.

[0097] Another exemplary embodiment of the disclosure provides a low-orbit satellite deorbit control system based on a particle swarm algorithm, including: a construction unit configured to construct an objective function according to a position parameter and a perturbation acceleration of each of particles in a particle swarm as well as a preset Lagrangian multiplier and a penalty factor; an updating unit configured to calculate an objective function value of each particle in the particle swarm based on the objective function, determine an objective fitness of each particle and a population fitness of the particle swarm according to a calculation result, update the position parameter and the velocity of each particle, and obtain a population optimal fitness at a maximum number of iterations and the position parameters of particles corresponding to the population optimal fitness; and a processing unit configured to update the Lagrangian multiplier and the penalty factor according to the population optimal fitness, obtain a current number of iterations and an objective deviation value based on the population optimal fitness, compare the current number of iterations and the objective deviation value with a preset iteration condition and a preset convergence condition, and determine an objective population fitness and an objective position parameter according to a comparison result to obtain an objective deorbit mode.

[0098] The detailed description of each unit in the aforementioned control system can refer to the description of each step in the control method, which will not be repeatedly described. The aforementioned control system can achieve the same function as the control method.

[0099] On the basis of the above-mentioned embodiments, an embodiment of the disclosure provides a computer readable storage medium storing one or more readable programs adapted to cause a processor is executable by one or more processors to implement the steps of a low-orbit satellite deorbit control method based on a particle swarm algorithm as described in the first embodiment, for example.

[0100] Optionally, the aforementioned storage medium may be a non-transitory computer-readable storage medium. For example, the non-transitory computer-readable storage medium may be ROM, random access memory (RAM), CD-ROM, magnetic tape, floppy disk, and optical data storage device, etc.

[0101] Additionally, another embodiment of the disclosure further provides a computer device, including a processor and a memory.

[0102] The memory is configured to store computer programs therein.

[0103] The processor is configured to execute the computer programs stored in the memory to implement a low-orbit satellite deorbit control method based on a particle swarm algorithm as described in the first embodiment, for example.

[0104] The aforementioned processor may be a general-purpose processor, including central processing unit (CPU), network processor (NP), etc.; it may also be digital signal processor (DSP), application specific integrated circuit (ASIC), field-programmable gate array (FPGA) or other programmable logic devices, discrete gates or transistor logic devices, discrete hardware components.

[0105] Accordingly, the above description is only considered as one embodiment of the disclosure, rather than limitation of the invention in any way. Although the invention has been described in the above with the preferred embodiments, the invention is not limited to those. Those skilled in the art, without departing from the scope of the invention, can make some equivalent changes or modifications using the content disclosed above, which will fall within the scope of the invention.