A method for modeling nash equilibrium of on-orbit servicing of impulse thrust spacecraft under illumination constraint

Abstract

The application discloses a kind of illumination constraints under the Nash equilibrium modeling method of impulse thrust spacecraft on-orbit servicing, solve the problems of theoretical fault, illumination constraint modeling missing and the insufficient coupling of multi-physical field constraints in impulse thrust spacecraft on-orbit servicing modeling, belong to the field of spacecraft illumination constraint optimization;Including: setting the target scene corresponding to task demand;The dynamics model under target scene is composed of illumination model and relative motion state equation;Zero-sum game model is jointly constituted by Hamilton function, co-state equation and cross-section condition;The Nash equilibrium model is composed of the dynamics model under target scene and zero-sum game model, and the Nash equilibrium solution of Nash equilibrium model is found by using adaptive differential evolution algorithm;The application improves the efficiency and success rate of task planning.

Description

Technical Field

[0001] This invention belongs to the field of spacecraft illumination constraint optimization technology, and relates to a Nash equilibrium modeling method for on-orbit servicing of pulse thrust spacecraft under illumination constraints. Background Technology

[0002] With the large-scale development of space activities, the field of on-orbit servicing faces increasingly complex dynamic decision-making challenges. When both servicing spacecraft have maneuverability, traditional classical control theory, nonlinear control theory, and unilateral optimal control theory have significant limitations: classical control theory is difficult to handle dynamic interactions among multiple agents, nonlinear control theory lacks adaptability in decision-making under conflicting objectives, and unilateral optimal control theory cannot balance conflicting interests among multiple parties.

