Atmospheric differential drag phase adjustment control method for low-orbit micro-nano satellite constellation

Abstract

The application discloses a method for adjusting and controlling the phase of a low-orbit micro-nano satellite constellation by means of atmospheric differential drag, and belongs to the field of on-orbit phase control of a constellation. The method aims at solving the problem of low phase control precision of the existing atmospheric differential drag control satellite constellation method. The method comprises the following steps: establishing a single-satellite atmospheric drag dynamics model to obtain the corresponding relationship between the single-satellite control input and the control command; determining the state constraint condition of other single satellites based on a reference satellite, combining the energy constraint and the task constraint, and adopting a simulated annealing algorithm to perform constellation slot allocation, and taking the shortest control time as the target to obtain the optimal result of the constellation slot allocation; calculating an initial control command matrix; taking the lowest adaptive function value and meeting the target control precision as the control target, optimizing the initial control command matrix to obtain an optimal control command matrix, so as to obtain the optimal satellite control input, and performing the angular acceleration control of the single satellite, thereby realizing the phase adjustment of the same-orbit constellation. The application is used for the phase adjustment control of the constellation.

Description

TECHNICAL FIELD

[0001] The application relates to an atmospheric differential drag phase adjustment control method for a low-orbit micro-nano satellite constellation, and belongs to the field of on-orbit phase control of the constellation. BACKGROUND

[0002] With the development of micro-satellites, low-orbit micro-satellite constellations have become a development trend. The development of satellite Internet is promoted, and the needs of a space infrastructure system can be met. In the phase control of a low-orbit micro-satellite constellation, the mass, cost, propellant consumption and atmospheric drag received by a low-orbit spacecraft become main problems.

[0003] Under the current technical level, considering the adaptability of a carrier rocket and the implementability of a batch production project, the deployment of the constellation is generally implemented by using the same type of carrier rocket and the same launch site, and the launch site launch cycle constraints affect the pace of constellation deployment. Meanwhile, considering the on-orbit operation reliability factors of the satellite and the need to enter a working orbit as soon as possible after entering the orbit to save the on-orbit life, a deployment mode of directly entering a predetermined orbital plane is generally adopted, and long-term large-scale orbital plane drift is not performed. In addition, for a near-earth low-orbit satellite system, considering that the adjustment direction of the phase before and after entering the orbit may change with the launch strategy, in order to avoid frequent changes of launch elements, the carrier rocket usually directly aims at the deployment of the satellite predetermined orbital height, and the satellite enters the orbit and relies on its own orbit control to perform the orbit transfer process, which will consume a large amount of remote control operation, human and time costs. In addition, due to the orbital characteristics of the low-orbit satellite, the spacecraft is greatly affected by the atmospheric perturbation during the deployment of the constellation, and additional propellant consumption is needed to eliminate this influence.

[0004] After Leonard et al. proposed the method of using atmospheric differential drag to control the motion of a satellite, researchers have improved the control model, optimized the control method and researched in different fields. Research shows that atmospheric differential drag control is an effective method for formation control, and the on-orbit test of the American Dove satellite constellation proves its engineering feasibility. However, due to imperfect on-orbit situation analysis, there are differences between real data and simulation, and the control precision cannot meet the requirements in engineering practice when used in engineering practice. SUMMARY

[0005] In view of the low phase control precision of the existing atmospheric differential drag control satellite constellation method, the application provides an atmospheric differential drag phase adjustment control method for a low-orbit micro-nano satellite constellation.

[0006] The atmospheric differential drag phase adjustment control method for a low-orbit micro-nano satellite constellation provided by the application comprises,

[0007] For a same-orbit micro-nano satellite constellation operating in a low-earth orbit, a single-satellite atmospheric drag dynamics model is established, and the corresponding relationship between the single-satellite control input and the control command is obtained.

[0008] The reference satellite is determined in the satellite constellation, the state constraint condition of other single satellites is determined based on the reference satellite, the constellation slot allocation is carried out by using the simulated annealing algorithm in combination with the energy constraint and the task constraint of the single satellite, the constellation slot allocation optimization is carried out with the shortest control time as the target, and the optimal result of the constellation slot allocation is obtained.

[0009] The initial control command matrix is calculated based on the optimal result of the constellation slot allocation, the adaptive function is set, the initial control command matrix is optimized with the lowest adaptive function value and the target control precision as the control target, and the optimal control command matrix is obtained.

[0010] According to the corresponding relationship between the single satellite control input and the control command, the optimal satellite control input is obtained from the optimal control command matrix, the angular acceleration control of the single satellite in the constellation is carried out, the angular acceleration difference of different satellites in the constellation is generated, and the in-orbit phase adjustment of the same-orbit constellation is realized.

[0011] According to the atmospheric differential drag phase adjustment control method of the low-orbit micro-nano satellite constellation, the single satellite atmospheric drag dynamics model is:

[0012]

[0013] In the formula, θ k is the absolute in-orbit phase of the satellite in the discrete time period k, is the absolute in-orbit angular velocity of the satellite in the discrete time period k, Δt is the control time step, k represents the kth discrete time, u k is the control command in the discrete time period k, B k is the control input in the discrete time period k.

[0014] According to the formula (1), the corresponding relationship between the single satellite control input and the control command is:

[0015]

[0016] In the formula, θ is the absolute in-orbit angular acceleration of the satellite in the high-resistance state in the discrete time period k, is the absolute in-orbit angular acceleration of the satellite in the low-resistance state in the discrete time period k, the high-resistance state is the satellite attitude when the satellite windward area is maximum, the low-resistance state is the satellite attitude when the satellite windward area is minimum, k=0, 1, 2, 3, …, end; in the formula, end represents the number of discrete time periods.

[0017] According to the atmospheric differential drag phase adjustment control method of the low-orbit micro-nano satellite constellation, the reference satellite is taken as the No. 1 satellite in the satellite constellation, and the state constraint condition of other single satellites is represented as:

[0018]

[0019] where θ i,0 is the initial absolute phase of the i-th satellite, i = 2, 3, 4, …, N, where N is the number of satellites in the constellation; θ 1,0 is the initial absolute phase of the 1st satellite, θ i ′ ,0 is the initial relative phase of the i-th satellite with respect to the 1st satellite, θ i,end is the final absolute phase of the i-th satellite after the end segment of the discrete-time control, θ 1,end is the final absolute phase of the 1st satellite after the end segment of the discrete-time control, θ i ′ ,f is the target relative phase of the i-th satellite with respect to the 1st satellite.

[0020] The atmospheric differential drag phase adjustment control method of the low-orbit micro-nano satellite constellation according to the application, the energy constraint of a single satellite is:

[0021] u k = α, 0 < α < 1,

[0022] wherein α represents the attitude command of the satellite in the charging window;

[0023] The task constraint of a single satellite is:

[0024] u k = β, 0 < β < 1

[0025] wherein β represents the attitude command of the satellite in the task window;

[0026] The satellite charging window phase interval Φ E and the satellite task window phase interval Φ M are set, so that the on-orbit absolute angular acceleration of the satellite is between the on-orbit absolute angular acceleration in the high resistance state and the on-orbit absolute angular acceleration in the low resistance state:

[0027] The satellite charging window phase interval The satellite task window phase interval wherein is the minimum charging phase, is the maximum charging phase, is the minimum task phase, is the maximum task phase;

[0028] The control command u ik of the i-th satellite in the discrete-time segment k is obtained as:

[0029]

[0030] The low-orbit micro-nano satellite constellation atmospheric differential drag phase adjustment control method according to the application adopts the following method for constellation slot allocation by using the simulated annealing algorithm:

[0031] An annealing initial temperature T0, an annealing temperature lower limit T min , an annealing coefficient ξ are set, and an initial slot allocation vector S0 is generated:

[0032] S0 = [1 2... N],

[0033] The constellation numbers in S0 are arranged in correspondence with the target relative phase numbers;

[0034] The initial slot allocation vector S0 is optimized to obtain the optimal result of constellation slot allocation.

[0035] The application has the following advantages: the method considers satellite energy and task constraints, optimizes atmospheric differential drag control, and realizes phase control closer to engineering practice. The method can avoid consuming propulsion costs. It first establishes atmospheric drag control dynamics and kinematics models, considers satellite energy, on-orbit tasks and other practical engineering constraints, realizes time-optimal control of same-orbit constellation phase distribution based on slot allocation and atmospheric differential drag control. On this basis, the simulated annealing algorithm is used to optimize slot allocation, and a control command matrix is generated based on the allocation result, realizing time-optimal control of same-orbit constellation phase adjustment.

[0036] The method uses atmospheric differential drag control principle to control constellation phase, realizes considering atmospheric perturbation as dynamic application, converts energy required for orbit control into small torque required for attitude control, and provides a constellation phase adjustment method saving energy;

[0037] The method is based on the simulated annealing algorithm, constructs a satellite constellation slot allocator, and optimizes satellite constellation phase adjustment time.

[0038] The method considers the energy constraint conditions and task constraint conditions of the satellite constellation in space, and generates a control instruction matrix for the satellite constellation control under the constraint conditions, thereby being more practical in engineering. BRIEF DESCRIPTION OF DRAWINGS

[0039] Figure 1 is a schematic diagram of generating a control command matrix in the low-orbit micro-nano satellite constellation atmospheric differential drag phase adjustment control method described in the application;

[0040] Figure 2 is a schematic diagram of an idealized constraint window on a satellite orbit;

[0041] Figure 3 is a schematic diagram of the slot allocation optimization process;

[0042] Figure 4 is a constellation control time estimation comparison chart before and after slot assignment optimization;

[0043] Figure 5 is a schematic diagram of the command matrix optimization process under the constraint condition;

[0044] Figure 6 is a simulation schematic diagram of the relative phase of the same-orbit constellation changing over time at the end of the control time. DETAILED DESCRIPTION

[0045] The technical solutions in the embodiments of the present application will be clearly and completely described below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, rather than all the embodiments of the present application. Based on the embodiments in the present application, all other embodiments obtained by a person of ordinary skill in the art without creative work fall within the scope of protection of the present application.

[0046] It should be noted that the embodiments in the present application and the features in the embodiments can be combined with each other without conflict.

[0047] The present application will be further described below in combination with the drawings and specific embodiments, but is not limited to the present application.

[0048] Specific embodiment one, in combination Figure 1 As shown in the figure, the present application provides an atmospheric differential drag phase adjustment control method for a low-orbit micro-nano satellite constellation, comprising,

[0049] For a same-orbit micro-nano satellite constellation operating in a low earth orbit, a single-satellite atmospheric drag dynamics model is established to obtain the corresponding relationship between the control input and the control command of a single satellite;

[0050] A reference satellite is determined in the satellite constellation, the state constraint condition of other single satellites is determined based on the reference satellite, and then the simulated annealing algorithm is used for constellation slot assignment in combination with the energy constraint and the task constraint of the single satellite, and the constellation slot assignment optimization is performed with the shortest control time as the target to obtain the optimal result of the constellation slot assignment;

[0051] Based on the optimal result of the constellation slot assignment, an initial control command matrix is calculated, an adaptive function is set, the adaptive function value is lowest and the target control accuracy is met as the control target, and the initial control command matrix is optimized to obtain an optimal control command matrix;

[0052] According to the correspondence between the control input of a single satellite and the control command, the optimal satellite control input is obtained from the optimal control command matrix, and the angular acceleration of a single satellite in the constellation is controlled, so that angular acceleration differences are generated among different satellites in the constellation, thereby achieving phase adjustment of the same-orbit constellation.

[0053] This implementation scheme controls satellite attitude by designing a control command matrix, thereby implementing a phase adjustment control strategy. It provides energy and mission constraints for the constellation based on its actual on-orbit mission requirements. Based on the constellation's initial state vector, energy constraints, and imaging constraints, a simulated annealing algorithm is used to determine the target relative phase of each satellite in the constellation and allocate constellation slots. Based on the slot allocation results, a simulated annealing algorithm is then used to generate a command matrix for constellation phase adjustment.

[0054] To perform atmospheric differential drag control on a co-orbital constellation, we first establish a single satellite dynamic model under atmospheric differential drag control. According to Kepler's law, for a satellite moving at a constant speed on a circular orbit, the velocity v along the orbit in the inertial system is:

[0055]

[0056] Where μ E is the Earth's gravitational constant, R E is the radius of the Earth, h g is the orbital altitude. Atmospheric differential drag control controls the satellite through atmospheric dynamic pressure. In the exponential atmospheric model, the atmospheric density ρ can be expressed as:

[0057]

[0058] Where ρ0, H and h0 are constants in the model and are related to the orbital altitude. The atmospheric dynamic pressure q on the satellite is:

[0059]

[0060] For a mass m, the drag coefficient is C d For a satellite with a frontal area of ​​A, the satellite's ballistic coefficient is B. C =m / C d A. Use v x represents the satellite's velocity along the track. When the satellite moves in the atmosphere, the equation of motion along the track is:

[0061]

[0062] The above formula is written in differential form as follows:

[0063]

[0064] where R is the distance of the spacecraft from the center of gravity, in a circular orbit β0is the angle between the satellite velocity vector and the horizon, β0= 0 since the satellite is in a circular orbit; then we have:

[0065]

[0066] Next, to illustrate the satellite's dynamics model under atmospheric differential drag control, an ideal satellite on the same circular orbit as the satellite discussed is established, which always maintains uniform motion on the orbit and the orbit height is constant. Let v y be the normal velocity of the satellite in the inertial system, according to the characteristics of the ideal satellite, the normal relative velocity between the satellite and the ideal satellite is also v y . According to the Hill relative motion equation, the relative motion solution of a satellite under atmospheric differential drag control with respect to the ideal satellite in the orbital plane within a time step Δt is:

[0067]

[0068] Taking the satellite's entry point into the orbit as the phase zero point, let the satellite's absolute phase on the orbit be θ, then the angular acceleration is:

[0069]

[0070] From the above formula, it can be seen that by changing the satellite's windward area A, the ballistic coefficient B C can be changed, and then the satellite's absolute angular acceleration on the orbit

[0071] Further, the satellite's phase adjustment process is divided into several discrete time periods with a time step of Δt, and the satellite's windward area A remains unchanged in each time period, then the single-satellite atmospheric drag dynamics model in the along-track direction satisfies:

[0072]

[0073] where θ k is the satellite's absolute phase on the orbit in the discrete time period k, is the satellite's absolute angular velocity on the orbit in the discrete time period k, Δt is the control time step, which can be regarded as a small quantity, k represents the kth discrete time period, u k is the control command in the discrete time period k, u k ∈ [0, 1], B k is the control input in the discrete time period k;

[0074] According to formula (1), the corresponding relationship between the single-satellite control input and the control command is:

[0075]

[0076] wherein is the absolute angular acceleration of the satellite in the high drag state in the discrete time period k, is the absolute angular acceleration of the satellite in the low drag state in the discrete time period k, the high drag state is the satellite attitude when the satellite has the maximum windward area, the low drag state is the satellite attitude when the satellite has the minimum windward area, k = 0, 1, 2, 3, … end; wherein end represents the number of discrete time periods.

[0077] A strategy is designed to control the angular acceleration of the satellite by changing the satellite attitude and controlling the area-mass ratio of the satellite, so as to realize the phase adjustment of the satellite constellation. By controlling the angular acceleration of each satellite in the constellation, the angular acceleration difference of different satellites in the constellation is generated, so as to realize the phase adjustment of the satellite constellation. By uploading the control command to the satellite, the current attitude of the satellite is changed, and the value of the control command represents the area-mass ratio of the satellite in different attitudes.

[0078] The control command of each satellite in the constellation in each control time period during the phase adjustment process can form a command matrix U.

[0079] In the phase adjustment process of the satellite constellation:

[0080] 1) A common reference satellite is selected in the constellation, the rest of the satellites in the constellation are non-reference satellites, and the satellites in the constellation are numbered, wherein the reference satellite is No. 1;

[0081] 2) The initial relative motion state of each non-reference satellite relative to the reference satellite is determined;

[0082] 3) The target relative phase of each non-reference satellite relative to the reference satellite is determined.

[0083] The reference satellite is taken as No. 1 satellite in the satellite constellation, and the state constraint condition of the other single satellites is expressed as:

[0084]

[0085] wherein θ i,0 is the initial absolute phase of the i-th satellite, i = 2, 3, 4, …, N, wherein N is the number of satellites in the constellation; θ 1,0 is the initial absolute phase of the No. 1 satellite, θ i ′ ,0 is the initial relative phase of the i-th satellite relative to the No. 1 satellite, θ i,end is the final absolute phase of the i-th satellite after the control of the end discrete time periods, θ 1,end is the final absolute phase of the No. 1 satellite after the control of the end discrete time periods, θ i ′ ,fTarget relative phase of satellite i with respect to satellite 1.

[0086] The process of constellation phase control is composed of two parts, slot assignment and command generation, as shown in Figure 1 . The slot assignment part gives the slot assignment scheme that can make the phase adjustment time shortest, and then the command generation part generates the control command matrix for phase adjustment in the optimal time based on the assignment result.

[0087] Considering the actual on-orbit situation of the low-orbit satellite constellation, energy and task constraints are given.

[0088] When the satellite solar panel is oriented to the sun, the attitude is different from the high resistance state and the low resistance state, and the angular acceleration of the satellite motion is between and .

[0089] In this embodiment, when the satellite is in the light area and charging, the energy constraint of a single satellite is:

[0090] u k = α, 0 < α < 1,

[0091] wherein α represents the attitude command of the satellite in the charging window;

[0092] When the low-orbit satellite is on-orbit, it needs to perform on-orbit tasks, which usually have requirements for the satellite attitude. Taking the imaging task as an example, during the imaging task, the on-board camera is required to be oriented to the ground. Similar to charging, the attitude during the execution of the on-orbit task is also different from the high resistance state and the low resistance state, and the task constraint of a single satellite is:

[0093] u k = β, 0 < β < 1

[0094] wherein β represents the attitude command of the satellite in the task window;

[0095] According to the above analysis, the satellite needs a continuous charging time in the light area and a continuous on-orbit task time determined by the task. Taking the execution of remote sensing imaging or ground communication tasks in a specific latitude area as an example, without considering the precession of the orbital plane, the absolute phase of the task window in such tasks can be assumed to be determined. Set the satellite charging window phase interval Φ E and the satellite task window phase interval Φ M , so that the on-orbit absolute angular acceleration of the satellite is between the on-orbit absolute angular acceleration in the high resistance state and the on-orbit absolute angular acceleration in the low resistance state: The design of the constraint window depends on the actual task requirements of the constellation, as shown in Figure 2 When the satellite is in the charging window or the task window, the satellite attitude is subject to the corresponding constraints.