[0003] Differential game theory, by constructing an asymmetric optimization framework, offers a new approach to characterizing complex interactive on-orbit servicing scenarios (such as debris removal and satellite repair). While continuous thrust models can be described by differential equations in existing research, key technical bottlenecks remain in the on-orbit servicing of pulsed thrust spacecraft. The discrete control characteristics of pulsed thrust (such as short pulse duration and adherence to Kepler's laws between pulses) are fundamentally different from the continuous-time dynamics framework of traditional differential game theory, resulting in an underdeveloped theoretical framework and solution methods in this field. Furthermore, most existing studies neglect the impact of illumination conditions on the mission, while visible light imaging quality depends on illumination constraints, directly affecting the efficiency and success rate of on-orbit servicing missions.

[0004] In summary, the existing technology has the following technical problems:

[0005] The theoretical framework of the pulse thrust model is incomplete: traditional differential game theory relies on continuous-time dynamics and continuous control inputs, while the discrete control characteristics of pulse maneuvers (such as sudden changes in velocity increments and Kepler's law constraints during free motion) prevent the direct application of existing theories. Furthermore, the optimal strategy solution for continuous systems with superimposed discrete control has not yet formed a complete system, and is mainly solved through intelligent optimization algorithms. For example, in the invention publication CN117610255A, application number 202311551557.X, entitled "A Strategy Solution Method for Track Chase Game Based on Dynamic Game Tree Search," extensive research was conducted on pulse game tasks, a task model was established, and a minimax search method was used to solve the strategies of both sides. However, because a heuristic algorithm was used, the solution obtained was a feasible solution, not a unique one.

[0006] Illumination constraints are often overlooked: Existing research on on-orbit servicing mission planning and control strategy optimization has not fully considered the impact of illumination conditions on visible light imaging quality, and lacks a systematic analysis of mission feasibility under illumination constraints. This can lead to situations where effective observations cannot be achieved in actual missions due to insufficient illumination. However, it is necessary to use visible light cameras to obtain prior information about the target when carrying out on-orbit servicing missions, but current research on on-orbit servicing problems under illumination constraints is limited. Summary of the Invention

[0007] To address the theoretical gaps (such as the conflict between discrete control and continuous dynamics frameworks), lack of illumination-constrained modeling, and insufficient multiphysics constraint coupling in existing on-orbit servicing modeling of pulsed thrust spacecraft, this invention proposes a Nash equilibrium modeling method for on-orbit servicing of pulsed thrust spacecraft under illumination constraints. This method constructs a dynamic model that integrates the target scenario and a zero-sum game model, achieving systematic modeling of orbital dynamics, imaging feasibility, and the interaction of strategies between the two parties during pulsed maneuvers. It provides a model framework that combines theoretical rigor and engineering applicability for solving the optimal decision-making problem of on-orbit servicing of pulsed spacecraft under illumination constraints, thereby improving the efficiency and success rate of mission planning.

[0008] The objective of this invention is specifically achieved through the following technical solutions:

[0009] This invention discloses a Nash equilibrium modeling method for on-orbit servicing of pulse thrust spacecraft under illumination constraints, the method comprising:

[0010] Step 1: Set the target scenario corresponding to the mission requirements; Under the target scenario, construct a solar position model based on the ecliptic inclination and the average longitude of the sun, and construct an illumination model based on the solar position model; Simultaneously, in the virtual dynamic point LVLH coordinate system, the distance between the virtual origin and the spacecraft is much smaller than the distance to the Earth's center, and only considering the J2 term perturbation of the Earth's non-spherical shape, construct the relative motion state equation through the perturbation acceleration vectors of the service spacecraft and the target spacecraft; The illumination model and the relative motion state equation constitute the dynamic model under the target scenario;

[0011] Step 2: Construct the objective function using the position vectors of the service spacecraft and the target spacecraft at the terminal moment; construct the constraint conditions based on the orbit change position constraint, velocity increment constraint, time constraint, and terminal illumination constraint; construct Hamiltonian functions for the service spacecraft and the target spacecraft respectively; construct the costate equation using the partial derivatives of the Hamiltonian function with respect to the state vector; obtain the transverse condition based on the objective function, constraint conditions, and interior point function; and construct a zero-sum game model using the Hamiltonian function, costate equation, and transverse condition.

[0012] Step 3: The Nash equilibrium model is composed of the dynamic model and the zero-sum game model under the target scenario, and the Nash equilibrium solution of the Nash equilibrium model is found by using the adaptive differential evolution algorithm.

[0013] In step one, the target scenario corresponding to the task requirements is set as follows:

[0014] Both the service spacecraft and the target spacecraft achieve velocity increments through instantaneous pulses. The pulse timing and duration are the same, and each velocity increment satisfies the velocity increment constraint and the interval satisfies the time constraint.

[0015] Under the constraints of terminal illumination, the camera automatically points to the center of the target, satisfying the observation alignment conditions.

[0016] In step one, the solar position model is as follows:

[0017] ;in,

[0018] ;

[0019] ;

[0020] ;

[0021] In the formula, This provides information about the Sun's position in the geocentric inertial frame of reference. For the terminal time, The ecliptic inclination, d represents the average longitude of the sun; d represents the cumulative number of days, calculated from 12:00:00 on January 1, 2000. It is the angle closest to the point.

[0022] In step one, the lighting model is as follows:

[0023] ;in,

[0024] ;

[0025] In the formula, The illumination angle at the terminal moment. This represents the relative position vector between the target spacecraft and the sun at the terminal moment. This represents the relative position vector between the target spacecraft and the service spacecraft at the terminal moment; denoted as the position vector of the target spacecraft at the terminal moment.

[0026] In step one, the equation of state for relative motion is:

[0027] ;in,

[0028] ;

[0029] In the formula, Let be the acceleration vector of spacecraft i along the x, y, and z axes in the virtual moving point LVLH coordinate system, where i = P represents the serving spacecraft, i = E represents the target spacecraft, and T is the transpose; The average angular velocity of the orbit. , Here, is the gravitational parameter of the central celestial body, and 'a' is the semi-major axis of the orbit. , Let be the state vectors of spacecraft i along the x, y, and z axes in the virtual moving point LVLH coordinate system. Let be the position vectors of spacecraft i along the x, y, and z axes in the virtual moving point LVLH coordinate system. ; Let be the velocity vectors of spacecraft i along the x, y, and z axes in the virtual moving point LVLH coordinate system. ; Let x be the perturbation acceleration vector of spacecraft i along the x, y, and z axes in the virtual moving point LVLH coordinate system;

[0030] R represents the coefficients of the second-order harmonic term. e The average radius of the Earth The distance from the center of the virtual moving point. Indicates the inclination angle of the reference orbit. The latitude argument of the virtual moving point.

[0031] In step two, the objective function is:

[0032] ;

[0033] In the formula, Let u be the objective function. P u E These are the velocity increment sequences of the service spacecraft and the target spacecraft at different times, respectively. To serve the spacecraft's position vector at the terminal moment, denoted as the position vector of the target spacecraft at the terminal moment.

[0034] In step two, the constraints are:

[0035] ;

[0036] ;

[0037] ;

[0038]

[0039] In the formula, Defined as the orbital position constraint for spacecraft i; where the colon ":" indicates a definition. Let be the position vector of spacecraft i after its k-th orbital maneuver. Let be the position vector of spacecraft i before the k-th orbit change, where k is the orbit change sequence number;

[0040] Defined as the velocity increment constraint for spacecraft i. Let be the velocity vector of spacecraft i after its k-th orbital maneuver. Let be the velocity vector of spacecraft i before its k-th orbital maneuver. For the maximum speed increment, This represents the maximum number of orbital changes.

[0041] Defined as the time constraint for spacecraft i. For the mission window time of spacecraft i, For terminal time;

[0042] K is defined as the terminal illumination constraint. This is the critical value for the illumination angle.

[0043] In step two, the Hamiltonian function constructed is:

[0044] ;

[0045] In the formula, H is the Hamiltonian function. The Hamiltonian function constructed for spacecraft i, Let i be the costate variable of spacecraft i. Let t be the state transition equation for spacecraft i, where t is the time variable.

[0046] In step two, the costate equation is:

[0047] ;

[0048] The interior-point function is: ;

[0049] ;

[0050] The transverse condition includes the transverse condition at the interior point and the transverse condition at the terminal point; wherein, the transverse condition at the interior point is:

[0051] ;

[0052] ;

[0053] The transverse condition at the terminal is:

[0054] ;

[0055] in, ;

[0056] Calculations show that, ;

[0057] In the formula, Let i be the rate of change of the costate variable at position i. Let i be the rate of change of the velocity costate variable of spacecraft i;

[0058] Let i be an interior point function of spacecraft i. For the Lagrange multipliers constrained by the orbital change position, For Lagrange multipliers constrained by velocity increment, For time-constrained Lagrange multipliers; Lagrange multipliers constrained by terminal illumination;

[0059] Let be the state vector of spacecraft i after its k-th orbital maneuver. Let i be the state vector of spacecraft i before its k-th orbital maneuver. For the velocity increment of spacecraft i; Let i be the costate variable before the k-th orbital maneuver of spacecraft i. Let be the costate variable of spacecraft i after its k-th orbital maneuver.

[0060] This is the terminal function, representing the state values ​​of the serving spacecraft and the target spacecraft at the termination time. To serve the spacecraft's time constraints, Lagrange multipliers To serve the time constraints of spacecraft, For the time-constrained Lagrange multipliers of the target spacecraft, For the time constraints of the target spacecraft, Let i be the costate variable of spacecraft i at the terminal moment. Let be the state vector of spacecraft i at the terminal moment.

[0061] In step three, the method for finding the Nash equilibrium solution of the Nash equilibrium model using the adaptive differential evolution algorithm is as follows:

[0062] S1, pulse timing and Encoded as a velocity increment sequence; the initial decision variables for the velocity increment sequence for: ,in, , ;

[0063] S2, initial values ​​of randomly constructed costate variables As parameters to be optimized, and the initial decision variables Jointly construct optimization variables ;in, These are the initial values ​​of the position covariance and velocity covariance of spacecraft i, respectively.

[0064] S3, the deviation is calculated using a zero-sum game model based on the Nash equilibrium model. The calculation method is as follows:

[0065] ;

[0066] In the formula, , , , , These are, respectively, the costate continuity deviation at the pulse moment, the costate and target gradient deviation at the terminal moment, the velocity increment constraint deviation, the time constraint deviation, and the terminal illumination constraint deviation;

[0067] S4, construct the total deviation function based on the deviation. The total deviation function is:

[0068] ;

[0069] In the formula, Let M1, M2, M3, M4, and M5 be the optimization variables for minimizing the total deviation function; M1, M2, M3, M4, and M5 are respectively... , , , , The corresponding deviation affects the weight;

[0070] S5, when minimizing the objective variable If the error precision is less than or equal to the current error, the iteration stops, and the current optimized variable is output as the Nash equilibrium solution; otherwise, the adaptive differential evolution algorithm (ADE) is used to optimize the variables in the total deviation function. Perform iterative updates until the convergence condition is met, and output the current optimization variable as the Nash equilibrium solution; if convergence is not achieved even after reaching the upper limit of the number of iterations, output the current optimization variable as the approximate Nash equilibrium solution; the iterative update method for the optimization variable is as follows:

[0071] ;

[0072] In the formula, For the optimized variables after iterative updates, For the current optimization variables, for Iteratively update the gradient. To update the step size iteratively.

[0073] The beneficial effects of this invention are:

[0074] 1. The constructed dynamic model breaks through the limitation of traditional models that ignore illumination conditions. It transforms the relative positional relationship between the sun, the target, and the servicing spacecraft into terminal illumination constraints. Coupled with the discrete dynamic characteristics of pulse thrust (such as the Keplerian free motion phase and the upper limit of pulse increment) and the J2 term perturbation of the non-spherical Earth, a unified model framework including orbital dynamics, imaging feasibility, and control constraints is constructed. The dynamic model integrates the illumination model and the equation of state of relative motion, making the on-orbit servicing mission planning more in line with the actual scenario. It provides clear boundary conditions and prerequisite constraints for subsequent model construction, ensuring that the model not only conforms to the physical characteristics of pulse thrust but also meets the engineering requirements of on-orbit servicing, thus improving the model's relevance and practicality.

[0075] 2. Objective functions are constructed for the opposing objectives of the servicing spacecraft (minimizing terminal distance) and the target spacecraft (maximizing terminal distance), and constraint conditions are formed based on orbital position constraints, velocity increment constraints, time constraints, and terminal illumination constraints. Hamiltonian functions are constructed for both the servicing and target spacecraft. Co-state equations are constructed from the partial derivatives of the Hamiltonian functions with respect to the state vector, and transverse conditions are obtained based on the objective function, constraint conditions, and interior point functions. The Hamiltonian function, co-state equations, and transverse conditions together form a zero-sum game model, transforming the strategy interaction into a multi-point boundary value problem with impulse interior point constraints, overcoming the deficiency of traditional unilateral optimization in balancing the conflicting interests of both parties. At the same time, by deriving the co-state equations and transverse conditions, a unified solution framework for discrete impulse control and dynamics is established, ensuring the optimality of the strategy.

[0076] 3. The Nash equilibrium model, composed of the dynamic model and the zero-sum game model under the target scenario, improves the engineering applicability of mission planning: By integrating the illumination constraint and the pulse thrust model, the on-orbit service mission planning is more in line with the actual scenario. When the illumination angle meets the constraint at the terminal moment (e.g., ≤90°), the clarity of visible light imaging can be ensured; and accurate decision-making under multiple constraints is achieved: The multi-physics coupling model can simultaneously meet the orbital dynamic constraints, pulse increment limits and illumination conditions, and has good adaptability.

[0077] 4. An adaptive differential evolution algorithm (ADE) is employed to find the Nash equilibrium solution of the Nash equilibrium model. By encoding decision variables and dynamically adjusting mutation factors and crossover probabilities, the algorithm efficiently solves multi-point boundary value problems. This addresses the problem that existing technologies suffer from a lack of theoretical basis for optimal strategy solutions and difficulty in guaranteeing solution validity due to the conflict between the discrete control characteristics of pulse thrust (such as abrupt velocity increment changes and Kepler constraints in the free motion phase) and the continuous time frame of traditional differential games. Furthermore, it overcomes the shortcomings of existing technologies in solving pulse strategies under multiple constraints, where traditional intelligent optimization algorithms suffer from slow convergence, premature convergence, and poor adaptability to high-dimensional constraints, resulting in low solution efficiency and insufficient solution accuracy. This approach achieves a globally efficient search for pulse strategies, significantly improving solution speed and accuracy while ensuring compliance with illumination constraints and pulse increment limits. Attached Figure Description

[0078] The present invention will now be described in further detail with reference to the accompanying drawings and embodiments.

[0079] Figure 1 This is a schematic diagram illustrating the changes in relative distance and speed provided in an embodiment of the present invention.

[0080] Figure 2 This is a schematic diagram of the change in illumination angle provided in an embodiment of the present invention. Detailed Implementation

[0081] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of the present invention.

[0082] This invention provides a Nash equilibrium modeling method for on-orbit servicing of pulsed thrust spacecraft under illumination constraints, the method comprising:

[0083] Step 1: Set the target scenario corresponding to the mission requirements; Under the target scenario, construct a solar position model based on the ecliptic inclination and the average longitude of the sun, and construct an illumination model based on the solar position model; Simultaneously, in the virtual dynamic point LVLH coordinate system, the distance between the virtual origin and the spacecraft is much smaller than the distance to the Earth's center, and only considering the J2 term perturbation of the Earth's non-spherical shape, construct the relative motion state equation through the perturbation acceleration vectors of the service spacecraft and the target spacecraft; The illumination model and the relative motion state equation constitute the dynamic model under the target scenario;

[0084] The virtual moving point LVLH coordinate system is constructed with the virtual origin of the circular track as the center, the x-axis as the radius vector, the y-axis along the tangent, and the z-axis perpendicular to the track plane.

[0085] Step 2: Construct the objective function using the position vectors of the service spacecraft and the target spacecraft at the terminal moment; construct the constraint conditions based on the orbit change position constraint, velocity increment constraint, time constraint, and terminal illumination constraint; construct Hamiltonian functions for the service spacecraft and the target spacecraft respectively; construct the costate equation using the partial derivatives of the Hamiltonian function with respect to the state vector; obtain the transverse condition based on the objective function, constraint conditions, and interior point function; and construct a zero-sum game model using the Hamiltonian function, costate equation, and transverse condition.

[0086] In on-orbit servicing scenarios, the servicing spacecraft and the target spacecraft have conflicting objectives. By defining conflicting objective functions for both and combining them with constraints, the strategy solution for both can be transformed into an equilibrium solution problem with multi-point boundary constraints, thus achieving systematic modeling of the process.

[0087] Step 3: The Nash equilibrium model is composed of the dynamic model and the zero-sum game model under the target scenario, and the Nash equilibrium solution of the Nash equilibrium model is found by using the adaptive differential evolution algorithm.

[0088] In step one, the target scenario corresponding to the task requirements is set as follows:

[0089] Both the service spacecraft and the target spacecraft achieve velocity increments through instantaneous pulses. The pulse timing and duration are the same, and each velocity increment satisfies the velocity increment constraint and the interval satisfies the time constraint.

[0090] Under the constraints of terminal illumination, the camera automatically points to the center of the target, satisfying the observation alignment conditions.

[0091] In step one, the solar position model is as follows:

[0092] ;in,

[0093] ;

[0094] ;

[0095] ;

[0096] In the formula, This provides information about the Sun's position in the geocentric inertial frame of reference. For the terminal time, The ecliptic inclination, d represents the average longitude of the sun; d represents the cumulative number of days, calculated from 12:00:00 on January 1, 2000. It is the angle closest to the point.

[0097] In step one, the lighting model is as follows:

[0098] ;in,

[0099] ;

[0100] In the formula, The illumination angle at the terminal moment. This represents the relative position vector between the target spacecraft and the sun at the terminal moment. This represents the relative position vector between the target spacecraft and the service spacecraft at the terminal moment; denoted as the position vector of the target spacecraft at the terminal moment.

[0101] In step one, the equation of state for relative motion is:

[0102] ;in,

[0103] ;

[0104] In the formula, Let be the acceleration vector of spacecraft i along the x, y, and z axes in the virtual moving point LVLH coordinate system, where i = P represents the serving spacecraft, i = E represents the target spacecraft, and T is the transpose; The average angular velocity of the orbit. , Here, is the gravitational parameter of the central celestial body, and 'a' is the semi-major axis of the orbit. , Let be the state vectors of spacecraft i along the x, y, and z axes in the virtual moving point LVLH coordinate system. Let be the position vectors of spacecraft i along the x, y, and z axes in the virtual moving point LVLH coordinate system. ; Let be the velocity vectors of spacecraft i along the x, y, and z axes in the virtual moving point LVLH coordinate system. ; Let x be the perturbation acceleration vector of spacecraft i along the x, y, and z axes in the virtual moving point LVLH coordinate system;

[0105] R represents the coefficients of the second-order harmonic term. e The average radius of the Earth The distance from the center of the virtual moving point. Indicates the inclination angle of the reference orbit. The latitude argument of the virtual moving point.

[0106] In step two, the objective function is:

[0107] ;

[0108] In the formula, Let u be the objective function. P u E These are the velocity increment sequences of the service spacecraft and the target spacecraft at different times, respectively. To serve the spacecraft's position vector at the terminal moment, denoted as the position vector of the target spacecraft at the terminal moment.

[0109] The magnitude of the miss rate in the objective function not only reflects the effectiveness of the task but also the superiority or inferiority of the strategies employed by both sides. Therefore, the miss rate becomes a key indicator for measuring the outcome of the contest between the two sides, directly affecting the success or failure of the task.

[0110] In step two, the constraints are:

[0111] ;

[0112] ;

[0113] ;

[0114]

[0115] In the formula, Defined as the orbit change position constraint for spacecraft i, it solves the problem of requiring the position at the orbit change point to be continuous since the position cannot be instantaneously moved; where the colon ":" indicates the definition; Let be the position vector of spacecraft i after its k-th orbital maneuver. Let be the position vector of spacecraft i before the k-th orbit change, where k is the orbit change sequence number;

[0116] Defined as the velocity increment constraint for spacecraft i, it indicates that the velocity will change abruptly before and after the pulse orbit change. However, due to the limitations of the spacecraft's thrust system, the velocity impulse applied in a single pulse cannot be limited to an infinitely large value. Let be the velocity vector of spacecraft i after its k-th orbital maneuver. Let be the velocity vector of spacecraft i before its k-th orbital maneuver. For the maximum speed increment, This represents the maximum number of orbital changes.

[0117] Defined as the time constraint for spacecraft i, this addresses the issue that the mission cannot continue indefinitely; therefore, the terminal time must be less than the mission window time interval to ensure compliance with the mission's time constraints. For the mission window time of spacecraft i, For terminal time;

[0118] K is defined as the terminal illumination constraint, which solves the problem of requiring the service spacecraft's observation angle of the target spacecraft to be within the observable range at the terminal moment. This is the critical value for the illumination angle.

[0119] In step two, the Hamiltonian function constructed is:

[0120] ;

[0121] In the formula, H is the Hamiltonian function. The Hamiltonian function constructed for spacecraft i, Let i be the costate variable of spacecraft i. Let t be the state transition equation for spacecraft i, where t is the time variable.

[0122] In step two, the costate equation is:

[0123] ;

[0124] The interior-point function is: ;

[0125] ;

[0126] The transverse condition includes the transverse condition at the interior point and the transverse condition at the terminal point; wherein, the transverse condition at the interior point is:

[0127] ;

[0128] ;

[0129] The transverse condition at the terminal is:

[0130] ;

[0131] in, ;

[0132] Calculations show that, ;

[0133] In the formula, Let i be the rate of change of the costate variable at position i. Let i be the rate of change of the velocity costate variable of spacecraft i;

[0134] Let i be an interior point function of spacecraft i. For the Lagrange multipliers constrained by the orbital change position, For Lagrange multipliers constrained by velocity increment, For time-constrained Lagrange multipliers; Lagrange multipliers constrained by terminal illumination;

[0135] Let be the state vector of spacecraft i after its k-th orbital maneuver. Let i be the state vector of spacecraft i before its k-th orbital maneuver. For the velocity increment of spacecraft i; Let i be the costate variable before the k-th orbital maneuver of spacecraft i. Let be the costate variable of spacecraft i after its k-th orbital maneuver.

[0136] This is the terminal function, representing the state values ​​of the serving spacecraft and the target spacecraft at the termination time. To serve the spacecraft's time constraints, Lagrange multipliers To serve the time constraints of spacecraft, For the time-constrained Lagrange multipliers of the target spacecraft, For the time constraints of the target spacecraft, Let i be the costate variable of spacecraft i at the terminal moment. Let be the state vector of spacecraft i at the terminal moment.

[0137] In step three, the method for finding the Nash equilibrium solution of the Nash equilibrium model using the adaptive differential evolution algorithm is as follows:

[0138] S1, pulse timing and Encoded as a velocity increment sequence; the initial decision variables for the velocity increment sequence for: ,in, , The velocity increments of spacecraft i along the x, y, and z axes in the virtual moving point LVLH coordinate system;

[0139] S2, initial values ​​of randomly constructed costate variables As parameters to be optimized, and the initial decision variables Jointly construct optimization variables ;in, These are the initial values ​​of the position covariance and velocity covariance of spacecraft i, respectively.

[0140] S3, the deviation is calculated using a zero-sum game model based on the Nash equilibrium model. The calculation method is as follows:

[0141] ;

[0142] In the formula, , , , , These are, respectively, the costate continuity deviation at the pulse moment, the costate and target gradient deviation at the terminal moment, the velocity increment constraint deviation, the time constraint deviation, and the terminal illumination constraint deviation;

[0143] S4, construct the total deviation function based on the deviation. The total deviation function is:

[0144] ;

[0145] In the formula, Let M1, M2, M3, M4, and M5 be the optimization variables for minimizing the total deviation function; M1, M2, M3, M4, and M5 are respectively... , , , , The corresponding deviation affects the weight;

[0146] S5, when minimizing the objective variable If the error precision is less than or equal to the current error, the iteration stops, and the current optimized variable is output as the Nash equilibrium solution; otherwise, the adaptive differential evolution algorithm (ADE) is used to optimize the variables in the total deviation function. Perform iterative updates until the convergence condition is met, and output the current optimization variable as the Nash equilibrium solution; if convergence is not achieved even after reaching the upper limit of the number of iterations, output the current optimization variable as the approximate Nash equilibrium solution; the iterative update method for the optimization variable is as follows:

[0147] ;

[0148] In the formula, For the optimized variables after iterative updates, For the current optimization variables, for Iteratively update the gradient. To update the step size iteratively.

[0149] To verify the effectiveness of the technical solution of the present invention, a specific example is provided:

[0150] This example sets constraints based on the requirements of actual space missions. Specific parameter settings for the experimental scenario are shown in Table 1:

[0151] Table 1

[0152]

[0153] Indicates the orbital state of the spacecraft, where , representing the phase of the spacecraft in the orbital plane, t f For simulation time, T is the maximum pulse amplitude that the spacecraft can apply. iThe minimum pulse interval is given, where t represents the scene time. The constraint value of the illumination angle at the terminal moment is used. When the actual illumination angle is less than the critical value, the constraint condition is considered to be met, indicating that the target spacecraft is on the illuminated surface, which can ensure that the service spacecraft can establish effective observation conditions through the optical sensor.

[0154] The servicing spacecraft conducted an on-orbit servicing mission in geostationary orbit with an initial phase difference of 5°. The mission scenario was set to 12:00 on January 4, 2023, lasting for 20,000 seconds. The pulse amplitudes of the servicing and target spacecraft were 50 m / s and 40 m / s, respectively. With limited maneuvering advantages, both spacecraft employed the optimal maneuvers obtained from the solution, resulting in a reduction in their relative distance of 409 km, ultimately achieving equilibrium.

[0155] Figure 1 shows the changes in the relative distance, relative velocity, and pulse maneuver timing between the service spacecraft and the target spacecraft during the mission.

[0156] Figure 2 shows the evolution of the terminal illumination angle. The initial illumination angle is about 130° (outside the favorable observation range), and it continues to decrease as the mission progresses, eventually meeting the terminal illumination constraint of "illumination angle ≤ 90°". This ensures that the service spacecraft can establish favorable observation conditions for the target spacecraft.

[0157] When the maneuverability of the servicing spacecraft is increased, its maneuver advantage relative to the target spacecraft is significantly enhanced, and the equilibrium state tilts towards the servicing spacecraft.

[0158] In summary, the simulation results demonstrate that the proposed Nash equilibrium modeling and solution method exhibits excellent convergence performance and engineering applicability. Simulation results show that the dynamic model of the target scenario, composed of the integrated illumination model and the relative motion state equations, along with the Nash equilibrium model based on a zero-sum game model using Hamiltonian functions, can effectively address the coupling problem between discrete control characteristics and continuous dynamics caused by pulse thrust.

[0159] The adaptive differential evolution algorithm (ADE) is used to solve the problem. By simultaneously optimizing the pulse sequence and the initial values ​​of the costate variables in the decision space, and constructing the total deviation function as the convergence criterion, the efficient search for the Nash equilibrium solution is achieved.

[0160] In the established low-Earth orbit scenario, the technical solution disclosed in this invention can achieve stable convergence within a preset time, with a relative distance decrease of 409 km and an illumination angle satisfying the observation constraint of ≤90°. This verifies the effectiveness and reliability of the method in simultaneously handling multiple constraints such as orbital dynamics, pulse increment limits, and illumination conditions. The experimental results fully demonstrate that the modeling and solution strategy provided by this invention possesses both theoretical rigor and engineering practicality, providing effective decision support for on-orbit servicing missions of pulse-thrust spacecraft under complex constraints.

[0161] The beneficial effects of the embodiments of the present invention are:

[0162] 1. The constructed dynamic model breaks through the limitation of traditional models that ignore illumination conditions. It transforms the relative positional relationship between the sun, the target, and the servicing spacecraft into terminal illumination constraints. Coupled with the discrete dynamic characteristics of pulse thrust (such as the Keplerian free motion phase and the upper limit of pulse increment) and the J2 term perturbation of the non-spherical Earth, a unified model framework including orbital dynamics, imaging feasibility, and control constraints is constructed. The dynamic model integrates the illumination model and the equation of state of relative motion, making the on-orbit servicing mission planning more in line with the actual scenario. It provides clear boundary conditions and prerequisite constraints for subsequent model construction, ensuring that the model not only conforms to the physical characteristics of pulse thrust but also meets the engineering requirements of on-orbit servicing, thus improving the model's relevance and practicality.

[0163] 2. Objective functions are constructed for the opposing objectives of the servicing spacecraft (minimizing terminal distance) and the target spacecraft (maximizing terminal distance), and constraint conditions are formed based on orbital position constraints, velocity increment constraints, time constraints, and terminal illumination constraints. Hamiltonian functions are constructed for both the servicing and target spacecraft. Co-state equations are constructed from the partial derivatives of the Hamiltonian functions with respect to the state vector, and transverse conditions are obtained based on the objective function, constraint conditions, and interior point functions. The Hamiltonian function, co-state equations, and transverse conditions together form a zero-sum game model, transforming the strategy interaction into a multi-point boundary value problem with impulse interior point constraints, overcoming the deficiency of traditional unilateral optimization in balancing the conflicting interests of both parties. At the same time, by deriving the co-state equations and transverse conditions, a unified solution framework for discrete impulse control and dynamics is established, ensuring the optimality of the strategy.

[0164] 3. The Nash equilibrium model, composed of the dynamic model and the zero-sum game model under the target scenario, improves the engineering applicability of mission planning: By integrating the illumination constraint and the pulse thrust model, the on-orbit service mission planning is more in line with the actual scenario. When the illumination angle meets the constraint at the terminal moment (e.g., ≤90°), the clarity of visible light imaging can be ensured; and accurate decision-making under multiple constraints is achieved: The multi-physics coupling model can simultaneously meet the orbital dynamic constraints, pulse increment limits and illumination conditions, and has good adaptability.

[0165] 4. An adaptive differential evolution algorithm (ADE) is employed to find the Nash equilibrium solution of the Nash equilibrium model. By encoding decision variables and dynamically adjusting mutation factors and crossover probabilities, the algorithm efficiently solves multi-point boundary value problems. This addresses the problem that existing technologies suffer from a lack of theoretical basis for optimal strategy solutions and difficulty in guaranteeing solution validity due to the conflict between the discrete control characteristics of pulse thrust (such as abrupt velocity increment changes and Kepler constraints in the free motion phase) and the continuous time frame of traditional differential games. Furthermore, it overcomes the shortcomings of existing technologies in solving pulse strategies under multiple constraints, where traditional intelligent optimization algorithms suffer from slow convergence, premature convergence, and poor adaptability to high-dimensional constraints, resulting in low solution efficiency and insufficient solution accuracy. This approach achieves a globally efficient search for pulse strategies, significantly improving solution speed and accuracy while ensuring compliance with illumination constraints and pulse increment limits.

[0166] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A Nash equilibrium modeling method for on-orbit servicing of pulsed thrust spacecraft under illumination constraints, characterized in that, The method includes: Step 1: Set the target scenario corresponding to the mission requirements; Under the target scenario, construct a solar position model based on the ecliptic inclination and the average longitude of the sun, and construct an illumination model based on the solar position model; Simultaneously, in the virtual dynamic point LVLH coordinate system, the distance between the virtual origin and the spacecraft is much smaller than the distance to the Earth's center, and only considering the J2 term perturbation of the Earth's non-spherical shape, construct the relative motion state equation through the perturbation acceleration vectors of the service spacecraft and the target spacecraft; The illumination model and the relative motion state equation constitute the dynamic model under the target scenario; Step 2: Construct the objective function using the position vectors of the service spacecraft and the target spacecraft at the terminal moment; construct the constraint conditions based on the orbit change position constraint, velocity increment constraint, time constraint, and terminal illumination constraint; construct Hamiltonian functions for the service spacecraft and the target spacecraft respectively; construct the costate equation using the partial derivatives of the Hamiltonian function with respect to the state vector; obtain the transverse condition based on the objective function, constraint conditions, and interior point function; and construct a zero-sum game model using the Hamiltonian function, costate equation, and transverse condition. Step 3: The Nash equilibrium model is composed of the dynamic model and the zero-sum game model under the target scenario, and the Nash equilibrium solution of the Nash equilibrium model is found by using the adaptive differential evolution algorithm.

2. The method as described in claim 1, characterized in that, In step one, the target scenario corresponding to the task requirements is set as follows: Both the service spacecraft and the target spacecraft achieve velocity increments through instantaneous pulses. The pulse timing and duration are the same, and each velocity increment satisfies the velocity increment constraint and the interval satisfies the time constraint. Under the constraints of terminal illumination, the camera automatically points to the center of the target, satisfying the observation alignment conditions.

3. The method as described in claim 2, characterized in that, In step one, the solar position model is as follows: ;in, ； ； ； In the formula, This provides information about the Sun's position in the geocentric inertial frame of reference. For the terminal time, The ecliptic inclination, d represents the average longitude of the sun; d represents the cumulative number of days, calculated from 12:00:00 on January 1, 2000. It is the angle closest to the point.

4. The method as described in claim 3, characterized in that, In step one, the lighting model is as follows: ;in, ； In the formula, The illumination angle at the terminal moment. This represents the relative position vector between the target spacecraft and the sun at the terminal moment. This represents the relative position vector between the target spacecraft and the service spacecraft at the terminal moment; denoted as the position vector of the target spacecraft at the terminal moment.

5. The method as described in claim 4, characterized in that, In step one, the equation of state for relative motion is: ;in, ； In the formula, Let be the acceleration vector of spacecraft i along the x, y, and z axes in the virtual moving point LVLH coordinate system, where i = P represents the serving spacecraft, i = E represents the target spacecraft, and T is the transpose; The average angular velocity of the orbit. , Here, is the gravitational parameter of the central celestial body, and 'a' is the semi-major axis of the orbit. , Let be the state vectors of spacecraft i along the x, y, and z axes in the virtual moving point LVLH coordinate system. Let be the position vectors of spacecraft i along the x, y, and z axes in the virtual moving point LVLH coordinate system. ; Let be the velocity vectors of spacecraft i along the x, y, and z axes in the virtual moving point LVLH coordinate system. ; Let x be the perturbation acceleration vector of spacecraft i along the x, y, and z axes in the virtual moving point LVLH coordinate system; R represents the coefficients of the second-order harmonic term. e The average radius of the Earth The distance from the center of the virtual moving point. Indicates the inclination angle of the reference orbit. The latitude argument of the virtual moving point.

6. The method as described in claim 5, characterized in that, In step two, the objective function is: ； In the formula, Let u be the objective function. P u E These are the velocity increment sequences of the service spacecraft and the target spacecraft at different times, respectively. To serve the spacecraft's position vector at the terminal moment, denoted as the position vector of the target spacecraft at the terminal moment.

7. The method as described in claim 6, characterized in that, In step two, the constraints are: ； ； ； In the formula, Defined as the orbital position constraint for spacecraft i; where the colon ":" indicates a definition. Let be the position vector of spacecraft i after its k-th orbital maneuver. Let be the position vector of spacecraft i before the k-th orbit change, where k is the orbit change sequence number; Defined as the velocity increment constraint for spacecraft i. Let be the velocity vector of spacecraft i after its k-th orbital maneuver. Let be the velocity vector of spacecraft i before its k-th orbital maneuver. For the maximum speed increment, This represents the maximum number of orbital changes. Defined as the time constraint for spacecraft i. For the mission window time of spacecraft i, For terminal time; K is defined as the terminal illumination constraint. This is the critical value for the illumination angle.

8. The method as described in claim 7, characterized in that, In step two, the Hamiltonian function constructed is: ； In the formula, H is the Hamiltonian function. The Hamiltonian function constructed for spacecraft i, Let i be the costate variable of spacecraft i. Let t be the state transition equation for spacecraft i, where t is the time variable.

9. The method as described in claim 8, characterized in that, In step two, the costate equation is: ； The interior-point function is: ; ； The transverse condition includes the transverse condition at the interior point and the transverse condition at the terminal point; wherein, the transverse condition at the interior point is: ； ； The transverse condition at the terminal is: ； in, ; Calculations show that, ; In the formula, Let i be the rate of change of the costate variable at position i. Let i be the rate of change of the velocity costate variable of spacecraft i; Let i be an interior point function of spacecraft i. For the Lagrange multipliers constrained by the orbital change position, For Lagrange multipliers constrained by velocity increment, For time-constrained Lagrange multipliers; Lagrange multipliers constrained by terminal illumination; Let be the state vector of spacecraft i after its k-th orbital maneuver. Let i be the state vector of spacecraft i before its k-th orbital maneuver. For the velocity increment of spacecraft i; Let i be the costate variable before the k-th orbital maneuver of spacecraft i. Let be the costate variable of spacecraft i after its k-th orbital maneuver. This is the terminal function, representing the state values ​​of the serving spacecraft and the target spacecraft at the termination time. To serve the spacecraft's time constraints, Lagrange multipliers To serve the time constraints of spacecraft, For the time-constrained Lagrange multipliers of the target spacecraft, For the time constraints of the target spacecraft, Let i be the costate variable of spacecraft i at the terminal moment. Let be the state vector of spacecraft i at the terminal moment.

10. The method as described in claim 9, characterized in that, In step three, the method for finding the Nash equilibrium solution of the Nash equilibrium model using the adaptive differential evolution algorithm is as follows: S1, pulse timing and Encoded as a velocity increment sequence; the initial decision variables for the velocity increment sequence for: ,in, , ; S2, initial values ​​of randomly constructed costate variables As parameters to be optimized, and the initial decision variables Jointly construct optimization variables ;in, These are the initial values ​​of the position covariance and velocity covariance of spacecraft i, respectively. S3, the deviation is calculated using a zero-sum game model based on the Nash equilibrium model. The calculation method is as follows: ； In the formula, , , , , These are, respectively, the costate continuity deviation at the pulse moment, the costate and target gradient deviation at the terminal moment, the velocity increment constraint deviation, the time constraint deviation, and the terminal illumination constraint deviation; S4, construct the total deviation function based on the deviation. The total deviation function is: ； In the formula, Let M1, M2, M3, M4, and M5 be the optimization variables for minimizing the total deviation function; M1, M2, M3, M4, and M5 are respectively... , , , , The corresponding deviation affects the weight; S5, when minimizing the objective variable If the error precision is less than or equal to the current error, the iteration stops, and the current optimized variable is output as the Nash equilibrium solution; otherwise, the adaptive differential evolution algorithm (ADE) is used to optimize the variables in the total deviation function. Perform iterative updates until the convergence condition is met, and output the current optimization variable as the Nash equilibrium solution; if convergence is not achieved even after reaching the upper limit of the number of iterations, output the current optimization variable as the approximate Nash equilibrium solution; the iterative update method for the optimization variable is as follows: ； In the formula, For the optimized variables after iterative updates, For the current optimization variables, for Iteratively update the gradient. To update the step size iteratively.

Images