[0096] Satellite charging window phase interval Satellite mission window phase interval wherein is the minimum charging phase, is the maximum charging phase, is the minimum mission phase, is the maximum mission phase;

[0097] The energy and mission constraints affect the phasing control effect of constellation atmospheric differential drag control because the angular acceleration of the satellite in the high and low resistance states within the constraint window is smaller than the relative angular acceleration of other satellites in the constellation, which prolongs the time required for phase adjustment.

[0098] The control command u of satellite i at discrete time k under multi-resistance state phase control is obtained according to the different resistance states of the satellite ik is:

[0099]

[0100] In a set of low-orbit micro-nano satellite constellation in the same orbit, a satellite is confirmed as a reference satellite, and the phases of other satellites are considered as relative phases with the reference satellite. Each non-reference satellite forms a satellite pair with the reference satellite.

[0101] Further, the method of using simulated annealing algorithm for constellation slot allocation is:

[0102] Set the initial annealing temperature T0, the lower limit of annealing temperature T min , the annealing coefficient ξ, and generate the initial slot allocation vector S0:

[0103] S0 = [1 2... N],

[0104] The constellation numbers in S0 are arranged in correspondence with the target relative phase numbers; N is the number of satellites in the constellation.

[0105] The initial slot allocation vector S0 is optimized to obtain the optimal result of constellation slot allocation.

[0106] The physical meaning of the slot allocation vector is that if the m0th bit in S0 is n0, then the target relative phase of the n0th satellite with respect to the reference satellite is m0.

[0107] Slot refers to the corresponding relationship between the satellite and the target relative phase, and slot allocation refers to the corresponding order of the satellites in the constellation and the target phases of the constellation. The time required for phase adjustment of a satellite pair is related to the initial relative state and the target relative phase of the satellite pair. Therefore, in order to shorten the total control time of constellation phase adjustment, it is necessary to optimize the slot allocation to obtain the slot allocation with the shortest control time. Before optimizing the slot allocation, an initial slot allocation vector must be generated.

[0108] Randomly generate the initial slot allocation S0, that is, randomly number the satellites in the constellation, and the No. 1 satellite is the reference satellite. The constellation number is arranged in one-to-one correspondence with the target relative phase number.

[0109] In this embodiment, the method for calculating the phase adjustment control time according to the initial slot allocation vector S0 is:

[0110] Each non-reference satellite is combined with the reference satellite to form a satellite pair;

[0111] Initialize the current annealing temperature T = T0, and the current old slot allocation vector S old = S0;

[0112] Calculate the old slot allocation vector S old In the case of each satellite pair, the phase adjustment control time of the satellite pair, for each satellite pair, the relative angular acceleration of the satellite pair is

[0113]

[0114] In the formula, the subscript sat corresponds to the non-reference satellite in the satellite pair, the subscript ref corresponds to the reference satellite in the satellite pair, is the on-orbit absolute angular acceleration of the non-reference satellite in the satellite pair, is the on-orbit absolute angular acceleration of the reference satellite in the satellite pair; q is the atmospheric dynamic pressure acting on the satellite, is the orbit average semi-major axis, B C is the ballistic coefficient of the satellite,

[0115] In the satellite pair, the change in orbital height caused by the differential atmospheric drag control is a small amount compared to the orbital height. It can be assumed that the orbital height of the satellite does not change during the phase adjustment process, the atmospheric dynamic pressure acting on the satellite is constant, the orbit average semi-major axis of the reference satellite and the non-reference satellite is the same, and formula (5) is simplified as:

[0116]

[0117] Calculate the shortest control time of the satellite pair phase adjustment:

[0118] The change in orbital height during the phase adjustment process is very small, and the satellite motion period changes very little. It can be considered that the satellite maintains a constant motion period. Therefore, based on the specific constraint conditions designed, it is considered that the phase adjustment total time of the phase proportion of the constraint window on the orbit and the running time of the satellite in the constraint window is approximately equal.

[0119] ​In order to simplify the calculation of the shortest time of phase adjustment of the satellite pair, the phase adjustment process is simplified as the following extreme case: when the reference satellite is in the high resistance state, the non-reference satellite is in the low resistance state; when the reference satellite is in the low resistance state, the non-reference satellite is in the high resistance state; and the two satellites in the satellite pair are simultaneously located in the same type of constraint window. In this extreme case, the phase adjustment process of the satellite pair is least affected by constraints, the controllability of the satellite pair is highest, and it is the shortest time of phase adjustment of the satellite pair. On the basis of simplification, the shortest time of phase adjustment of the satellite pair is calculated, and the control process is divided into three sections A, B and C according to time; in section A, the resistance states of the reference satellite and the non-reference satellite are opposite; in section B, the resistance states of the reference satellite and the non-reference satellite are opposite to those in section A; in section C, the reference satellite and the non-reference satellite are located in the constraint window, and the relative phase is unchanged; the constraint window includes the satellite charging window and the satellite task window; the running time t C of the satellite pair in section C is:

[0120]

[0121] wherein t A is the running time of the satellite pair in section A, t B is the running time of the satellite pair in section B;

[0122] The initial relative phase of the satellite pair is set as θ0′, and the initial relative angular velocity of the satellite pair is The target relative phase of the satellite pair is θ f ′, and the target relative angular velocity of the satellite pair is The relative angular acceleration of the satellite pair in section A is expressed as The relative angular acceleration of the satellite pair in section B is expressed as

[0123] If the non-reference satellite is in the high resistance state and the reference satellite is in the low resistance state in section A, then:

[0124]

[0125] wherein is the on-orbit absolute angular acceleration of the satellite in the high resistance state in the satellite pair, is the on-orbit absolute angular acceleration of the satellite in the low resistance state in the satellite pair;

[0126] If the non-reference satellite is in the low resistance state and the reference satellite is in the high resistance state in section A, then:

[0127]

[0128] According to the kinematic principle, it is obtained that:

[0129]

[0130] wherein Δθ′=θ f ′-θ0′,

[0131] Substitute and into formula (6) and (7) respectively to solve, get four groups of t A and t B solution, select only one group of positive real number solution as t A and t B The final solution of, and determine the real state of the satellite pair in section A and section B;

[0132] Thus the phase adjustment control time t is obtained:

[0133] t=t A +t B +t C .

[0134] Further, the method for optimizing the initial slot allocation vector S0 is:

[0135] Calculate the fitness function f(S old ) corresponding to the current old slot allocation vector S old , f(S old ) takes the shortest control time as the optimization goal, and the fitness function f(S old ) is the constellation phase adjustment control time under the current old slot allocation vector S old .

[0136] Assume that the correspondence between the current old slot allocation vector S old and the slot number is:

[0137] S old =[s1 s2 s3...s N ],

[0138] In the formula, s i (i=1,2,...,N) is the slot number, each slot number s i corresponds to a constellation number, assuming that the constellation number corresponding to the slot number s i is I, then the target relative phase of the Ith constellation in the slot s i is I, I=2,3,4,...,N; wherein the slot s1 corresponds to the first satellite, and the target relative phase is 0, s1≡1.

[0139] The constellation phase adjustment time takes the maximum value of the control time of each satellite pair in the constellation, so the fitness function f(S old ) is:

[0140] f(S old )=max{t(s2),t(s3),...,t(s N},

[0141] where t(s N ) is the slot s N corresponding to the phase adjustment control time of the satellite pair composed of the satellite and the reference satellite; the smaller the value of f(s old ), the shorter the time required for constellation phase adjustment, and the better the solution adaptability.

[0142] Iterative optimization is performed according to the following process:

[0143] Randomly exchange two slot positions in the current old slot allocation vector S old , assuming that the slot and the slot are exchanged, then the position relationship after the exchange is:

[0144]

[0145] where m1 and m2 are two randomly generated exchange slot numbers;

[0146] The new slot allocation vector S new obtained after the exchange is taken as a new solution generated by disturbance; the fitness function f(s new ) is calculated.

[0147] According to the Metropolis principle, the acceptance probability P of the fitness function f(s new ) is calculated:

[0148]

[0149] If the new slot allocation vector S new is accepted, S old =S new and f(s old )=f(s new ) are performed, then annealing processing is performed, and T decreases with the iteration number. The current annealing temperature T is reduced to λT, λ∈(0, 1), λ is a temperature adjustment coefficient, if T>T min , then the slot position exchange step is returned, and the next round of iterative calculation is performed, until T≤T min , the loop is ended.

[0150] The final obtained slot allocation vector S new is taken as the optimal constellation slot allocation result S best , and the final obtained phase adjustment control time t is taken as the shortest control time t best :

[0151] t best =f(s best ).

[0152] The above process can be expressed by pseudo code as follows:

[0153]

[0154] Based on the obtained optimal slot allocation result, the approximate ratio of high-resistance state time and low-resistance state time of each satellite in the constellation is obtained to obtain a better command matrix as an initial solution.

[0155] In the embodiment, the method for calculating the initial control command matrix is as follows:

[0156] For the non-reference satellite, the approximate ratio p of the high-resistance state running time and the low-resistance state running time of the satellite pair is calculated i :

[0157]

[0158] The approximate ratio p1 of the high-resistance state running time and the low-resistance state running time of the reference satellite in the satellite pair is:

[0159]

[0160] Determine the control time step Δt;

[0161] Since the satellite attitude is limited by the constraint condition, the time when the satellite is in the constraint window cannot be determined, and the concept of high-low resistance matrix U HL is introduced. Ideally, U HL is the high-resistance state and low-resistance state part in the command matrix U, but because the high-resistance state and low-resistance state part of each satellite in the command matrix is not necessarily equal in number of columns, there may be a difference in the right end.

[0162] First, based on the optimal results of constellation slot allocation and the calculation results of the approximate ratio, the initial high-low resistance matrix U HL0 is established. The initial high-low resistance matrix U HL0 has N rows and M columns, and:

[0163]

[0164] In the formula, v 0N is the control command of the Nth satellite;

[0165] The elements in the high-low resistance matrix U HL will be used as the elements in the high-low resistance state of the command matrix in turn. When the satellite is located in the charging window and the task window, the control instruction in the command matrix is determined by formula (4), so the optimization of the command matrix can be regarded as the optimization of the high-resistance state and low-resistance state part, that is, the optimization of U HL .

[0166] The initial high-low resistance matrix U HL0controlling the constellation, when the satellite is in the constraint window, in U HL0 supplementing the control command calculated according to formula (4) with the corresponding position in U best to obtain an initial control command matrix U0.

[0167] The method for obtaining the current control command matrix comprises:

[0168] setting the current high-low resistance matrix U HL = U HL0 , and the current loop number n = 1, and setting the maximum loop number as n max .

[0169] initializing a loop variable, setting the time t0 = 0, j = 1, i = 1, wherein j is the control command number of the i-th satellite;

[0170] judging whether the i-th satellite in the constellation is located in the charging window or the task window in the current discrete time period;

[0171] The iterative calculation process for obtaining the current control command matrix comprises:

[0172] if the current satellite is located in the satellite charging window, u k = a; if the current satellite is located in the satellite task window, u k = b;

[0173] if the current satellite is not in the constraint window, i.e. neither in the charging window nor in the task window, taking the j-th control command of the current satellite in the current high-low resistance matrix U HL as the control command of the i-th satellite, and then setting j = j + 1; when j exceeds the column width of the current high-low resistance matrix U HL , i.e. exceeds the dimension, making the corresponding control command randomly 0 or 1 to obtain the current control command matrix U;

[0174] controlling the i-th satellite according to the corresponding control command in the current control command matrix U according to formula (4) and formula (5), and updating t k+1 = t k + Δt;

[0175] judging whether the current t k+1 is less than f(S best ); if yes, returning to perform the iterative calculation of the next control command matrix; until the control of the current satellite is completed; otherwise, recording the absolute phase and absolute angular velocity of the current satellite at the last time, and updating the loop variable, t0 = 0, j = 1, i = i + 1, and continuing to control the next satellite until the control of all satellites in the satellite constellation is completed.

[0176] Finally, the method to obtain the optimal control command matrix is:

[0177] Perform fitness function calculation and define fitness function g(U):

[0178]

[0179] The control effect of the command matrix is ​​judged by the function value of the adaptability function. The lower the adaptability function g(U), the better the control effect obtained in the shortest control time and the better the control command matrix.

[0180] The optimal control command matrix is ​​obtained by iterating iteratively as follows:

[0181] Randomly select the current high and low resistance matrix U HL An element in , changes the original control command u of the current element old , get the new control command u of the current element new :

[0182] u new =1-u old ,

[0183] All new control commands u new Get the new high and low resistance matrix U HL,new ; where u old For U HL Any element in u new For the new high and low resistance matrix U HL,new Neutralize old Elements at the same position.

[0184] Calculate the new high and low resistance matrix U according to the method of obtaining the current control command matrix HL,new The corresponding new control command matrix Unew and the corresponding adaptability function g(U new );

[0185] According to the Metropolis principle, the acceptance of the fitness function g(U new ) probability W:

[0186]

[0187] If accepted, U HL =U HL,new , U=U new , g(U)=g(U new ); and determine whether the current satellite has reached the target relative phase, set the preset accuracy of the target relative phase to h, and determine whether the satellite has reached the target phase:

[0188] (θi,end-θ1,end)-θi′,f|<h(8),

[0189] If all satellite pairs satisfy formula (8), it is considered that the control accuracy condition is satisfied, the loop can be ended, and the current control command matrix U is taken as the optimal control command matrix;

[0190] If all satellite pairs do not satisfy formula (8), but the current has reached the maximum loop number n max , the current control command matrix U is taken as the optimal control command matrix, and the loop is combined; if the current loop number has not reached the maximum loop number n max , the loop iteration of the optimal control command matrix is returned to continue, until the ending condition is reached.

[0191] Finally, the optimal command matrix is simulated to obtain the simulation results.

[0192] The specific implementation process can be represented by pseudo code as follows:

[0193]

[0194]

[0195] Next, a satellite constellation model composed of ten satellites is taken as an example to simulate the method. The constellation simultaneously enters a near-earth circular orbit at an altitude of 440 km, and is expected to achieve on-orbit phase uniformity. Under the energy and imaging constraint conditions, Φ M = [50° 100°], Φ E = [0° 50°], α = 0.4, and β = 0.6. Taking the satellite with the smallest absolute angular velocity as the reference satellite, the initial relative angular velocities between each satellite and the reference satellite caused by the on-orbit impulse are shown in Table 1, and the maximum relative angular velocity is 5.501 deg / day.

[0196] Table 1 Initial relative angular velocity of simulation constellation

[0197]

[0198] According to the empirical exponential atmospheric density model, when alt = 440 km, ρ0 = 3.725 × 10 -12 kg / m 3 , h0 = 400 km, and H = 59.4 km, the atmospheric density is calculated as ρ = 1.900 × 10 -12 kg / m 3 . The ballistic coefficients in the low resistance state and the high resistance state are respectively taken as There are

[0199] First, an unoptimized initial slot allocation Y0 is randomly generated. Except for the reference satellite, the target slots for the remaining satellites are arranged in the default order used for relative state estimation, i.e., S0 = [1 2 ... 9 10]. The target slot numbering is shown in Table 2 below.

[0200] Table 2 Target slot numbers

[0201]

[0202] Under the initial slot allocation method, the estimated constellation control time is 44.399 days. The number of iterations is set to 314 times. After slot allocation optimization, the optimal slot allocation method is S best =[18107 6539 24], the corresponding constellation control time estimate is 39.647 days. Compared with the initial slot allocation method, the control time estimate is shortened by 4.752 days. Figure 3 The variation of the estimated control time during the annealing optimization process is shown, and the estimated control time eventually converges to the minimum value. Figure 4 Phase modulation estimation with two different slot allocations is shown, and the impact of slot allocation on control time estimation can be seen.

[0203] The satellites in the constellation are renumbered based on the slot allocation results, and an initial high- and low-resistance matrix is ​​generated. Due to the low altitude and short period of the satellite orbit, a relatively short discrete control step size is required to ensure that the control conditions meet the actual control requirements. Consider a control time step of Δt = 0.004 days.

[0204] The result of reference state estimation randomly generates the initial high and low resistance matrix U HL0 , the corresponding fitness function value g(U HL0 )=10328deg 2 , does not meet the control requirements. In the simulated annealing optimization process, the number of cycles is set to 2×10 6 times, such as Figure 5 As shown in the figure, the optimal high and low resistance matrix U is obtained after optimization HL,best and the corresponding optimal command matrix U best , and obtain the corresponding optimal value of fitness function g(U best )=1×10- 4 deg 2 . After inspection, it meets the control accuracy requirements.

[0205] The simulation results show that, under the control of the optimal command matrix, the constellation can complete the phase adjustment within the control time estimation value of the slot allocation, meet the control accuracy requirement of equation (8), realize the phase uniform distribution of the constellation, and continuously maintain the phase under the control of the optimal control matrix. After 39.980 days of phase adjustment control, the constellation reaches the target phase within 1 deg, and is stably maintained at the target phase for 0.067 days thereafter.

[0206] Figure 6 For the simulation results of the relative phase of the in-orbit constellation with respect to time at the end of the control time, it can be seen that, during the operation of the constellation, based on the atmospheric differential drag control, the multi-resistance state phase control method is used, and the controlled constellation realizes the control target of in-orbit phase uniform distribution. The absolute phase of each satellite in the constellation at the end of the control time and the relative phase with respect to the reference satellite are shown in Table 3, and the difference between the relative phase and the target relative phase is within 1 deg, meeting the control requirements.

[0207] Table 3 End state of optimal control phase adjustment considering constraints

[0208]

[0209]

[0210] In summary, the method of the application is used for phase deployment of phase uniform distribution of in-orbit low-orbit micro-nano satellite constellation, considers the actual in-orbit operation requirements of the constellation, and provides a method for realizing the phase control of the constellation by using atmospheric differential force drag force. In practice, it can be applied to the orbit deployment of low-orbit spacecraft constellation according to the specific tasks and mass ratio characteristics of the satellite.

[0211] Although the application has been described herein with reference to particular embodiments, it will be understood that the examples are merely examples of the principles and applications of the application. It should be understood that many modifications can be made to the exemplary embodiments, and other arrangements can be designed without departing from the spirit and scope of the application as defined in the appended claims. It should be understood that different dependent claims and features described herein can be combined in ways other than those described. It can also be understood that features described in connection with individual embodiments can be used in other described embodiments.

Claims

1. A method for controlling atmospheric differential drag phase adjustment for a low-orbit micro-nano satellite constellation, characterized in that: For a constellation of micro- and nano-satellites operating in low Earth orbit, a single-satellite atmospheric drag dynamics model is established to obtain the corresponding relationship between the control input and control command of a single satellite. A reference satellite is identified in the satellite constellation. Based on the reference satellite, the state constraints of other single satellites are determined. Then, a simulated annealing algorithm is used to allocate constellation slots, combining the energy constraints and mission constraints of the single satellite. The constellation slot allocation is optimized with the goal of minimizing the control time to obtain the optimal result. The initial control command matrix is ​​calculated based on the optimal result of constellation slot allocation. Then, an adaptability function is set, and the initial control command matrix is ​​optimized with the lowest adaptability function value and meeting the target control accuracy as the control goal to obtain the optimal control command matrix. Based on the correspondence between single-satellite control input and control commands, the optimal satellite control input is obtained from the optimal control command matrix, and the angular acceleration of a single satellite in the constellation is controlled, so that different satellites in the constellation have angular acceleration differences, thereby achieving phase adjustment of the same-orbit constellation. The atmospheric drag dynamics model for a single satellite is: Where θ k is the absolute phase of the satellite in orbit in discrete time period k, is the absolute angular velocity of the satellite in the discrete time period k, Δt is the control time step, k represents the kth discrete time, u k is the control command in discrete time period k, B k is the control input in discrete time period k; According to formula (1), the corresponding relationship between the single satellite control input and the control command is: In the formula is the absolute angular acceleration of the satellite in high-impedance state in discrete time period k, is the absolute angular acceleration of the satellite in the low-resistance state in discrete time period k. The high-resistance state is the satellite attitude when the satellite's frontal area is maximum, and the low-resistance state is the satellite attitude when the satellite's frontal area is minimum, k = 0, 1, 2, 3, ..., end; where end represents the number of discrete time periods.

2. The atmospheric differential drag phase adjustment control method for a low-orbit micro-nano satellite constellation according to claim 1 is characterized in that: Taking the reference satellite as satellite No. 1 in the satellite constellation, the state constraints of other single satellites are expressed as: Where θ i,0 is the initial absolute phase of satellite i, i = 2, 3, 4, ..., N, where N is the number of satellites in the constellation; θ 1,0 is the initial absolute phase of satellite No. 1, θ′ i,0 is the initial relative phase of satellite i relative to satellite 1, θ i,end is the final absolute phase of satellite i after the discrete time control of the end segment, θ 1,end is the final absolute phase reached by satellite 1 after the discrete time control of the end segment, θ′ i,f is the target relative phase of satellite i relative to satellite 1.

3. The atmospheric differential drag phase adjustment control method for a low-orbit micro-nano satellite constellation according to claim 2 is characterized in that: The energy constraint of a single satellite is: u k ＝α,0＜α＜1， Where α represents the attitude command when the satellite is within the charging window; The mission constraints of a single satellite are: u k ＝β,0＜β＜1 Where β represents the attitude command of the satellite within the mission window; Set the satellite charging window phase interval Φ E and the satellite mission window phase interval Φ M , so that the satellite's on-orbit absolute angular acceleration is between the high-resistance state on-orbit absolute angular acceleration and the low-resistance state on-orbit absolute angular acceleration: Satellite charging window phase interval Satellite mission window phase interval in is the minimum charging phase, is the maximum charging phase, is the minimum mission phase, is the maximum mission phase; Get the control command u of satellite i in discrete time period k ik for:

4. The atmospheric differential drag phase adjustment control method for a low-orbit micro-nano satellite constellation according to claim 3 is characterized in that: The method for allocating constellation slots using the simulated annealing algorithm is: Set the annealing initial temperature T0 and the annealing temperature lower limit T min , annealing coefficient ξ, generate the initial slot allocation vector S0: S0=[12...N], The constellation numbers in S0 are arranged corresponding to the target relative phase numbers; With the goal of minimizing the control time, the initial slot allocation vector S0 is optimized to obtain the optimal result of constellation slot allocation.

5. The atmospheric differential drag phase adjustment control method for a low-orbit micro-nano satellite constellation according to claim 4 is characterized in that: The method for calculating the phase adjustment control time according to the initial slot allocation vector S0 is: Each non-reference satellite is combined with a reference satellite to form a satellite pair; Initialize the current annealing temperature T = T0, the current old slot allocation vector S old =S0; For each satellite pair, the relative angular acceleration of the satellite pair for: Where the subscript sat corresponds to the non-reference satellite in the satellite pair, and ref corresponds to the reference satellite in the satellite pair. is the absolute angular acceleration of the satellite relative to the non-reference satellite, is the absolute angular acceleration of the reference satellite in orbit; q is the atmospheric dynamic pressure on the satellite, is the mean semi-major axis of the orbit, B C is the satellite ballistic coefficient, Assuming that the satellite orbit height remains unchanged during the phase adjustment process, the atmospheric dynamic pressure is constant, and the average semi-major axis of the reference satellite and the non-reference satellite is the same, then formula (5) is simplified to: Calculate the shortest control time of satellite phase adjustment: The control process is divided into three sections, A, B, and C, according to time. In section A, the resistance states of the reference satellite and the non-reference satellite are opposite. In section B, the resistance states of the reference satellite and the non-reference satellite are opposite to their resistance states in section A. In section C, both the reference satellite and the non-reference satellite are located in the constraint window, and the relative phase remains unchanged. The constraint window includes the satellite charging window and the satellite mission window. The operation time of the satellite pair in section C is t C for: Where t A is the A-segment operation time of the satellite pair, t B is the B-segment operation time of the satellite pair; Assume that the initial relative phase of the satellite pair is θ0′ and the initial relative angular velocity of the satellite pair is The relative phase of the satellite to the target is θ f ′, the relative angular velocity of the satellite to the target is The relative angular acceleration of the satellite pair in segment A is expressed as The relative angular acceleration of the satellite pair in segment B is expressed as If the non-reference satellite in segment A is in high impedance state and the reference satellite is in low impedance state, then: In the formula is the absolute angular acceleration of the satellite on orbit relative to the medium and high resistance satellite, is the absolute angular acceleration of the satellite on orbit relative to the low-resistance satellite; If the non-reference satellite in segment A is in low impedance state and the reference satellite is in high impedance state, then: According to the principle of kinematics: In the equation Δθ′=θ f ′-θ0′, Will and Substitute into formula (6) and (7) respectively to solve and obtain four groups of t A and t B The solution of t is to select a set of positive real number solutions as t A and t B The final solution of , and determine the true resistance state of the satellite pairs in segments A and B; The phase adjustment control time t is obtained as follows: t=t A +t B +t C 。 6. The atmospheric differential drag phase adjustment control method for a low-orbit micro-nano satellite constellation according to claim 5 is characterized in that: The method for optimizing the initial slot allocation vector S0 is: Calculate the current old slot allocation vector S old The corresponding fitness function f(S old ), f(S old ) takes the shortest control time as the optimization goal, and the fitness function f(S old ) allocates vector S for the current old slot old The constellation phase under allocation adjusts the control time; Assume that the current old slot allocation vector S old The corresponding relationship with the slot number is: S old =[s1 s2 s3...s N ], Where s i (i=1,2,...,N) is the slot number, each slot number s i Corresponding to a constellation number, assuming the slot number s i The corresponding constellation number is I, then slot s i The target relative phase of constellation I is phase I, where I = 2, 3, 4, ..., N; slot s1 corresponds to satellite 1, and its target relative phase is 0, s1≡1; Then the fitness function f(S old )for: f(S old )=max{t(s2),t(s3),...,t(s N )}, Where t(s N ) is slot s N Phase adjustment control time of a satellite pair consisting of a corresponding satellite and a reference satellite; Perform iterative optimization as follows: Randomly swap the current old slot allocation vector S old In the two slot positions, assuming that the slot and slots If the positions are exchanged, the position relationship after the exchange is: Where m1 and m2 are two randomly generated switch slot numbers; After the exchange, the new slot allocation vector S is obtained new , calculate the fitness function f(S new ); According to the Metropolis principle, calculate the fitness function f(S new )’s acceptance probability P: If a new slot allocation vector S is accepted new , then S old =S new ,f(S old )=f(S new ), annealing is performed to reduce the current annealing temperature to T = λT, λ∈(0,1), where λ is the temperature adjustment coefficient. If T>T min , then return to the slot position exchange step and perform the next round of iterative calculation until T≤T min , end the loop; The resulting slot allocation vector S new As the optimal result of constellation slot allocation S best , the final phase adjustment control time t is taken as the shortest control time t best : t best =f(S best )。 7. The atmospheric differential drag phase adjustment control method for a low-orbit micro-nano satellite constellation according to claim 6, characterized in that: The method for calculating the initial control command matrix is: Calculate the approximate ratio p of the high-impedance operation time and low-impedance operation time of non-reference satellite i in the satellite pair i : Then the approximate ratio p1 of the high-impedance state operation time and the low-impedance state operation time of the satellite centering reference satellite is: Based on the optimal result of constellation slot allocation and the calculation results of approximate ratio, the initial high and low resistance matrix U is established HL0 , initial high and low resistance matrix U HL0 There are N rows and M columns, and: Where v 0N Control command for satellite N; Using the initial high and low resistance matrix U HL0 Control the constellation. When the satellite is within the constraint window, HL0 The corresponding position in the equation (4) is supplemented with the control command calculated, combined with the shortest control time f(S best ), and obtain the initial control command matrix U0.

8. The atmospheric differential drag phase adjustment control method for a low-orbit micro-nano satellite constellation according to claim 7, characterized in that: Methods for obtaining the current control command matrix include: Make the current high and low resistance matrix U HL =U HL0 , and the current number of loops n = 1, set the maximum number of loops to n max ; Initialize loop variables so that time t0 = 0, j = 1, i = 1, where j is the control command number of satellite i; The iterative calculation process to obtain the current control command matrix is: If the current satellite is within the satellite charging window, u k =α; if the current satellite is within the satellite mission window, then u k =β; If the current satellite is not within the constraint window, the current high and low resistance matrix U HL The jth control command of the current satellite is used as the control command of satellite i, and then j=j+1; when j exceeds the current high and low resistance matrix U HL The column width is set to 0 or 1 randomly, and the current control command matrix U is obtained; According to formula (4) and formula (5), the corresponding control command in the current control command matrix U is used to control satellite i, and t is updated. k+1 =t k +Δt; Determine the current t k+1 Is it less than f(S best ), if so, return to perform the next round of iterative calculation of the control command matrix; until the control of the current satellite is completed; otherwise, record the absolute phase and absolute angular velocity of the current satellite at the last moment, and update the loop variables, t0=0, j=1, i=i+1, and continue to control the next satellite until the control of all satellites in the satellite constellation is completed.

9. The atmospheric differential drag phase adjustment control method for a low-orbit micro-nano satellite constellation according to claim 8, characterized in that: The method to obtain the optimal control command matrix is: Define the fitness function g(U): The lower the adaptability function g(U), the better the control effect obtained in the shortest control time and the better the control command matrix; The optimal control command matrix is ​​obtained by iterating iteratively as follows: Randomly select the current high and low resistance matrix U HL An element in , changes the original control command u of the current element old , get the new control command u of the current element new : in new =1-in old , All new control commands u new Get the new high and low resistance matrix U HL,new ; Calculate the new high and low resistance matrix U according to the method of obtaining the current control command matrix HL,new The corresponding new control command matrix U new and the corresponding fitness function g(U new ); According to the Metropolis principle, the acceptance of the fitness function g(U new ) probability W: If accepted, U HL =U HL,new , U=U new , g(U)=g(U new ); and determine whether the current satellite has reached the target relative phase, and set the preset accuracy of the target relative phase to h: |(θ i,end -θ 1,end )-θ′ i,f |<h (8), If all satellite pairs satisfy formula (8), the current control command matrix U is taken as the optimal control command matrix; If all satellite pairs do not satisfy formula (8), but the maximum number of cycles n has been reached max , then the current control command matrix U is used as the optimal control command matrix; if the current number of cycles does not reach the maximum number of cycles n max , then return to continue the loop iteration of the optimal control command matrix until the end condition is reached.

Images