Track optimization method, detector and electronic equipment

Abstract

The invention relates to a trajectory optimization method, a detector and electronic equipment, which are used for improving the trajectory accuracy of the detector in a set time period. The method comprises the following steps: taking the earliest time of a set time period as first time, taking an integral step length as a first integral step length, taking an initial motion state as a first motion state, and iteratively executing: determining a first duty ratio function value and a second duty ratio function value based on the first time, duty ratio reference time and the first integral step length; determining a second integral step length based on the first duty ratio function value, the second duty ratio function value and the threshold value; determining a second motion state of a second time based on the first motion state and a second integration step size, and taking the second motion state as a new first motion state, the second integration step size as a new first integration step size and the second time as a new first time until the second time reaches the latest time of the set time period; a trajectory is generated based on the initial motion state and the second motion state of each iterative process.

Description

Technical Field

[0001] The present application relates to the field of aerospace technology, and in particular to a trajectory optimization method, a detector, and electronic equipment. Background Art

[0002] Due to its high specific impulse, electric propulsion can significantly reduce fuel consumption and increase payload mass compared to chemical propulsion, enabling probes to reach more distant deep space targets and increasing the returns of space exploration. Consequently, it has become increasingly widely used in deep space exploration missions in recent years. NASA's (National Aeronautics and Space Administration) DeepSpace-1 pioneered the use of electric propulsion technology for interplanetary transfer. The Dawn mission successfully visited Vesta and Ceres using electric propulsion technology. The Japan Aerospace Exploration Agency's (JAXA) Hayabusa and Hayabusa2 asteroid sample return missions both employed ion electric propulsion engines. With the success of these missions, an increasing number of deep space exploration missions are considering electric propulsion systems.

[0003] The optimal trajectory optimization problem for electric propulsion fuel is a typical bang-bang control problem. The homology method is typically used to solve the optimal fuel control problem, obtaining initial values for the co-state. Based on these initial values, the time-varying relationship between the probe's motion state and the time-varying relationship between the probe's motion state and the time-varying relationship between the probe's motion state and the time-varying relationship are determined. However, in actual electric propulsion engineering mission design, multiple constraints must be considered, significantly increasing the number of thruster on / off switching cycles in the probe. This non-smooth and discontinuous nature of thruster switching can affect the accuracy of the probe's motion state, resulting in a significant deviation between the probe's ideal trajectory and its actual trajectory, making it difficult to obtain the probe's optimal trajectory. Summary of the Invention

[0004] Embodiments of the present application provide a trajectory optimization method, a detector, and an electronic device for improving the accuracy of a generated trajectory of a detector within a set time period.

[0005] In a first aspect, the present application provides a trajectory optimization method, the method comprising:

[0006] Obtaining a co-state initial value, setting an integration step and a duty cycle reference time, and determining an initial motion state of the detector at the earliest time in a set time period based on the co-state initial value;

[0007] Set the earliest time as the first time, set the integration step as the first integration step, and the initial motion state as the first motion state, and trigger the following iteration:

[0008] determining a first duty cycle function value based on the first time and the duty cycle reference time, and determining a second duty cycle function value based on the first time, the first integration step, and the duty cycle reference time;

[0009] comparing a product of the first duty cycle function value and the second duty cycle function value with a set threshold to obtain a first comparison result, determining whether a thrust switching point exists for the detector within a target time interval based on the first comparison result, and determining a second integration step size based on the presence of the thrust switching point, wherein the target time interval is determined based on the first time and the first integration step size;

[0010] Determining a second motion state of the detector at a second time based on the first motion state and the second integration step, using the second motion state as a new first motion state, the second integration step as a new first integration step, and the second time as a new first time, wherein the second time is determined based on the first time and the second integration step;

[0011] When the second time does not reach the latest time of the set time period, the iteration is triggered again, otherwise the iteration ends;

[0012] Based on the initial motion state and the second motion state in each iteration, a trajectory of the detector within a set time period is generated.

[0013] In an embodiment of the present application, the earliest time of the set time period is taken as the first time, the integration step is set as the first integration step, and the initial motion state of the detector at the earliest time determined based on the initial value of the co-state is taken as the first motion state, triggering the following iteration: determining the first duty cycle function value of the detector at the first time based on the first time and the duty cycle reference time; determining the second duty cycle function value of the detector at the latest time in the target time interval based on the first time, the first integration step and the duty cycle reference time; judging whether the detector has a thrust switching point within the target time interval based on the relationship between the product of the first duty cycle function value and the second duty cycle function value and the set threshold value, Improve the accuracy of judging the existence of a thrust switching point due to duty cycle constraints; determine the second integral step based on the existence of the thrust switching point, and determine the second motion state of the detector at the second time based on the first motion state and the second integral step, that is, by adjusting the integral step so that the single-step integration reaches the thrust switching point, thereby avoiding the thrust non-smoothness caused by the power on / off switching and improving the accuracy of the motion state of the detector; use the second motion state as the new first motion state, the second integral step as the new first integral step, and the second time as the new first time, and trigger iteration again if the second time does not reach the latest time of the set time period, otherwise the iteration ends. This application generates the optimal trajectory of the detector within a set time period based on the initial motion state and the second motion state in each iteration process.

[0014] In a possible embodiment, determining whether the probe has a thrust switching point within a target time interval based on the first comparison result includes:

[0015] If the first comparison result indicates that the product value is less than or equal to the set threshold, determining that the detector is at the thrust switching point within the target time interval;

[0016] If the first comparison result indicates that the product value is greater than the set threshold, it is determined that the thrust switching point does not exist in the detector within the target time interval.

[0017] In an embodiment of the present application, based on the size relationship between the product of the first duty cycle function value and the second duty cycle function value and the set threshold, it is determined whether the detector has a thrust switching point within the target time interval, thereby improving the accuracy of the judgment of the existence of the thrust switching point caused by the duty cycle constraint.

[0018] In a possible embodiment, determining the second integration step size based on the existence of the thrust switching point includes:

[0019] If it is determined that the detector is at the thrust switching point within the target time interval, the first duty cycle function value is compared with the set threshold, and the second duty cycle function value is compared with the set threshold; if the first duty cycle function value is greater than or equal to the set threshold, and the second duty cycle function value is less than or equal to the set threshold, the second integration step size is determined based on the first duty cycle function value and the set thrust duration; if the first duty cycle function value is less than or equal to the set threshold, and the second duty cycle function value is greater than or equal to the set threshold, the second integration step size is determined based on the first duty cycle function value, the set thrust duration, and the set forced coasting duration;

[0020] If it is determined that the thrust switching point does not exist in the detector within the target time interval, the first integration step is used as the second integration step.

[0021] In an embodiment of the present application, when it is determined that the detector has a thrust switching point within the target time interval, in order to make the single-step integration reach the thrust switching point, the second integration step is determined based on the first duty cycle function value; when it is determined that the detector does not have a thrust switching point within the target time interval, there is no need to adjust the integration step, and therefore, the first integration step is used as the second integration step.

[0022] In a possible embodiment, before determining the first duty cycle function value based on the first time and the duty cycle reference time, the method further includes:

[0023] determining a first cosine value of a Sun-Earth-Probe SEP angle / a Sun-Probe-Earth SPE angle at the first time based on the first motion state;

[0024] comparing the first cosine value with a set cosine threshold to obtain a second comparison result, and determining whether a thrust switching point exists for the detector within the target time interval based on the second comparison result;

[0025] Based on the existence of the thrust switching point, it is determined whether to update the first integration step size.

[0026] In this embodiment of the present application, a first cosine value of the SEP / SPE angle at a first time is determined based on a first motion state. Based on the magnitude relationship between this first cosine value and a set cosine threshold, a determination is made as to whether the probe has a thrust switching point within a target time interval. This improves the accuracy of determining the presence of a thrust switching point due to a solar eclipse shutdown constraint. Furthermore, based on the presence of the thrust switching point, the present application determines whether to update the first integration step size to ensure single-step integration up to the thrust switching point, thereby avoiding thrust non-smoothing caused by power-on / power-off switching.

[0027] In a possible embodiment, the setting of the cosine threshold includes a first cosine threshold and a second cosine threshold, the first cosine threshold is greater than the second cosine threshold, and the determining, based on the second comparison result, whether the detector has a thrust switching point within the target time interval includes:

[0028] If the second comparison result indicates that the first cosine value is greater than or equal to the first cosine threshold, predicting a second cosine value of the SEP angle / SPE angle at a third time based on the first motion state, the first cosine value, and the first integration step, and comparing the second cosine value with the first cosine threshold; if the second cosine value is less than the first cosine threshold, determining that the probe has the thrust switching point within the target time interval; and if the second cosine value is greater than or equal to the first cosine threshold, determining that the probe does not have the thrust switching point within the target time interval, wherein the third time is determined based on the first time and the first integration step;

[0029] If the second comparison result indicates that the first cosine value is greater than the second cosine threshold, and the first cosine value is less than the first cosine threshold, determining whether the probe has the thrust switching point within the target time interval based on the first time and the set solar transit shutdown constraint interval;

[0030] If the first comparison result indicates that the first cosine value is less than or equal to the second cosine threshold, it is determined that the probe does not have the thrust switching point within the target time interval.

[0031] In an embodiment of the present application, based on the size relationship between the first cosine value and the first cosine threshold and the second cosine threshold, it is determined whether the detector has a thrust switching point within the target time interval, thereby improving the accuracy of the judgment of the existence of the thrust switching point caused by the solar eclipse shutdown constraint.

[0032] In a possible embodiment, determining whether the probe has the thrust switching point within the target time interval based on the first time and the set solar eclipse shutdown constraint interval includes:

[0033] If the first time is within the set solar eclipse shutdown constraint interval, determining that the probe does not have the thrust switching point within the target time interval;

[0034] If the first time is not within the set solar transit shutdown constraint interval, a third motion state of the probe at a fourth time is determined based on the first motion state, the set solar transit early shutdown duration, and the first time. Based on the third motion state and the first integration step, a third cosine value of the SEP angle / SPE angle at a fifth time is predicted, and the third cosine value is compared with the first cosine threshold. If the third cosine value is greater than or equal to the first cosine threshold, it is determined that the probe has the thrust switching point within the target time interval. If the third cosine value is less than the first cosine threshold, it is determined that the probe does not have the thrust switching point within the target time interval. The fourth time is determined based on the set solar transit early shutdown duration and the first time, and the fifth time is determined based on the fourth time and the first integration step.

[0035] In an embodiment of the present application, when the first cosine value is greater than the second cosine threshold and the first cosine value is less than the first cosine threshold, based on the first time and the set solar eclipse shutdown constraint interval, it is determined whether the detector has a thrust switching point within the target time interval, thereby improving the accuracy of the determination of the existence of the thrust switching point caused by the solar eclipse shutdown constraint.

[0036] In a possible embodiment, the determining whether to update the first integration step size based on the existence of the thrust switching point includes:

[0037] If it is determined that the detector is at the thrust switching point within the target time interval, determining the thrust switching point using a set root method, and updating the first integration step size based on the thrust switching point and the first time;

[0038] If it is determined that the detector does not have the thrust switching point within the target time interval, the first integration step size is not updated.

[0039] In an embodiment of the present application, when it is determined that the detector has a thrust switching point within the target time interval, in order to make the single-step integration reach the thrust switching point, the set root method is used to determine the thrust switching point within the target time interval, and the first integration step is updated based on the thrust switching point and the first time; when it is determined that the detector does not have a thrust switching point within the target time interval, the first integration step is not updated.

[0040] In a possible embodiment, after determining the thrust switching point by using the set root method, the method further includes:

[0041] If the thrust switching point is the time of sunrise, the duty cycle reference time is updated based on the thrust switching point and the set sunrise delay power-on time.

[0042] In an embodiment of the present application, when the thrust switching point is the sunrise time, the duty cycle reference time is updated based on the thrust switching point and the set sunrise delay start time, so that the periodic duty cycle glide starts from the most recent start time, thereby ensuring the optimal thruster on / off profile curve, thereby saving fuel consumption.

[0043] In a possible embodiment, before determining the first duty cycle function value based on the first time and the duty cycle reference time, the method further includes:

[0044] determining a first available power of the detector at the first time based on the first motion state and a set detector system operating power;

[0045] predicting a second available power of the detector at a third time based on the first motion state, the first available power, and the first integration step;

[0046] comparing the first available power with a set power threshold to obtain a third comparison result, comparing the second available power with the set power threshold to obtain a fourth comparison result, and determining whether the thrust switching point exists within the target time interval based on the third and fourth comparison results;

[0047] Based on the existence of the thrust switching point, it is determined whether to update the first integration step size.

[0048] In an embodiment of the present application, based on the detector's first available power, first motion state, and first integration step size at a first time, the detector's second available power at a third time is predicted. Based on the magnitude relationship between the first available power and a set power threshold, and the magnitude relationship between the second power and the set power threshold, a determination is made as to whether the detector has a thrust switching point within a target time interval. This improves the accuracy of determining the presence of a thrust switching point due to power-grading constraints. Furthermore, based on the presence of a thrust switching point, the present application determines whether to update the first integration step size to ensure single-step integration to the thrust switching point, thereby avoiding thrust non-smoothing caused by power-on and power-off switching.

[0049] In a possible embodiment, the set power threshold includes a first power threshold and a second power threshold, the first power threshold is greater than the second power threshold, and the determining, based on the third comparison result and the fourth comparison result, whether the detector has the thrust switching point within the target time interval includes:

[0050] If the third comparison result indicates that the first available power is greater than or equal to the first power threshold, and the fourth comparison result indicates that the second available power is greater than or equal to the first power threshold, determining that the probe does not have the thrust switching point within the target time interval;

[0051] If the third comparison result indicates that the first available power is greater than or equal to the first power threshold, and the fourth comparison result indicates that the second available power is less than the first power threshold, determining that the detector is at the thrust switching point within the target time interval;

[0052] If the third comparison result indicates that the first available power is less than or equal to the second power threshold, and the fourth comparison result indicates that the second available power is less than or equal to the second power threshold, determining that the probe does not have the thrust switching point within the target time interval;

[0053] If the third comparison result indicates that the first available power is less than or equal to the second power threshold, and the fourth comparison result indicates that the second available power is greater than the second power threshold, determining that the detector is at the thrust switching point within the target time interval;

[0054] If the third comparison result indicates that the first available power is greater than the second power threshold and less than the first power threshold, the second available power is compared with the first available power and the third available power respectively to obtain a fifth comparison result. Based on the fifth comparison result, it is determined whether the probe has a thrust switching point within the target time interval, and the third available power is the available power of the probe at a sixth time, and the sixth time is earlier than the first time.

[0055] In an embodiment of the present application, based on the size relationship between the first available power and the set power threshold, and the size relationship between the second power and the set power threshold, it is determined whether the detector has a thrust switching point within the target time interval, thereby improving the accuracy of the judgment of the existence of the thrust switching point caused by the power grade constraint.

[0056] In a possible embodiment, determining, based on the fifth comparison result, whether the detector has a thrust switching point within the target time interval includes:

[0057] If the fifth comparison result indicates that the second available power is greater than or equal to the first available power, and the second available power is less than or equal to the third available power, determining that the probe does not have the thrust switching point within the target time interval;

[0058] If the fifth comparison result indicates that the second available power is greater than the third available power, determining that the detector is at the thrust switching point within the target time interval;

[0059] If the fifth comparison result indicates that the second available power is less than the first available power, it is determined that the probe is at the thrust switching point within the target time interval.

[0060] In an embodiment of the present application, when the first available power is greater than the second power threshold and less than the first power threshold, based on the size relationship between the second available power and the first available power and the third available power, it is determined whether the detector has a thrust switching point within the target time interval, thereby improving the accuracy of the determination of the existence of the thrust switching point caused by the power grade constraint.

[0061] In a possible embodiment, the determining whether to update the first integration step size based on the existence of the thrust switching point includes:

[0062] If it is determined that the detector is at the thrust switching point within the target time interval, determining the thrust switching point using a set root method, and updating the first integration step size based on the thrust switching point and the first time;

[0063] If it is determined that the detector does not have the thrust switching point within the target time interval, the first integration step size is not updated.

[0064] In an embodiment of the present application, when it is determined that the detector has a thrust switching point within the target time interval, in order to make the single-step integration reach the thrust switching point, the set root method is used to determine the thrust switching point within the target time interval, and the first integration step is updated based on the thrust switching point and the first time; when it is determined that the detector does not have a thrust switching point within the target time interval, the first integration step is not updated.

[0065] In a possible embodiment, after updating the first integration step size based on the thrust switching point and the first time, the method further includes:

[0066] determining a fourth motion state of the probe at the thrust switching point based on the first motion state and the first integration step size, and determining a fourth available power of the probe at the thrust switching point based on the fourth motion state and a set probe system operating power;

[0067] determining a second co-state initial value of the detector at the thrust switching point based on a first co-state initial value of the detector at the first time and the first integration step;

[0068] determining, based on a correspondence between each available power interval and a thruster specific impulse, a first thruster specific impulse corresponding to the available power interval to which the fourth available power belongs, and a second thruster specific impulse corresponding to the available power interval to which the first available power belongs;

[0069] determining a first switching function value based on the fourth motion state, the second co-state initial value, and the first thruster specific impulse; determining a second switching function value based on the first motion state, the first co-state initial value, and the second thruster specific impulse;

[0070] If the second switching function value is greater than a set homotopy parameter and the first switching function value is less than the set homotopy parameter, the duty cycle reference time is updated based on the thrust switching point.

[0071] In an embodiment of the present application, when the second switching function value of the detector at the first time is greater than the set homology parameter, and the first switching function value of the detector at the thrust switching point is less than the set homology parameter, it indicates that a power-on point exists. Therefore, the duty cycle reference time is updated based on the thrust switching point, so that the periodic duty cycle gliding starts from the most recent power-on time, thereby ensuring the optimal thruster power-on and power-off profile curve, thereby saving fuel consumption.

[0072] In a possible embodiment, the initial co-state value is the initial co-state value of the i-th iteration in the process of solving the detector fuel optimization problem using the homotopy algorithm, 0≤i≤N-1, N is a set number of homotopy steps, and the iteration is triggered again when the second time does not reach the latest time of the set time period. Otherwise, after the iteration ends, the method further includes:

[0073] Based on the second motion state in the last iteration and the set motion state of the detector at the latest time in the set time period, the co-state initial value is updated, and the updated co-state initial value is used as the co-state initial value of the (i+1)th iteration.

[0074] In a second aspect, the present application provides a trajectory optimization device, comprising:

[0075] An acquisition module is used to acquire an initial value of the co-state, set an integration step and a duty cycle reference time, and determine the initial motion state of the detector at the earliest time in a set time period based on the initial value of the co-state;

[0076] An iteration module, comprising a first determination module, a comparison module, a second determination module, and a judgment module, wherein the iteration module is configured to set the earliest time as the first time, set the integration step as the first integration step, and the initial motion state as the first motion state, and trigger the following iteration:

[0077] The first determining module is configured to determine a first duty cycle function value based on the first time and the duty cycle reference time, and determine a second duty cycle function value based on the first time, the first integration step, and the duty cycle reference time;

[0078] the comparison module being configured to compare a product of the first duty cycle function value and the second duty cycle function value with a set threshold value to obtain a first comparison result, determine whether a thrust switching point exists for the detector within a target time interval based on the first comparison result, and determine a second integration step size based on the presence of the thrust switching point, wherein the target time interval is determined based on the first time and the first integration step size;

[0079] The second determination module is configured to determine a second motion state of the detector at a second time based on the first motion state and the second integration step, use the second motion state as a new first motion state, the second integration step as a new first integration step, and the second time as a new first time, wherein the second time is determined based on the first time and the second integration step;

[0080] The judgment module is configured to trigger iteration again if the second time does not reach the latest time of the set time period, otherwise the iteration ends;

[0081] A generation module is used to generate a trajectory of the detector within a set time period based on the initial motion state and the second motion state in each iterative process.

[0082] In a third aspect, the present application provides a probe comprising a thruster and a controller, wherein the controller is configured to execute the steps included in the method according to any one of the first aspects.

[0083] In a fourth aspect, the present application provides an electronic device, comprising:

[0084] a memory for storing program instructions;

[0085] The processor is configured to call the program instructions stored in the memory and execute the steps included in any one of the methods of the first aspect according to the obtained program instructions.

[0086] In a fifth aspect, the present application provides a computer-readable storage medium, wherein the computer-readable storage medium stores a computer program, wherein the computer program includes program instructions, and when the program instructions are executed by a computer, the computer executes any one of the methods described in the first aspect.

[0087] In a sixth aspect, the present application provides a computer program product, comprising: a computer program code, which, when executed on a computer, enables the computer to execute any one of the methods described in the first aspect. BRIEF DESCRIPTION OF THE DRAWINGS

[0088] Figure 1 A flow chart of a trajectory optimization method provided in an embodiment of the present application;

[0089] Figure 2 A detailed flow chart of a trajectory optimization method provided in an embodiment of the present application;

[0090] Figure 3 A flow chart of a method for detecting a solar transit constraint thrust switching point provided in an embodiment of the present application;

[0091] Figure 4 A flow chart of a method for detecting power gear switching time provided in an embodiment of the present application;

[0092] Figure 5 A flow chart of a duty cycle reference time updating method provided in an embodiment of the present application;

[0093] Figure 6 Flowchart for implementing the trajectory optimization method provided in the embodiment of the present application;

[0094] Figure 7 A graph showing the change in thrust ratio over time during the homotopy iteration process provided in an embodiment of the present application;

[0095] Figure 8 A schematic diagram of the Earth-Mars transfer flight trajectory and push-slide arc provided in an embodiment of the present application;

[0096] Figure 9 A graph showing thrust variation over time provided in an embodiment of the present application;

[0097] Figure 10 A graph showing the change of thruster specific impulse over time provided in an embodiment of the present application;

[0098] Figure 11 A graph showing the thrust ratio and SEP angle over time provided in an embodiment of the present application;

[0099] Figure 12 A schematic diagram of a solar eclipse shutdown arc provided in an embodiment of the present application;

[0100] Figure 13 A graph showing the switching function changing over time provided in an embodiment of the present application;

[0101] Figure 14A structural diagram of a trajectory optimization device provided in an embodiment of the present application;

[0102] Figure 15 A structural diagram of a detector provided in an embodiment of the present application;

[0103] Figure 16 A structural diagram of an electronic device provided in an embodiment of the present application. DETAILED DESCRIPTION

[0104] In order to make the purpose, technical solutions and advantages of the present application clearer, the technical solutions in the embodiments of the present application will be clearly and completely described below in conjunction with 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 of the embodiments. Based on the embodiments in the present application, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of this application. Unless there is a conflict, the embodiments in the present application and the features in the embodiments can be combined with each other in any way. In addition, although a logical order is shown in the flowchart, in some cases, the steps shown or described can be performed in an order different from that here.

[0105] The terms "first" and "second" in the specification and claims of this application and the above-mentioned drawings are used to distinguish different objects, rather than to describe a specific order. In addition, the term "comprising" and any of its variations are intended to cover non-exclusive protection. For example, a process, method, system, product or device that includes a series of steps or units is not limited to the listed steps or units, but optionally also includes steps or units that are not listed, or optionally also includes other steps or units inherent to these processes, methods, products or devices. "Multiple" in this application can mean at least two, for example, two, three or more, and the embodiments of this application are not limited thereto.

[0106] The following description of exemplary embodiments of the present application is made in conjunction with the accompanying drawings, which include various details of the embodiments of the present application to facilitate understanding, and they should be considered as merely exemplary. Therefore, those of ordinary skill in the art should recognize that various changes and modifications can be made to the embodiments described herein without departing from the scope disclosed in this application. Similarly, for the sake of clarity and conciseness, the description of well-known functions and structures is omitted in the following description. It should be noted that in the embodiments of the present application, certain software, components, models and other existing solutions in the industry may be mentioned, which should be considered as exemplary, and their purpose is only to illustrate the feasibility of the implementation of the technical solution of the present application, but it does not mean that the applicant has or will necessarily use the solution.

[0107] In the technical solution of this application, the acquisition, transmission, storage, and use of data comply with the requirements of relevant national laws and regulations.

[0108] Before introducing the trajectory optimization method provided by the embodiment of the present application, in order to facilitate understanding, the technical background of the embodiment of the present application is first introduced in detail below.

[0109] Space probe: also known as space probe, deep space probe or cosmic probe, is an unmanned spacecraft and the main tool for space exploration that explores the moon and celestial bodies and space beyond the moon. It is divided into lunar probes, planetary and interplanetary probes, small celestial body probes, etc. according to the objects of exploration.

[0110] The optimal trajectory optimization problem for electric propulsion fuel is a typical bang-bang control problem. Homotopy methods are typically used to solve the optimal fuel control problem, obtaining initial values for the co-state. Based on these initial values, the time-varying relationship of the probe's motion state is determined, and the optimal trajectory of the probe is generated based on this relationship. However, in practical electric propulsion mission design, multiple constraints must be considered, significantly increasing the number of thruster on / off switching cycles in the probe. This non-smooth and discontinuous nature of thruster on / off switching can affect the accuracy of the probe's motion state, leading to significant deviations between the ideal and actual trajectory states, making it difficult to obtain the probe's optimal trajectory.

[0111] The above constraints include duty cycle constraints, power bin constraints, and solar eclipse shutdown constraints. The following is a detailed introduction to each of these three constraints:

[0112] Due to duty cycle constraints, the probe needs to plan periodic glide arcs when performing orbit determination, ground communications, or other exploration operations. The electric propulsion duty cycle is usually defined as the proportion of time used for thrust within a deterministic thrust cycle (i.e., a thrust cycle determined only by optimality conditions). Some key operations such as navigation and communication are usually performed in these glide arcs. For example, for a continuous power-on arc, the time starts from the power-on moment, and after 20 days of power-on, it needs to be shut down for 1 day. In this case, the electric propulsion duty cycle is 20 / 21.

[0113] To accommodate power-tier constraints, engine thrust varies in stages with input power. The thrust of the electric propulsion engine is affected by the power of the solar array. As the probe's distance from the sun (the probe-sun distance) increases, the generated power decreases. Depending on the input power, the engine operates at different operating points, each corresponding to a set of thrust and specific impulse values. Specific impulse describes the thruster's jet efficiency and is defined as the impulse generated per unit amount of propellant.

[0114] Due to the solar eclipse shutdown constraint, the spacecraft cannot effectively receive telemetry information during the solar eclipse, which is not conducive to spacecraft status monitoring. Therefore, the spacecraft is usually switched to protection mode, during which the electric propulsion engines are shut down. Solar eclipse shutdown is usually implemented a few days before the solar eclipse and resumes normal mode a few days after the solar eclipse.

[0115] In order to solve the above problems, the present application proposes a trajectory optimization method, a detector and an electronic device to improve the accuracy of the generated trajectory of the detector within a set time period.

[0116] Reference below Figure 1 The flowchart of a trajectory optimization method is shown to illustrate the technical solution provided by the embodiment of the present application:

[0117] Step 101 : Obtain the initial value of the co-state, set the integration step and the duty cycle reference time, and determine the initial motion state of the detector at the earliest time in the set time period based on the initial value of the co-state.

[0118] The integration step size and time period can be set based on actual conditions. The initial co-state value can be obtained by using the homotopy algorithm to solve the detector fuel optimization problem, or it can be the initial co-state value of the i-th iteration of the homotopy algorithm, where 0 ≤ i ≤ N-1, where N is the number of homotopy steps. The initial co-state value of the 0th iteration is determined when the detector energy is optimal. The detector's motion state includes the position vector, velocity vector, and mass.

[0119] In the embodiment of the present application, the duty cycle reference time is usually initialized to t ref = 0, and is updated according to the new start-up time obtained by the switching detection during the dynamic integration process, so that the periodic duty cycle coasting starts from the most recent start-up time, thereby ensuring the optimal thruster on / off profile curve. The duty cycle reference time t in this application is ref The power-on time can be updated based on the power-on time determined by the zero-crossing point of the power-on / off function, the power-on time caused by power switching, and the power-on time after solar eclipse.

[0120] In step 102 , the earliest time is taken as the first time, the integration step is set as the first integration step, the initial motion state is set as the first motion state, and the following iteration is triggered.

[0121] Step 103 : determining a first duty cycle function value based on the first time and the duty cycle reference time, and determining a second duty cycle function value based on the first time, the first integration step, and the duty cycle reference time.

[0122] Optionally, determining the first duty cycle function value based on the first time and the duty cycle reference time includes determining the first duty cycle function value based on the first time, the duty cycle reference time, a set thrust duration, and a set forced coasting duration, wherein the set thrust duration is a thrust duration of a thrust arc within the duty cycle period.

[0123] Optionally, determining the second duty cycle function value based on the first time, the first integration step and the duty cycle reference time includes: determining the second duty cycle function value based on the first time, the first integration step, the duty cycle reference time, the set thrust duration and the set forced glide duration.

[0124] In step 104, the product of the first duty cycle function value and the second duty cycle function value is compared with a set threshold value to obtain a first comparison result. Based on the first comparison result, it is determined whether the detector has a thrust switching point within the target time interval, and the second integration step is determined based on the existence of the thrust switching point.

[0125] The target time interval is determined based on the first time and the first integration step. For example, the first time is t k , the first integration step is h, and the target time interval is [t k ,t k The threshold value may be set according to actual conditions. For example, the threshold value may be set to 0.

[0126] Optionally, judging whether the detector has a thrust switching point within the target time interval based on the first comparison result includes: if the first comparison result represents a product value less than or equal to a set threshold, determining that the detector has a thrust switching point within the target time interval; if the first comparison result represents a product value greater than the set threshold, determining that the detector does not have a thrust switching point within the target time interval.

[0127] In the embodiment of the present application, the second integration step size is determined based on the existence of the thrust switching point in the above step 104, including the following three cases:

[0128] In the first case, if it is determined that the detector has a thrust switching point within the target time interval, the first duty cycle function value is greater than or equal to the set threshold, and the second duty cycle function value is less than or equal to the set threshold, then the second integration step is determined based on the first duty cycle function value and the set thrust duration.

[0129] In the second case, if it is determined that the detector has a thrust switching point within the target time interval, the first duty cycle function value is less than or equal to the set threshold, and the second duty cycle function value is greater than or equal to the set threshold, then the second integration step size is determined based on the first duty cycle function value, the set thrust duration, and the set forced glide duration.

[0130] In the third case, if it is determined that there is no thrust switching point for the detector within the target time interval, the first integration step will be used as the second integration step.

[0131] Step 105 , based on the first motion state and the second integration step, determine the second motion state of the detector at the second time, use the second motion state as the new first motion state, the second integration step as the new first integration step, and the second time as the new first time.

[0132] The second time is determined based on the first time and the second integration step. For example, the first time is t k , the second integration step is h′, the second time is t k +h′.

[0133] Optionally, based on the first motion state and the second integration step, determining the second motion state of the detector at the second time includes: based on the first motion state and the second integration step, using a 56-order variable step Runge-Kutta numerical integration method to determine the second motion state of the detector at the second time.

[0134] In the embodiment of the present application, the following formula can be used to determine the detector at t k +h′ time state x k+1 :

[0135] x k+1 =RK 56 (@f(x,t,u),x k ,t k ,t k +h′); Formula (1)

[0136] Among them, RK 56 is the 56th order variable step size Runge-Kutta numerical integration method, @f is the dynamic integral equation, x k The detector at t k The motion state at the moment, h′ is the second integration step, and u is the thrust ratio.

[0137] In this embodiment of the present application, if the second time determined based on the first time and the second integration step is later than the latest time of the set time period, the second integration step is updated using the latest time of the set time period and the first time to ensure that the integration reaches the latest time of the set time period, thereby obtaining the second motion state of the detector at the latest time of the set time period. For example, t f is the latest time of the set time period, t k is the first time, and the updated second integral step length h′=t f -t k .

[0138] Step 106: When the second time does not reach the latest time of the set time period, the iteration is triggered again; otherwise, the iteration ends.

[0139] In the embodiment of the present application, the initial value of the co-state is the initial value of the co-state at the i-th iteration in the process of solving the detector fuel optimization problem using the homotopy algorithm, 0≤i≤N-1, N is the set number of homotopy steps, and the iteration is triggered again if the second time does not reach the latest time of the set time period. Otherwise, after the iteration ends, the method further includes: updating the initial value of the co-state based on the second motion state in the last iteration and the set motion state of the detector at the latest time of the set time period, and using the updated initial value of the co-state as the initial value of the co-state for the i+1-th iteration. The number of homotopy steps can be set according to actual conditions.

[0140] Step 107 : generating a trajectory of the detector within a set time period based on the initial motion state and the second motion state in each iterative process.

[0141] In the embodiment of the present application, the above steps 103-104 are specific steps of the duty cycle constraint thrust switching point (switch switching point) detection method.

[0142] In an embodiment of the present application, the earliest time of the set time period is taken as the first time, the integration step is set as the first integration step, and the initial motion state of the detector at the earliest time determined based on the initial value of the co-state is taken as the first motion state, triggering the following iteration: determining the first duty cycle function value of the detector at the first time based on the first time and the duty cycle reference time; determining the second duty cycle function value of the detector at the latest time in the target time interval based on the first time, the first integration step and the duty cycle reference time; judging whether the detector has a thrust switching point within the target time interval based on the relationship between the product of the first duty cycle function value and the second duty cycle function value and the set threshold value, Improve the accuracy of judging the existence of a thrust switching point due to duty cycle constraints; determine the second integral step based on the existence of the thrust switching point, and determine the second motion state of the detector at the second time based on the first motion state and the second integral step, that is, by adjusting the integral step so that the single-step integration reaches the thrust switching point, thereby avoiding the thrust non-smoothness caused by the power on / off switching and improving the accuracy of the motion state of the detector; use the second motion state as the new first motion state, the second integral step as the new first integral step, and the second time as the new first time, and trigger iteration again if the second time does not reach the latest time of the set time period, otherwise the iteration ends. This application generates the optimal trajectory of the detector within a set time period based on the initial motion state and the second motion state in each iteration process.

[0143] In the embodiment of this application, Figure 2A detailed flow chart of a trajectory optimization method provided in an embodiment of the present application is shown in FIG. Figure 2 As shown, it at least includes the following steps 201-214:

[0144] Step 201 : Obtain the initial value of the co-state, set the integration step and the duty cycle reference time, and determine the initial motion state of the detector at the earliest time in the set time period based on the initial value of the co-state.

[0145] In step 202 , the earliest time is taken as the first time, the integration step is set as the first integration step, and the initial motion state is set as the first motion state.

[0146] Step 203 : determining a first duty cycle function value based on the first time and the duty cycle reference time, and determining a second duty cycle function value based on the first time, the first integration step, and the duty cycle reference time.

[0147] In the embodiment of the present application, the thruster duty cycle constraint is modeled as a periodic forced glide constraint, and the duty cycle switching function S is defined as DC =t on -mod(tt ref ,t on +t off ). When the spacecraft is forced to glide, S DC Is a negative value. Among them, t is the current time, t ref is the duty cycle reference time (duty cycle constraint reference time), t on and t off are the thrust duration (thrust duration) and the forced coasting duration of the thrust arc within the duty cycle, respectively. mod(.) is the remainder function.

[0148] The first time is t k , the first integration step is h as an example, the present application can use the following formula to determine the detector at t k The duty cycle function value S at the moment DC (t k ):

[0149] S DC (t k )=t on -mod(t k -t ref ,t on +t off ); formula (2)

[0150] Among them, t k For the first time, t ref is the duty cycle reference time, t on is the thrust duration, t offTo enforce the glide duration, mod(.) represents the remainder function.

[0151] This application can use the following formula to determine the detector at t k +h time duty cycle function value S DC (t k +h):

[0152] S DC (t k + h)=t on -mod(t k +ht ref ,t on +t off ); formula (3)

[0153] Among them, t k is the first time, h is the first integration step, t ref is the duty cycle reference time, t on is the thrust duration, t off To enforce the glide duration, mod(.) represents the remainder function.

[0154] Step 204 : Determine the product of the first duty cycle function value and the second duty cycle function value.

[0155] In the embodiment of the present application, the product value S can be determined using the following formula:

[0156] S=S DC (t k )*S DC (t k +h); Formula (4)

[0157] Among them, S DC (t k ) is the detector at t k The duty cycle function value at the moment, S DC (t k +h) is the detector at t k The duty cycle function value at time +h, where h is the first integration step.

[0158] Step 205 , determining whether the product value is greater than a set threshold, if so, executing step 206 , otherwise executing step 208 .

[0159] Step 206: Determine that the probe does not have a thrust switching point within the target time interval.

[0160] For example, if the threshold is set to 0, the product value S = S DC (t k )*S DC (t k+h)>0, it means that the detector is at t k time and t k +h is in the same on / off state, so the detector is in the target time interval [t k ,t k +h], there is no thrust switching point.

[0161] Step 207: Use the first integration step as the second integration step.

[0162] Step 208 : Determine whether the detector has a thrust switching point within the target time interval.

[0163] For example, if the threshold is set to 0, the product value S = S DC (t k )*S DC (t k +h)≤0, it means that the detector is at t k time and t k +h is in different on / off states at the moment, that is, in the target time interval [t k ,t k +h], there is a switch state change, so the detector is in the target time interval [t k ,t k +h] within the thrust switching point.

[0164] Step 209 : If the first duty cycle function value is greater than or equal to the set threshold and the second duty cycle function value is less than or equal to the set threshold, determine a second integration step size based on the first duty cycle function value and the set thrust duration.

[0165] In the embodiment of the present application, the second integration step length h′ can be determined using the following formula:

[0166] h ′ =t on -S DC (t k ); formula (5)

[0167] Among them, t on To set the thrust duration, S DC (t k ) is the detector at t k The duty cycle function value at time .

[0168] Taking the threshold value as 0 as an example, if the first duty cycle function value S DC (t k )≥0, and the second duty cycle function value S DC (t k +h)≤0, it means that the detector is at t k The detector is always in the power-on state.k +h is in the off state, so it is necessary to use the first duty cycle function value S DC (t k ) and set thrust duration t on , use the above formula (5) to determine the second integration step h′.

[0169] In step 210 , if the first duty cycle function value is less than or equal to the set threshold and the second duty cycle function value is greater than or equal to the set threshold, a second integration step size is determined based on the first duty cycle function value, the set thrust duration, and the set forced coasting duration.

[0170] In the embodiment of the present application, the second integration step length h′ can be determined using the following formula:

[0171] h ′ =t on +t off -S DC (t k ); formula (6)

[0172] Among them, t on To set the thrust duration, S DC (t k ) is the detector at t k Duty cycle function value at time t off To set the forced glide duration.

[0173] Taking the threshold value as 0 as an example, if the first duty cycle function value S DC (t k )≤0, and the second duty cycle function value S DC (t k +h)≥0, it means that the detector is at t k The detector is in the off state at all times. k +h is in the power-on state, therefore, it is necessary to use the first duty cycle function value S DC (t k ), set the thrust duration t on and set the forced coasting duration t off , use the above formula (6) to determine the second integration step h′.

[0174] Step 211 : Determine a second motion state of the detector at a second time based on the first motion state and the second integration step.

[0175] Step 212 , determining whether the second time has not reached the latest time of the set time period, if so, executing step 213 , otherwise executing step 214 .

[0176] Step 213: Use the second motion state as the new first motion state, the second integration step as the new first integration step, and the second time as the new first time.

[0177] Step 214 : Based on the initial motion state and the second motion state in each iteration, a trajectory of the detector within a set time period is generated.

[0178] In an embodiment of the present application, in addition to considering the duty cycle constrained thrust switching point detection during trajectory optimization, it is also necessary to consider the solar eclipse constrained thrust switching point detection. The SEP (Sun-Earth-Probe) angle is defined as the angle between the Earth-Sun vector and the Earth-probe vector, and the SPE (Sun-Probe-Earth) angle is defined as the angle between the probe-Sun vector and the probe-Earth vector. Its value range is [0°, 180°]. In this interval, the cosine function is monotonically decreasing. Based on this property, an analytical prediction formula is derived with the cosine value of the SEP angle / SPE angle as the research object, and the analytical prediction method and the Brent numerical root method are combined to accurately determine the thrust switching point. Therefore, the solar eclipse constrained thrust switching point detection is performed before step 203. Figure 3 This is a flow chart of a method for detecting a solar transit constraint thrust switching point provided in an embodiment of the present application, such as Figure 3 As shown, it at least includes the following steps 301-304:

[0179] Step 301: Determine a first cosine value of a SEP angle / SPE angle at a first time based on a first motion state.

[0180] Specifically, the position vector r of the Earth in the heliocentric inertial system at the first time is determined using the JPL (Jet Propulsion Laboratory) ephemeris. SE , according to the detector heliocentric position vector r SP , calculate the position vector r of the probe relative to the earth EP =r SP -r SE , the position vector r of the sun relative to the earth ES =-r SE Based on the position vector r EP and position vector r ES Determine the first cosine value of the SEP angle at the first time. SP , calculate the position vector r of the sun relative to the detector PS =-r SP , according to the position vector r of the earth in the heliocentric inertial system at the first time SE and the detector's heliocentric position vector rSP , calculate the position vector r of the earth relative to the sun PE =r SE -r SP Based on the position vector r PS and position vector r PE A first cosine value of the SPE angle at a first time is determined.

[0181] Step 302: Compare the first cosine value with a set cosine threshold to obtain a second comparison result.

[0182] The above-mentioned set cosine threshold includes a first cosine threshold and a second cosine threshold, and the first cosine threshold is greater than the second cosine threshold. The first cosine threshold and the second cosine threshold can be set according to actual conditions.

[0183] Step 303: Based on the second comparison result, determine whether the probe has a thrust switching point within the target time interval.

[0184] Step 304: Based on the existence of the thrust switching point, determine whether to update the first integration step size.

[0185] In the embodiment of the present application, in step 303, based on the second comparison result, determining whether the probe has a thrust switching point within the target time interval includes the following three situations:

[0186] In the first case, if the second comparison result indicates that the first cosine value is greater than or equal to the first cosine threshold, then based on the first motion state, the first cosine value and the first integration step, the second cosine value of the SEP angle / SPE angle at the third time is predicted. If the second cosine value is less than the first cosine threshold, it is determined that the thrust switching point of the detector exists within the target time interval; if the second cosine value is greater than or equal to the first cosine threshold, it is determined that the thrust switching point of the detector does not exist within the target time interval.

[0187] The third time is determined based on the first time and the first integration step. For example, the first time is t k , the first integration step is h, the third time is t k +h.

[0188] Optionally, based on the first motion state, the first cosine value and the first integration step, the second cosine value of the SEP angle / SPE angle at the third time is predicted, including: determining the first and second derivatives of the first cosine value based on the first motion state; and predicting the second cosine value of the SEP angle / SPE angle at the third time using an analytical prediction method based on the first cosine value, the first integration step, the first and second derivatives of the first cosine value.

[0189] In the embodiment of the present application, taking the SEP angle as an example, the analytical prediction method is to use Taylor's formula to estimate the function value of the k+1 step when the integration step h is small enough and the value of the k step is known. Therefore, t k The cosine value c of the SEP angle at time +h k,h You can use t k The cosine value c of the SEP angle at the moment k To estimate:

[0190]

[0191] Among them, c k,h t k The cosine value of the SEP angle at time +h, c k t k The cosine of the SEP angle at the moment, t k The first derivative of the cosine of the SEP angle at time, t k The second derivative of the cosine of the SEP angle at time , where h is the first integration step.

[0192] Due to t k The SEP angle at the moment is vector r EP and r ES The angle between them is:

[0193]

[0194] Among them, θ k t k SEP angle at the moment, c k t k The cosine value of the SEP angle at the moment, r EP is the position vector of the probe relative to the Earth, r ES is the position vector of the sun relative to the earth, r EP For r EP The module, r EP =‖r EP ‖, ‖.‖ is the modulo operation, r ES For r ES The module, r ES =‖r ES ‖.

[0195] Among them, the first-order derivative of the cosine function of the SEP angle with respect to time can be obtained analytically:

[0196]

[0197] Among them, c k t kThe cosine of the SEP angle at the moment, c k The first derivative of r EP is the position vector of the probe relative to the Earth, r ES is the position vector of the sun relative to the earth, r EP For r EP The module, r ES For r ES The model, For r EP The first derivative of For r ES The first derivative of .

[0198] Cosine function of SEP angle c k Second derivative with respect to time Can be parsed to obtain:

[0199]

[0200] Among them, c k t k The cosine of the SEP angle at the moment, c k The first derivative of c k The second derivative of r EP is the position vector of the probe relative to the Earth, r ES is the position vector of the sun relative to the earth, r EP For r EP The module, r ES For r ES The model, For r EP The first derivative of For r ES The first derivative of For r EP The second derivative of For r ES The second derivative of .

[0201] In the second case, if the second comparison result indicates that the first cosine value is greater than the second cosine threshold, and the first cosine value is less than the first cosine threshold, then based on the first time and the set solar eclipse shutdown constraint interval, it is determined whether the detector has a thrust switching point within the target time interval.

[0202] The solar transit shutdown constraint interval can be set according to actual conditions, or it can be determined by the last solar transit constraint thrust switching point detection. For example, the initial costate value is the initial costate value of the i-th iteration in the process of solving the probe fuel optimization problem using the homotopy algorithm. When performing the solar transit constraint thrust switching point detection, the solar transit shutdown constraint interval is set to the solar transit shutdown constraint interval determined during the i-1-th iteration, 0≤i≤N-1, and N is the set number of homotopy steps. When the initial costate value is the initial costate value of the 0-th iteration in the process of solving the probe fuel optimization problem using the homotopy algorithm, the solar transit shutdown constraint interval is set to the user setting. When the initial costate value is the costate value obtained after solving the probe fuel optimization problem using the homotopy algorithm, the solar transit shutdown constraint interval is set to the solar transit shutdown constraint interval determined during the N-1-th iteration.

[0203] In the second case, based on the first time and the set solar eclipse shutdown constraint interval, it is determined whether the probe has a thrust switching point within the target time interval, including the following two sub-cases:

[0204] In the first sub-case, if the first time is within the set solar eclipse shutdown constraint interval, it is determined that the probe does not have a thrust switching point within the target time interval.

[0205] In the second sub-case, if the probe is not in the set solar eclipse shutdown constraint interval at the first time, the third motion state of the probe at the fourth time is determined based on the first motion state, the set solar eclipse early shutdown duration and the first time. Based on the third motion state and the first integration step, the third cosine value of the SEP angle / SPE angle at the fifth time is predicted, and the third cosine value is compared with the first cosine threshold. If the third cosine value is greater than or equal to the first cosine threshold, it is determined that the probe has a thrust switching point within the target time interval. If the third cosine value is less than the first cosine threshold, it is determined that the probe does not have a thrust switching point within the target time interval.

[0206] The fourth time is determined based on the set early sun shutdown time and the first time. The set early sun shutdown time can be set according to actual conditions. For example, the first time is t k , set the shutdown time in advance to Δh in , the fourth time is t k +Δh in The fifth time is determined based on the fourth time and the first integration step. For example, the first time is t k , set the shutdown time in advance to Δh in , the first integration step is h, the fifth time is t k +Δh in +h.

[0207] In an embodiment of the present application, based on the first motion state, the set early solar eclipse shutdown time and the first time, the third motion state of the detector at the fourth time is determined, including: based on the first motion state, the set early solar eclipse shutdown time and the first time, using the 56th-order variable step-size Runge-Kutta numerical integration method to determine the third motion state of the detector at the fourth time.

[0208] In the embodiment of the present application, the following formula can be used to determine the detector at t k +Δh in The state of motion at the moment x k+Δh :

[0209] x k+Δh =RK 56 (@f(x,t,u),x k ,t k ,t k +Δh in ); formula (11)

[0210] Among them, RK 56 is the 56th order variable step size Runge-Kutta numerical integration method, @f is the dynamic integral equation, x k The detector at t k The state of motion at the moment, Δh in is the shutdown time before solar eclipse, u is the thrust ratio, and u=0 here.

[0211] In an embodiment of the present application, based on the third motion state and the first integration step, the third cosine value of the SEP angle / SPE angle at the fifth time is predicted, including: based on the third motion state, determining the fourth cosine value of the SEP angle / SPE angle at the fourth time, the first-order derivative and the second-order derivative of the fourth cosine value; based on the fourth cosine value, the first integration step, the first-order derivative and the second-order derivative of the fourth cosine value, using an analytical prediction method to predict the third cosine value of the SEP angle / SPE angle at the fifth time.

[0212] In a third case, if the first comparison result indicates that the first cosine value is less than or equal to the second cosine threshold, it is determined that the probe does not have a thrust switching point within the target time interval.

[0213] In the embodiment of the present application, the determination of whether to update the first integration step size based on the existence of the thrust switching point in step 304 includes:

[0214] If it is determined that there is a thrust switching point within the target time interval for the detector, the set root-finding method is used to determine the thrust switching point, and based on the thrust switching point and the first time, the first integration step size is updated; if it is determined that there is no thrust switching point within the target time interval for the detector, the first integration step size is not updated. Among them, the above set root-finding method can be a root-finding method such as Brent numerical root-finding method, Newton iteration method, secant method, etc.

[0215] Specifically, if the thrust switching point determined by the set root-finding method is the time t when exiting the solar eclipse out , then based on the set delay start-up duration Δh out and the thrust switching point t out , the end time t end = t out + Δh out ; if t end > t k and |t end - t k | < h, where h is the first integration step size, it indicates the end of the solar eclipse shutdown constraint. Therefore, based on the end time t end of the solar eclipse shutdown constraint and the first time t k , the first integration step size h = t end - t k .

[0216] If the thrust switching point determined by the set root-finding method is the start time t beg of the solar eclipse shutdown constraint, if t beg > t k and |t beg - t k | < h, where |.| is the absolute value and h is the first integration step size, it indicates entering the solar eclipse shutdown constraint. Therefore, based on the start time t beg of the solar eclipse shutdown constraint and the first time t k , the first integration step size h = t beg - t k .

[0217] In the embodiment of the present application, after determining the thrust switching point by using the set root-finding method, it further includes: if the thrust switching point is the time when exiting the solar eclipse, the duty cycle reference time is updated based on the thrust switching point and the set delay start-up duration when exiting the solar eclipse. Specifically, the thrust switching point is t out , and the set delay start-up duration when exiting the solar eclipse is Δh out , then the duty cycle reference time t ref = t out + Δh out . Among them, the set delay start-up duration when exiting the solar eclipse can be set according to the actual situation.

[0218] This embodiment of the present application provides a method for detecting the thrust switching point of a solar transit constraint, combining an analytical prediction method with the Brent numerical root method to determine the thrust switching point. The following describes the method in detail using the SEP angle as an example:

[0219] Note k The motion state of the detector at time x k , t k The SEP angle at the moment is θ k , t k The cosine of the SEP angle at the moment is c k =cos(θ k ), t k The predicted SEP angle cosine value at time +h is c k,h , t k The cosine of the true SEP angle at time +h is c k+1 ; Set the latest time (terminal time) of the time period to t f , the integration step is h, and the maximum integration step is H max The SEP angle corresponding to the first cosine value is θ sto , the SEP angle corresponding to the second cosine value is 2θ sto , where θ sto It can be set according to actual conditions, usually θ sto Less than 5 degrees. The first cosine value is recorded as c sto =cos(θ sto ), the second cosine value is recorded as c 2,sto =cos(2θ sto ), when c k ≥c sto The detector enters the solar eclipse. Define the solar eclipse state flag n sto , contains three states: the first state, n sto =0 means it is not within the solar eclipse shutdown constraint; the second state, n sto =1 means it is not within the solar eclipse constraint, but within the solar eclipse shutdown constraint; the third state, n sto =2 means it is within the solar eclipse constraint.

[0220] The time of solar eclipse is recorded as t in The time of the solar eclipse is t out The corresponding solar eclipse shutdown constraint start time is recorded as t beg The end time of the Heri Ling shutdown constraint is t end ; The duration of early shutdown at the beginning of the day is recorded as Δh in The delay time after the sun sets is recorded as Δh out .

[0221] The solar eclipse status detection includes solar eclipse entry detection and solar eclipse exit detection. The detection method includes the following steps:

[0222] Step a, when t k = 0, the initialization time of solar eclipse is t in =INF, time of solar eclipse t out =INF, solar eclipse status flag n sto = 0, the starting time of solar eclipse shutdown constraint is recorded as t beg =INF, the end time of the solar eclipse shutdown constraint is recorded as t end =INF, where INF represents infinity.

[0223] Step b, calculate t using JPL ephemeris k The position vector r of the Earth in the heliocentric inertial system at time SE , according to the detector heliocentric position vector r SP , calculate the position vector r of the probe relative to the earth EP =r SP -r SE , the position vector r of the sun relative to the earth ES =-r SE .

[0224] Step c, according to r EP and r ES Calculate t k The cosine value c of the SEP angle at the moment k , first-order derivative and the second-order derivative

[0225] Step d, the cosine value c k and the first cosine threshold c sto and the second cosine threshold c 2,sto Compare and judge the detector in [t k ,t k +h] interval, and based on the existence of the thrust switching point, determine whether to update the integration step size h, specifically including:

[0226] If c k ≥c sto , then n sto =2, indicating that the detector is k Always within the solar eclipse constraint, based on the cosine value c k , first-order derivative and the second-order derivative Using the analytical prediction method to predict one step forward, we get t k The cosine value c of the predicted SEP angle at time +h k,h , which includes the following two situations:

[0227] ①If c k,h <csto , it indicates that the detector has exited the solar eclipse at time t k +h. Therefore, there is a thrust switching point within the interval [t k , t k +h]. Thus, the Brent numerical root-finding method is used to search for the exact thrust switching point, i.e., the solar eclipse exit time t k , within the interval [t k +h]. Based on the set solar eclipse delay start-up duration Δh out and the thrust switching point t out , the solar eclipse shutdown constraint end time t out is determined as t end = t out +Δh out . The duty cycle constraint reference time t ref is updated as t end . If t end >t k and |t end - t k | < h, where h is the integration step size, then based on the solar eclipse shutdown constraint end time t end and the first time t k , the integration step size h is updated as h = t end - t k .

[0228] ② If c k,h ≥c sto , it indicates that the detector is still within the solar eclipse constraint at time t k +h. It is determined that there is no thrust switching point within the interval [t k , t k +h]. Therefore, the integration step size h is not updated.

[0229] If c k >c 2,sto and c k <c sto , it indicates that the detector is not within the solar eclipse constraint at time t k . At this time, there are two cases:

[0230] ① If t k >t beg and t k <t end , it indicates that the detector is within the solar eclipse shutdown constraint at time t k . Then n sto = 1. It is determined that there is no thrust switching point within the interval [t k , t k +h]. Therefore, the integration step size h is not updated.

[0231] ② Otherwise, it indicates that the detector is at t kIf the time is not within the solar eclipse shutdown constraint, then n sto = 0, and it is necessary to calculate the start time t beg of the solar eclipse shutdown constraint. Let the thrust ratio u = 0, and according to the motion state x k of the detector at time t k and the set early shutdown duration Δh in for entering the solar eclipse, perform numerical integration to obtain the motion state x k of the detector at time t in + Δh k+Δh = RK 56 (@f(x,t,u),x k ,t k ,t k + Δh in ). According to the motion state x k+Δh , calculate the cosine value c k of the SEP angle corresponding to the time t in + Δh k+Δh , the first derivative and the second derivative . According to the integration step size h, the cosine value c k+Δh , the first derivative and the second derivative , use the analytical prediction method to predict one step forward to obtain the predicted cosine value of the SEP angle at time t k + Δh in + h as c k+Δh,h . If c k+Δh,h ≥ c sto , it indicates that the detector has entered the solar eclipse at time t k + Δh in + h, that is, there is a thrust switching point in the interval [t k ,t k + h]. Therefore, use the Brent numerical root-finding method to find the exact thrust switching point in the interval [t k ,t k + h], that is, the start time t beg of the solar eclipse shutdown constraint. If t beg > t k and |t beg - t k | < h, where h is the integration step size, then based on the start time t beg of the solar eclipse shutdown constraint and the first time t k , update the integration step size h = t beg - t k , and, based on the early shutdown duration for entering the solar eclipse recorded as Δh in and the start time t beg of the solar eclipse shutdown constraint, update the time for entering the solar eclipse to t in = tbeg +Δh in , the time of the sun rising and falling t out =INF.

[0232] If c k,h ≤c 2,sto , then n sto =0, indicating that the detector is at t k The time is not within the solar eclipse shutdown constraint, confirm [t k ,t k +h] interval, there is no thrust switching point, so the integration step h is not updated.

[0233] The above method is also applicable to the analytical prediction of SPE angle. It only needs to change r EP and r ES The relevant vector is replaced by r PE and r PS Related vectors.

[0234] In the embodiment of the present application, in the process of trajectory optimization, in addition to considering the duty cycle constraint thrust switching point detection and the sun eclipse constraint thrust switching point detection, it is also necessary to consider the power gear switching time detection. Therefore, the power gear switching time detection is performed before step 203. Figure 4 This is a flow chart of a method for detecting power gear switching time provided in an embodiment of the present application, such as Figure 4 As shown, it at least includes the following steps 401-405:

[0235] Step 401 : determining a first available power of the detector at a first time based on a first motion state and a set detector system operating power.

[0236] Among them, the detector system working power P is set L It refers to the operating power of the spacecraft system, which is the power of each subsystem of the probe (GNC (Guidance, Navigation and Control System), thermal control, payload, etc.)

[0237] In the embodiment of the present application, the thruster available power (abbreviated as available power) P is the solar cell array output power P S and detector system operating power P L The difference is:

[0238]

[0239] Among them, the solar cell array power P S It is a function of the distance r between the instrument and the day, and its expression is P S=ISγ. The distance r between the detector and the sun can be determined according to the motion state of the detector, and the solar radiation flux I = I ⊙ / r 2 , I ⊙ =1367W / m 2 Watts per square meter represents the amount of solar radiation per unit area at a distance of 1 AU (Astronomical Unit). S and γ represent the area and efficiency of the solar cell array, respectively.

[0240] Step 402 : predicting the second available power of the detector at a third time based on the first motion state, the first available power, and the first integration step.

[0241] Optionally, based on the first motion state, the first available power and the first integration step, the second available power of the detector at the third time is predicted, including: determining the first-order derivative and the second-order derivative of the first available power based on the first motion state; based on the first available power, the first integration step, the first-order derivative and the second-order derivative of the first available power, using an analytical prediction method to predict the second available power of the detector at the third time.

[0242] In the embodiment of the present application, the analytical prediction method of the available power P is as follows: if the integration time step h is small enough, and the value of the k-th step is known, the Taylor formula is used to estimate the function value of the k+1-th step. Therefore, t k Available power P at time +h k,h You can use t k Available power P at the moment k To estimate:

[0243]

[0244] Among them, P k The detector at t k Available power at the moment, P k The formula P k =I ⊙ Sγr -2 -P L OK, I ⊙ is the solar radiation energy per unit area at a distance of 1AU, S and γ represent the area and efficiency of the solar cell array respectively, r is the distance between the device and the sun, h is the integration step, The available power P k The first derivative of The available power P k The second derivative of .

[0245] Available power P k The first and second derivatives with respect to time can be expressed as:

[0246]

[0247] in, The available power P k The first derivative of The available power P k The second derivative of I ⊙ is the solar radiation energy per unit area at a distance of 1AU, S and γ represent the area and efficiency of the solar cell array respectively, r is the distance between the device and the sun, is the first-order derivative of the device-day distance r, is the second-order derivative of the device-day distance r.

[0248] The instrument-day distance r and its first-order derivative and second-order derivative can be expressed as:

[0249]

[0250] Among them, x is the x-axis coordinate of the detector in the heliocentric inertial system, y is the y-axis coordinate of the detector in the heliocentric inertial system, z is the z-axis coordinate of the detector in the heliocentric inertial system, and v is the velocity of the detector. is the first-order derivative of x, is the first-order derivative of y, is the first-order derivative of z.

[0251] Step 403: Compare the first available power with the set power threshold to obtain a third comparison result, and compare the second available power with the set power threshold to obtain a fourth comparison result.

[0252] The set power threshold includes a first power threshold and a second power threshold, and the first power threshold is greater than the second power threshold. The first power threshold and the second power threshold can be set according to actual conditions.

[0253] Step 404 , based on the third comparison result and the fourth comparison result, determining whether the probe has a thrust switching point within the target time interval;

[0254] Step 405: Based on the existence of the thrust switching point, determine whether to update the first integration step size.

[0255] In the embodiment of the present application, the determination of whether the probe has a thrust switching point within the target time interval based on the third comparison result and the fourth comparison result in step 404 includes the following five situations:

[0256] Case 1: If the third comparison result indicates that the first available power is greater than or equal to the first power threshold, and the fourth comparison result indicates that the second available power is greater than or equal to the first power threshold, it is determined that the detector does not have a thrust switching point within the target time interval.

[0257] Case 2: If the third comparison result indicates that the first available power is greater than or equal to the first power threshold, and the fourth comparison result indicates that the second available power is less than the first power threshold, it is determined that the probe has a thrust switching point within the target time interval.

[0258] Case 3: If the third comparison result indicates that the first available power is less than or equal to the second power threshold, and the fourth comparison result indicates that the second available power is less than or equal to the second power threshold, it is determined that the detector does not have a thrust switching point within the target time interval.

[0259] Case 4: If the third comparison result indicates that the first available power is less than or equal to the second power threshold, and the fourth comparison result indicates that the second available power is greater than the second power threshold, it is determined that the probe has a thrust switching point within the target time interval.

[0260] Case 5: If the third comparison result indicates that the first available power is greater than the second power threshold and less than the first power threshold, the second available power is compared with the first available power and the third available power respectively to obtain a fifth comparison result. Based on the fifth comparison result, it is determined whether the detector has a thrust switching point within the target time interval.

[0261] The third available power is the available power of the detector at the sixth time, which is earlier than the first time. For example, the first time is t k , the first integration step is h, the sixth time is t k -h.

[0262] In the embodiment of the present application, in the above-mentioned situation 5, based on the fifth comparison result, determining whether the probe has a thrust switching point within the target time interval includes the following three sub-situations:

[0263] In subcase 1, if the fifth comparison result indicates that the second available power is greater than or equal to the first available power, and the second available power is less than or equal to the third available power, it is determined that the probe does not have a thrust switching point within the target time interval.

[0264] In subcase 2, if the fifth comparison result indicates that the second available power is greater than the third available power, it is determined that the probe has a thrust switching point within the target time interval.

[0265] In subcase 3, if the fifth comparison result indicates that the second available power is less than the first available power, it is determined that the probe has a thrust switching point within the target time interval.

[0266] In the embodiment of the present application, the determination of whether to update the first integration step size based on the existence of the thrust switching point in step 405 includes:

[0267] If it is determined that the probe has a thrust switching point within the target time interval, the thrust switching point is determined using the set rooting method, and the first integral step is updated based on the thrust switching point and the first time; if it is determined that the probe does not have a thrust switching point within the target time interval, the first integral step is not updated. The set rooting method can be a rooting method such as Brent numerical rooting method, Newton iteration method, secant method, etc. Specifically, the thrust switching point of the probe within the target time interval is t k,m , the first time is t k , update the first integration step h = t k,m -t k .

[0268] The present application provides a method for detecting the power gear switching time, which combines an analytical prediction method with the Brent numerical root method to determine the thrust switching point, i.e., the power switching time. The specific detection process is as follows:

[0269] Step a, based on the detector at t k The state of motion at the moment x k And set the detector system working power P L , determine the detector at t k Available power P k Available power P k The first derivative of and available power P k The second derivative of

[0270] Step b, based on the integration step h, available power P k , available power P k The first derivative of and available power P k The second derivative of Use the analytical prediction method to predict one step forward and get the detector at t k Available power P at time +h k,h .

[0271] Step c: the available power P k and power analysis prediction value P k,h The first power threshold P0 and the second power threshold P N Compare and judge the detector in [t k ,t k +h] interval, and based on the existence of the thrust switching point, determine whether to update the integration step size h, specifically including:

[0272] ① If P k ≥P0, and P k,h ≥P0, it means that the detector is at t kThe available power at time t is greater than or equal to the maximum power, and the detector at t k +h also has an available power greater than or equal to the maximum power, that is, the power remains unchanged in the interval [t k ,t k +h]. Therefore, there is no thrust switching point in the interval [t k ,t k +h].

[0273] ② If P k ≥P0 and P k,h <P0, it means that the available power of the detector at time t k is greater than or equal to the maximum power, and the available power of the detector at t k +h is less than the maximum power, that is, the power decreases in the interval [t k ,t k +h]. Therefore, there is a thrust switching point in the interval [t k ,t k +h].

[0274] ③ If P k ≤P N , and P k,h ≤P N , it means that the available power of the detector at time t k is less than or equal to the minimum power, and the available power of the detector at t k +h is also less than or equal to the minimum power, that is, the power remains unchanged in the interval [t k ,t k +h]. Therefore, there is no thrust switching point in the interval [t k ,t k +h].

[0275] ④ If P k ≤P N , and P k,h >P N , it means that the available power of the detector at time t k is less than or equal to the minimum power, and the available power of the detector at t k +h is greater than the minimum power, that is, the power increases in the interval [t k ,t k +h]. Therefore, there is a thrust switching point in the interval [t k ,t k +h].

[0276] ⑤ If P N <P k <P0, P k,h ≥P k and P k,h≤P k-1 , it means that the available power of the detector at time t k is the intermediate power, and the power remains unchanged within the interval [t k , t k +h]. Therefore, there is no thrust switching point within the interval [t k , t k +h]. Here, P k-1 is the available power of the detector at time t k -h.

[0277] ⑥ If P N <P k <P0, and P k,h >P k-1 , it means that the power increases within the interval [t k , t k +h]. Therefore, there is a thrust switching point within the interval [t k , t k +h].

[0278] ⑦ If P N <P k <P0, and P k,h <P k , it means that the power decreases within the interval [t k , t k +h]. Therefore, there is a thrust switching point within the interval [t k , t k +h].

[0279] Step d, if it is determined that there is a thrust switching point within the interval [t k , t k +h], the Brent root-finding method is used to solve for the exact thrust switching point t k , t k +h], and the integration step size h = t k,m -t k,m -t k is updated.

[0280] Since the switch function value is related to the specific impulse, the specific impulse changes after the power switch, causing the switch function value to jump and triggering the on-off switch. Therefore, after updating the first integration step size based on the thrust switching point and the first time, the duty cycle reference time needs to be updated. Figure 5 This is the flowchart of a method for updating the duty cycle reference time provided by an embodiment of the present application. As Figure 5 shown, it at least includes the following steps 501-505:

[0281] Step 501: Based on the first motion state and the first integration step, determine the fourth motion state of the detector at the thrust switching point; based on the fourth motion state and the set detector system operating power, determine the fourth available power of the detector at the thrust switching point.

[0282] In an embodiment of the present application, determining the fourth motion state of the probe at the thrust switching point based on the first motion state and the first integration step size includes: determining the fourth motion state of the probe at the thrust switching point using a numerical integration method based on the first motion state and the first integration step size. The numerical integration method may be a 56th-order variable-step Runge-Kutta numerical integration method or another numerical integration method.

[0283] Step 502 : Determine a second co-state initial value of the detector at the thrust switching point based on the first co-state initial value and the first integration step of the detector at the first time.

[0284] In an embodiment of the present application, determining a second co-state initial value of the detector at a thrust switching point based on a first co-state initial value and a first integration step size of the detector at a first time includes: determining the second co-state initial value of the detector at the thrust switching point using a numerical integration method based on the first co-state initial value and the first integration step size of the detector at the first time. The numerical integration method may be a 56-order variable-step Runge-Kutta numerical integration method or another numerical integration method.

[0285] Step 503 : Based on the correspondence between each available power interval and the thruster specific impulse, determine a first thruster specific impulse corresponding to the available power interval to which the fourth available power belongs, and a second thruster specific impulse corresponding to the available power interval to which the first available power belongs.

[0286] The solar electric propulsion system has a limited number of operating points, each of which can be represented by a corresponding thrust value T max 、Thruster specific impulse I sp and the range of available power P. Assuming the number of operating points is N, T in the discrete finite power model max and I sp The value of can be expressed as:

[0287] T max =T1,I sp =I sp,1 ,when P∈[P1,+∞);

[0288] …

[0289] T max =T N-1 ,I sp =I sp,N-1 ,when P∈[P N-2 ,PN-1 );

[0290] T max =T N ,I sp =I sp,N ,when P∈[P N-1 ,P N );

[0291] T max =0,I sp =0,when P∈[0,P min ); Formula (19)

[0292] Among them, [P1,+∞),…,[P N-2 ,P N-1 ),[P N-1 ,P N ) represents the corresponding available power range, T1,…,T N-1 ,T N Indicates the maximum available thrust of the propulsion system corresponding to different available power ranges, I sp,1 ,…,I sp,N-1 ,I sp,N represents the corresponding thruster specific impulse, P represents the available power of the thruster, P min =P N-1 Indicates the minimum available thruster input power. When the thruster input power is less than P min , the detector's engine shuts down.

[0293] For example, the detector is at t k Available power P k Belongs to the available power range [P N-1 ,P N ), according to the corresponding relationship between each available power range and thruster specific impulse in the above formula (19), the corresponding thruster specific impulse I is determined sp =I sp,N The detector is at t k Available power P at time +h k,h Belongs to the available power range [P N-2 ,P N-1 ), according to the corresponding relationship between each available power range and thruster specific impulse in the above formula (19), the corresponding thruster specific impulse I is determined sp,h =I sp,N-1 .

[0294] Step 504 : Determine a first switching function value based on the fourth motion state, the second costate initial value, and the first thruster specific impulse; and determine a second switching function value based on the first motion state, the first costate initial value, and the second thruster specific impulse.

[0295] In the embodiment of the present application, the switch function value can be calculated as follows:

[0296]

[0297] Among them, ρ represents the switching function value, I sp is the thruster specific impulse, g0 is the gravity acceleration at sea level, m is the mass of the probe, and λ v is the co-state variable corresponding to the detector motion state v (velocity vector), λ m is the costate variable corresponding to the mass m of the detector, and λ0 represents the normalized multiplier.

[0298] Step 505 : If the second switching function value is greater than the set homotopy parameter and the first switching function value is less than the set homotopy parameter, then the duty cycle reference time is updated based on the thrust switching point.

[0299] In an embodiment of the present application, if the initial value of the co-state is the co-state value obtained after solving the detector fuel optimization problem using the homology algorithm, the homology parameter is set to 0; if the initial value of the co-state is the initial value of the co-state of the i-th iteration in the process of solving the detector fuel optimization problem using the homology algorithm, the homology parameter is set to the homology parameter corresponding to the i-th iteration, 0≤i≤N-1, N is the set number of homology steps.

[0300] For example, the detector is at t k The state of motion at the moment x k , the detector is at t k The initial value of the co-state at time λ(t k ), the detector is at t k The thruster specific impulse I corresponding to the moment sp , based on the motion state x k , initial value of co-state λ(t k ) and thruster specific impulse I sp , use the above formula (20) to determine the detector at t k The switching function value ρ at the moment k , the detector is at t k The state of motion at the moment x k,h , the detector is at t k The initial value of the co-state at time λ(t k +h), the detector is at t k +h moment, that is, the thrust switching point t k,m , the corresponding thruster specific impulse I sp,h , use the above formula (20) to determine the detector at t k The switching function value ρ at time +h k,h If ρ k >ε and ρ k,h <ε, indicating that the detector is [tk ,t k +h] interval has a thrust switching point, i.e., the start-up point, so the duty cycle reference time t is updated. ref =t k,m .

[0301] In an embodiment of the present application, a classical homology method is used to solve the fuel optimal control problem to obtain the initial value of the co-state when the detector fuel is optimal; based on the initial value of the co-state when the detector fuel is optimal, the trajectory of the detector within a set time period is generated. Figure 6 The implementation flow chart of the trajectory optimization method provided in the embodiment of the present application is as follows: Figure 6 As shown, it at least includes the following steps 601-604:

[0302] Step 601: Construct a two-point boundary value problem shooting equation according to the optimal control principle.

[0303] Assuming that the probe is only affected by the heliocentric gravity and the thrust of the solar electric propulsion system, its dynamic equation can be expressed as:

[0304]

[0305] Where r and v are the position and velocity vectors of the probe in the heliocentric inertial system, r = ‖r‖ is the distance from the probe to the sun (the probe-sun distance), and m is the mass of the probe. is the first-order derivative r of the position vector, is the first-order derivative of the velocity vector v, is the first-order derivative r of mass m, T max is the maximum thrust of the detector, u∈[0,1] is the thrust ratio, α is the thrust direction unit vector, I sp is the thruster specific impulse; μ is the solar gravitational constant, and g0 is the acceleration of gravity at the Earth’s sea level.

[0306] Initial and final state constraints:

[0307] r(t0)=r0,v(t0)=v0,m(t0)=m0,r(t f )=r f ,v(t f )=v f ,m(t f )>0; Formula (22)

[0308] Among them, t0 represents the initial time of the detector, that is, the earliest time of the set time period, t f The terminal time of the detector is the latest time of the set time period. The position, velocity vector and mass of the detector at the departure time are known. The position and velocity vector of the detector at the target time are known, and the mass of the detector at the target time is greater than zero.

[0309] The optimization goal is:

[0310]

[0311] Where J is the consumed fuel, λ0 is the normalization multiplier, λ0>0, and ε is the homotopy parameter that links energy optimization and fuel optimization. When the homotopy parameter ε=1, it is the energy optimal index, and when the homotopy parameter ε=0, it is the fuel optimal index. max is the maximum thrust of the detector, u∈[0,1] is the thrust ratio, t0 is the initial time of the detector, t f is the terminal moment of the detector, I sp is the thruster specific impulse, and g0 is the acceleration of gravity at sea level.

[0312] Hamiltonian function:

[0313]

[0314] Among them, λ r ,λ v and λ m Represent the detector states: position vector r, velocity vector v and co-state variables corresponding to mass m. r and v are the position and velocity vectors of the detector in the heliocentric inertial system, r = ‖r‖ is the distance from the detector to the sun, μ is the solar gravitational constant, α is the thrust direction unit vector, T max is the maximum thrust of the detector, I sp is the thruster specific impulse, g0 is the acceleration of gravity at sea level, u∈[0,1] is the thrust magnitude ratio, m is the mass of the probe, λ0 represents the normalized multiplier, λ0>0, and ε is the homotopy parameter.

[0315] According to optimal control theory, the optimal thrust vector direction α (thrust direction) is:

[0316] α=-λ v / ‖λ v ‖; Formula (25)

[0317] Among them, λ v is the co-state variable corresponding to the velocity vector of the probe, and α is the thrust direction of the thruster.

[0318] Substituting formula (25) into formula (24), the Hamiltonian function can be rewritten as:

[0319]

[0320] Among them, λ r ,λ v and λ mRepresent the detector states: position vector r, velocity vector v and co-state variables corresponding to mass m. r and v are the position and velocity vectors of the detector in the heliocentric inertial system, r = ‖r‖ is the distance from the detector to the sun, μ is the solar gravitational constant, T max is the maximum thrust of the detector, I sp is the thruster specific impulse, g0 is the acceleration of gravity at sea level, u∈[0,1] is the thrust magnitude ratio, m is the mass of the probe, λ0 represents the normalized multiplier, λ0>0, and ε is the homotopy parameter.

[0321] The control force ratio that minimizes the Hamiltonian function is related to the homotopy parameter ε, and its expression is:

[0322]

[0323] Where ρ represents the switching function value, ε represents the homotopy parameter, and u∈[0,1] is the thrust ratio. The thrust of the thruster can be calculated using formula (27).

[0324] In the embodiment of the present application, after the duty cycle constraint thrust switching point is introduced into the trajectory optimization method, the actual control force ratio (thrust magnitude ratio) u is rewritten from formula (27) as u=δ(t)*u.

[0325] Define δ DC is the control ratio under the duty cycle constraint, and its value is determined by the duty cycle switching function, that is,

[0326]

[0327] Among them, t represents the current time, t ref Represents the duty cycle reference time (duty cycle constraint reference time), t on and t off They represent the thrust arc duration and the forced coasting duration within the duty cycle, respectively. mod(.) represents the remainder function.

[0328] The covariate differential equation is:

[0329]

[0330] Among them, λ r ,λ v and λ m Represent the detector states: position vector r, velocity vector v and co-state variables corresponding to mass m, is the covariate variable λ r The first derivative of is the covariate variable λ v The first derivative of is the covariate variable λ mThe first derivative of , μ is the solar gravitational constant, T max is the maximum thrust of the probe, u∈[0,1] is the thrust ratio, r=‖r‖ is the distance from the probe to the sun, and m is the mass of the probe.

[0331] remember (represents a new vector obtained by combining multiple vectors). According to the co-state initial value normalization method, we construct λ(t) = λ(t) / ‖λ(t0)‖, then λ(t0) = 1. It should be noted that λ0 is a constant and does not change with time. is an 8-dimensional vector and T is the transpose.

[0332] by As the independent variable, the shooting equation of the energy optimal control problem is:

[0333] Φ e (λ(t0))=[r(t f )-r f ,v(t f )-v f ,λ m (t f ),‖λ(t0)‖-1] T =0; Formula (30)

[0334] Among them, Φ e (λ(t0)) is a vector function with λ(t0) as the independent variable, t0 represents the initial time of the detector, t f Indicates the terminal time of the detector.

[0335] Fix λ0, and As the independent variable, the target shooting equation of the fuel optimal control problem is:

[0336] Φ f (λ 2:8 (t0))=[r(t f )-r f ,v(t f )-v f ,λ m (t f )] T =0; Formula (31)

[0337] Among them, λ 2:8 (t0) is the 2nd to 8th component of λ(t0), Φ f (λ 2:8 (t0)) is the value of λ 2:8 (t0) is the vector function of the independent variable, t0 represents the initial time of the detector, t f represents the terminal moment of the detector, and T is the transpose.

[0338] Step 602 , using a particle swarm algorithm to estimate the initial value of the co-state and solve the energy optimal control problem, to obtain the initial value of the co-state when the detector energy is optimal.

[0339] Specifically, the particle swarm optimization (PSO) algorithm is used to solve the energy optimal control problem and obtain the initial co-state guess value p 0,pso ; Based on the initial co-state guess value p 0,pso , using Minpack-1 (a library for solving nonlinear equations and nonlinear least squares problems) to solve the energy optimal control problem (corresponding to equation (30), the initial value of the co-state when the detector energy is optimal is obtained 0,energy .

[0340] Minpack-1 is commonly used to solve nonlinear minimization problems. Using Minpack-1 to adjust the initial values of covariates can help converge to the optimal solution more quickly, avoiding local optimal solutions or numerical instability.

[0341] Step 603: Iteratively update the initial value of the co-state determined when the detector energy is optimal according to the set number of homology steps to obtain the initial value of the co-state when the detector fuel is optimal.

[0342] Specifically, a fixed number of homotopy steps N is used, and the initial value of the co-state p when the detector energy is optimal is used. 0,energy Start the homotopy process and use Minpack-1 to solve the homotopy coefficient ε i The optimal control problem (corresponding to formula (31)) is obtained by the initial value of the co-state p 0,i ; with p 0,i Calculate the homotopy coefficient ε for the initial guess value i+1 The optimal control problem of p 0,i+1 ; When ε i When it gradually decreases from 1 to 0, the initial value of the co-state when the detector fuel is optimal is obtained.

[0343] Step 604 : Based on the initial value of the co-state when the detector fuel is optimal, generate the trajectory of the detector within a set time period.

[0344] The optimal fuel control problem is a typical Bang-Bang control problem. The switching function determines the thrust switching time. Further consideration of duty cycle constraints, power step constraints, and solar eclipse shutdown constraints introduces new thrust switching times. Therefore, it is necessary to detect the thrust switching points introduced by these constraints separately. By integrating the switching function detection, duty cycle constraint thrust switching point detection, power step constraints, and solar eclipse shutdown constraints, a multi-constraint thrust switching time detection method is designed during the forward integration of orbital dynamics as follows:

[0345] Step 1, set the initial time \(t = 0\), \(k = 0\), and the terminal time is \(t\). k f .

[0346] Step 2, calculate the thrust and specific impulse of the detector at time \(t\) according to the motion state \(x\), and initialize the integration step size as \(h = H\). If \(t - t < h\), then \(h = t - t\). Then successively execute the following sub-steps: k k max f k f k .

[0347] Step 21, detect the power gear switching time in the interval \([t, t + h]\). If the thrust switching point \(t\) is detected, then update \(h = t - t\); if the thrust switching point is not detected, then \(h\) is not updated. k k switch switch k .

[0348] Step 22, detect the sun outage constraint thrust switching point in the interval \([t, t + h]\). If the detected thrust switching point is the start time \(t\) of the sun outage shutdown constraint, \(t > t\) and \(|t - t| < h\), then update \(h = t - t\); if the detected thrust switching point is the time \(t\) to get out of the sun outage, based on the sun outage delay startup duration \(\Delta h\) and the time \(t\) to get out of the sun outage, determine the end time \(t\) of the sun outage shutdown constraint \(= t+\Delta h\). If \(t > t\) and \(|t - t| < h\), then update \(h = t - t\); if the thrust switching point is not detected, then \(h\) is not updated. k k bef bef k bef k bef k . out out out end out out end k end k end k .

[0349] Step 23, detect the switching function in the interval \([t, t + h]\). If the thrust switching point \(t\) is detected, then update \(h = t - t\). k k switch switch k ​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​; If the thrust switching point is not detected, h is not updated.

[0350] The specific process of the above switch function detection is prior art and will not be described in detail here.

[0351] Step 24, in [t k ,t k +h] interval to detect the duty cycle constraint thrust switching point. If the thrust switching point t switch , then update h = t switch -t k ; If the thrust switching point is not detected, h is not updated.

[0352] In the embodiment of the present application, the execution order of the above steps 21, 22 and 23 can be adjusted according to the actual situation. k ,t k The duty cycle reference time during the duty cycle constraint thrust switching point detection in the [+h] interval may be updated during the power gear switching time detection, the sun eclipse constraint thrust switching point detection, and the switch function detection. Therefore, the duty cycle constraint thrust switching point detection needs to be executed after the power gear switching time detection, the sun eclipse constraint thrust switching point detection, and the switch function detection are executed.

[0353] Step 3: Use the 56th-order variable-step Runge-Kutta method to integrate forward one step: x k+1 =RK 56 (@f(x,t,u),x k ,t k ,t k +h), and get the new motion state x k+1 .

[0354] Step 4: Update parameter t k =t k +h,x k =x k+1 , k=k+1, if t k ≥t f , end the loop and output the motion state x at the terminal moment k ; Otherwise, proceed to step 2.

[0355] The above process completes a dynamic forward integration, that is, taking the motion state x0 at time t = 0 as the initial state, according to the dynamic equation (Formula (21)) and the optimal control force, the integration is completed to t = t f When solving the equation (Formula (30) or (31)) in target shooting, the above integration process will be called multiple times.

[0356] Take the optimal fuel transfer trajectory design from Earth to Mars as an example. Assume that the departure time from Earth is 2012-10-31T00:00:00 Beijing Time, the arrival time on Mars is 2015-10-01T00:00:00 Beijing Time, and the Earth-Mars transfer flight duration is 1065 days. The probe's initial mass m0 = 1000 kg (kilograms), and the solar panel area S = 26 m 2 (square meters), the conversion efficiency of the solar cell array is γ = 15%, and the power of other working systems acting on the detector is P L =500W (watts). The solar radiation flux at 1AU is 1367.0W / m 2 (Watts per square meter). The motion states of the probe at the initial moment and the terminal moment in the heliocentric mean ecliptic inertial coordinate system are shown in Table 1.

[0357] Table 1 Initial and final states of the probe (heliocentric mean ecliptic inertial coordinate system)

[0358] Initial state Terminal State Time / BJT 2012-10-31T00:00:00.0000 2015-10-01T00:00:00.0000 X(m) 117884866947.759 -170529806410.914 Y(m) 90341827889.912 178927997713.575 Z(m) -370789.406 7940962265.422 Vx(m / s) -18595.029 -16624.119 Vy(m / s) 23525.886 -14652.596 Vz(m / s) 1.074 100.459

[0359] The simulation analysis used two NSTAR (North Star) engine models with discrete operating points to optimize the optimal fuel trajectory for Earth-Mars transfer. The corresponding input power, thrust, and specific impulse for each discrete operating point are shown in Table 2. The solar eclipse constraint was set as follows: when the SEP angle is less than 5 degrees, the probe enters the solar eclipse arc, and the probe's electric propulsion engine is shut down from five days before to five days after the solar eclipse. The duty cycle constraint was set as follows: after 20 days of operation, the probe's electric thrusters must be shut down for one day.

[0360] Table 2 Discrete operating points of two NSTAR engine systems

[0361]

[0362]

[0363] In order to reduce the numerical sensitivity, the normalization operation is used in the solution process. The heliocentric gravitational constant GU=1.32712440018×10 20 m 3 / s 2 As the unit of normalized gravitational constant, the distance between the sun and the earth LU=1.49597870691×10 11 m is the normalized length unit; the initial mass of the detector m0 is the normalized mass unit, and the gravitational acceleration at sea level is g0 = 9.80665 m / s 2 The 56-order variable step size Runge-Kutta method is used as the numerical integrator, and the normalized minimum step size is set to H. min =1.0×10 -12 , maximum step length Hmax =4.0×10 -3 The convergence tolerance of MinPack-1 is set to 1.0×10 -13 The particle swarm algorithm is used to estimate the initial value of the co-state when the detector energy is optimal. The initial value of the co-state when the detector energy is optimal is shown in Table 3.

[0364] Table 3 Initial values of optimal energy co-state

[0365]

[0366] The effectiveness of the proposed method is verified by simulation data. This method is applicable to the design of optimal electric propulsion trajectory for deep space probes considering power binning, duty cycle constraints and solar transit constraints. Figure 7 The following is a graph showing the change of the thrust ratio over time during the homotopy iteration process provided in the embodiment of the present application, such as Figure 7 As shown in Figure 4, as the homotopy process proceeds, the thrust ratio u changes from 1 to 0, which realizes the solution of the fuel optimal control problem. The initial value of the co-state when the detector fuel is optimal is obtained as shown in Table 4. Figure 8 The schematic diagram of the Earth-Mars transfer flight trajectory and push-slide arc provided in the embodiment of the present application is as follows: Figure 8 As shown in , due to the consideration of duty cycle constraints, the thrust arc segment presents a piecewise continuous phenomenon. When the fuel is optimal, the thrust and specific impulse change curves with time are as follows: Figure 9 and Figure 10 As shown in Figure 1, due to the consideration of power grade constraints, the thrust and specific impulse change in a step-by-step manner with the daily distance of the vehicle. The curves of thrust ratio and SEP angle changing with time are shown in Figure 1. Figure 11 As shown, the glide arc constraints before and after the solar eclipse are as follows Figure 12 As shown in the figure, the method accurately detects the leading and trailing edges of the solar eclipse shutdown constraint time, where the solar eclipse time range is 939.5659 to 975.5007 days, and the shutdown time range before and after the solar eclipse is 934.5659 to 980.5007 days. Figure 13 The switching function of the embodiment of the present application is provided as a curve diagram over time, such as Figure 13 As shown in Figure 5, there are 8 zero crossing points in the current time period. The first 6 on / off switching are caused by the zero crossing of the switch function value, and the last two are caused by the sudden change of the switch function after the thrust shift. Table 5 shows the zero crossing time of the switch function. ref The fuel consumption during real-time update is 161.892655 kg, and the duty cycle is referenced to time t ref When it is always 0, the fuel consumption is 161.893559 kg. By updating the duty cycle time, the fuel consumption is saved to a certain extent.

[0367] Table 4 Initial values of optimal fuel co-state

[0368]

[0369] Table 5 Switching function zero crossing and transition point time

[0370] Serial number 1 2 3 4 5 6 7 8 Time / day 16.299 185.991 242.688 335.97 436.997 566.793 892.965 931.371

[0371] Based on the same inventive concept, the present application embodiment provides a trajectory optimization device, please refer to Figure 14 , the device comprises:

[0372] An acquisition module 141 is configured to acquire an initial value of the co-state, set an integration step size and a duty cycle reference time, and determine an initial motion state of the detector at the earliest time in a set time period based on the initial value of the co-state;

[0373] The iteration module 142 includes a first determination module 1421, a comparison module 1422, a second determination module 1423, and a judgment module 1424. The iteration module 142 is configured to use the earliest time as the first time, set the integration step as the first integration step, and the initial motion state as the first motion state, and trigger the following iteration:

[0374] The first determining module 1421 is configured to determine a first duty cycle function value based on the first time and the duty cycle reference time, and determine a second duty cycle function value based on the first time, the first integration step, and the duty cycle reference time;

[0375] The comparison module 1422 is configured to compare a product of the first duty cycle function value and the second duty cycle function value with a set threshold value to obtain a first comparison result, determine whether the detector has a thrust switching point within a target time interval based on the first comparison result, and determine a second integration step size based on the presence of the thrust switching point, where the target time interval is determined based on the first time and the first integration step size;

[0376] The second determining module 1423 is configured to determine a second motion state of the detector at a second time based on the first motion state and the second integration step, and use the second motion state as a new first motion state, the second integration step as a new first integration step, and the second time as a new first time, wherein the second time is determined based on the first time and the second integration step;

[0377] The judgment module 1424 is configured to trigger iteration again if the second time does not reach the latest time of the set time period, otherwise the iteration ends;

[0378] The generating module 143 is configured to generate a trajectory of the detector within a set time period based on the initial motion state and the second motion state in each iterative process.

[0379] As an optional implementation, the comparison module 1422 is configured to:

[0380] If the first comparison result indicates that the product value is less than or equal to the set threshold, determining that the detector is at the thrust switching point within the target time interval;

[0381] If the first comparison result indicates that the product value is greater than the set threshold, it is determined that the thrust switching point does not exist in the detector within the target time interval.

[0382] As an optional implementation, the comparison module 1422 is configured to:

[0383] If it is determined that the detector is at the thrust switching point within the target time interval, the first duty cycle function value is compared with the set threshold, and the second duty cycle function value is compared with the set threshold; if the first duty cycle function value is greater than or equal to the set threshold, and the second duty cycle function value is less than or equal to the set threshold, the second integration step size is determined based on the first duty cycle function value and the set thrust duration; if the first duty cycle function value is less than or equal to the set threshold, and the second duty cycle function value is greater than or equal to the set threshold, the second integration step size is determined based on the first duty cycle function value, the set thrust duration, and the set forced coasting duration;

[0384] If it is determined that the thrust switching point does not exist in the detector within the target time interval, the first integration step is used as the second integration step.

[0385] As an optional implementation manner, before determining the first duty cycle function value based on the first time and the duty cycle reference time, the first determining module 1421 is further configured to:

[0386] determining a first cosine value of a Sun-Earth-Probe SEP angle / a Sun-Probe-Earth SPE angle at the first time based on the first motion state;

[0387] comparing the first cosine value with a set cosine threshold to obtain a second comparison result, and determining whether a thrust switching point exists for the detector within the target time interval based on the second comparison result;

[0388] Based on the existence of the thrust switching point, it is determined whether to update the first integration step size.

[0389] As an optional implementation manner, the set cosine threshold includes a first cosine threshold and a second cosine threshold, the first cosine threshold is greater than the second cosine threshold, and the first determining module 921 is configured to:

[0390] If the second comparison result indicates that the first cosine value is greater than or equal to the first cosine threshold, predicting a second cosine value of the SEP angle / SPE angle at a third time based on the first motion state, the first cosine value, and the first integration step, and comparing the second cosine value with the first cosine threshold; if the second cosine value is less than the first cosine threshold, determining that the probe has the thrust switching point within the target time interval; and if the second cosine value is greater than or equal to the first cosine threshold, determining that the probe does not have the thrust switching point within the target time interval, wherein the third time is determined based on the first time and the first integration step;

[0391] If the second comparison result indicates that the first cosine value is greater than the second cosine threshold, and the first cosine value is less than the first cosine threshold, determining whether the probe has the thrust switching point within the target time interval based on the first time and the set solar transit shutdown constraint interval;

[0392] If the first comparison result indicates that the first cosine value is less than or equal to the second cosine threshold, it is determined that the probe does not have the thrust switching point within the target time interval.

[0393] As an optional implementation manner, the first determining module 1421 is configured to:

[0394] If the first time is within the set solar eclipse shutdown constraint interval, determining that the probe does not have the thrust switching point within the target time interval;

[0395] If the first time is not within the set solar transit shutdown constraint interval, a third motion state of the probe at a fourth time is determined based on the first motion state, the set solar transit early shutdown duration, and the first time. Based on the third motion state and the first integration step, a third cosine value of the SEP angle / SPE angle at a fifth time is predicted, and the third cosine value is compared with the first cosine threshold. If the third cosine value is greater than or equal to the first cosine threshold, it is determined that the probe has the thrust switching point within the target time interval. If the third cosine value is less than the first cosine threshold, it is determined that the probe does not have the thrust switching point within the target time interval. The fourth time is determined based on the set solar transit early shutdown duration and the first time, and the fifth time is determined based on the fourth time and the first integration step.

[0396] As an optional implementation manner, the first determining module 1421 is configured to:

[0397] If it is determined that the detector is at the thrust switching point within the target time interval, determining the thrust switching point using a set root method, and updating the first integration step size based on the thrust switching point and the first time;

[0398] If it is determined that the detector does not have the thrust switching point within the target time interval, the first integration step size is not updated.

[0399] As an optional implementation manner, after determining the thrust switching point by using the set root method, the first determining module 1421 is further configured to:

[0400] If the thrust switching point is the time of sunrise, the duty cycle reference time is updated based on the thrust switching point and the set sunrise delay power-on time.

[0401] As an optional implementation manner, before determining the first duty cycle function value based on the first time and the duty cycle reference time, the first determining module 1421 is further configured to:

[0402] determining a first available power of the detector at the first time based on the first motion state and a set detector system operating power;

[0403] predicting a second available power of the detector at a third time based on the first motion state, the first available power, and the first integration step;

[0404] comparing the first available power with a set power threshold to obtain a third comparison result, comparing the second available power with the set power threshold to obtain a fourth comparison result, and determining whether the thrust switching point exists within the target time interval based on the third and fourth comparison results;

[0405] Based on the existence of the thrust switching point, it is determined whether to update the first integration step size.

[0406] As an optional implementation manner, the set power threshold includes a first power threshold and a second power threshold, the first power threshold is greater than the second power threshold, and the first determining module 1421 is configured to:

[0407] If the third comparison result indicates that the first available power is greater than or equal to the first power threshold, and the fourth comparison result indicates that the second available power is greater than or equal to the first power threshold, determining that the probe does not have the thrust switching point within the target time interval;

[0408] If the third comparison result indicates that the first available power is greater than or equal to the first power threshold, and the fourth comparison result indicates that the second available power is less than the first power threshold, determining that the detector is at the thrust switching point within the target time interval;

[0409] If the third comparison result indicates that the first available power is less than or equal to the second power threshold, and the fourth comparison result indicates that the second available power is less than or equal to the second power threshold, determining that the probe does not have the thrust switching point within the target time interval;

[0410] If the third comparison result indicates that the first available power is less than or equal to the second power threshold, and the fourth comparison result indicates that the second available power is greater than the second power threshold, determining that the detector is at the thrust switching point within the target time interval;

[0411] If the third comparison result indicates that the first available power is greater than the second power threshold and less than the first power threshold, the second available power is compared with the first available power and the third available power respectively to obtain a fifth comparison result. Based on the fifth comparison result, it is determined whether the probe has a thrust switching point within the target time interval, and the third available power is the available power of the probe at a sixth time, and the sixth time is earlier than the first time.

[0412] As an optional implementation manner, the first determining module 1421 is configured to:

[0413] If the fifth comparison result indicates that the second available power is greater than or equal to the first available power, and the second available power is less than or equal to the third available power, determining that the probe does not have the thrust switching point within the target time interval;

[0414] If the fifth comparison result indicates that the second available power is greater than the third available power, determining that the detector is at the thrust switching point within the target time interval;

[0415] If the fifth comparison result indicates that the second available power is less than the first available power, it is determined that the probe is at the thrust switching point within the target time interval.

[0416] As an optional implementation manner, the first determining module 1421 is configured to:

[0417] If it is determined that the detector is at the thrust switching point within the target time interval, determining the thrust switching point using a set root method, and updating the first integration step size based on the thrust switching point and the first time;

[0418] If it is determined that the detector does not have the thrust switching point within the target time interval, the first integration step size is not updated.

[0419] As an optional implementation manner, after updating the first integration step size based on the thrust switching point and the first time, the first determining module 1421 is further configured to:

[0420] determining a fourth motion state of the probe at the thrust switching point based on the first motion state and the first integration step size, and determining a fourth available power of the probe at the thrust switching point based on the fourth motion state and a set probe system operating power;

[0421] determining a second co-state initial value of the detector at the thrust switching point based on a first co-state initial value of the detector at the first time and the first integration step;

[0422] determining, based on a correspondence between each available power interval and a thruster specific impulse, a first thruster specific impulse corresponding to the available power interval to which the fourth available power belongs, and a second thruster specific impulse corresponding to the available power interval to which the first available power belongs;

[0423] determining a first switching function value based on the fourth motion state, the second co-state initial value, and the first thruster specific impulse; determining a second switching function value based on the first motion state, the first co-state initial value, and the second thruster specific impulse;

[0424] If the second switching function value is greater than a set homotopy parameter and the first switching function value is less than the set homotopy parameter, the duty cycle reference time is updated based on the thrust switching point.

[0425] As an optional implementation, the initial co-state value is the initial co-state value of the i-th iteration in the process of solving the detector fuel optimization problem using the homotopy algorithm, 0≤i≤N-1, N is a set number of homotopy steps, and the iteration is triggered again when the second time does not reach the latest time of the set time period; otherwise, after the iteration ends, the judgment module 1424 is further used to:

[0426] Based on the second motion state in the last iteration and the set motion state of the detector at the latest time in the set time period, the co-state initial value is updated, and the updated co-state initial value is used as the co-state initial value of the (i+1)th iteration.

[0427] Based on the same inventive concept, an embodiment of the present application provides a detector. Since the detector is the detector in the method in the embodiment of the present invention, and the principle of solving the problem by the detector is similar to that of the method, the implementation of the detector can refer to the implementation of the method, and the repeated parts will not be repeated.

[0428] like Figure 15 As shown, the probe includes a thruster 151 and a controller 152, wherein:

[0429] The controller 152 is configured to perform the following steps:

[0430] Obtaining a co-state initial value, setting an integration step and a duty cycle reference time, and determining an initial motion state of the detector at the earliest time in a set time period based on the co-state initial value;

[0431] Set the earliest time as the first time, set the integration step as the first integration step, and the initial motion state as the first motion state, and trigger the following iteration:

[0432] determining a first duty cycle function value based on the first time and the duty cycle reference time, and determining a second duty cycle function value based on the first time, the first integration step, and the duty cycle reference time;

[0433] comparing a product of the first duty cycle function value and the second duty cycle function value with a set threshold to obtain a first comparison result, determining whether a thrust switching point exists for the detector within a target time interval based on the first comparison result, and determining a second integration step size based on the presence of the thrust switching point, wherein the target time interval is determined based on the first time and the first integration step size;

[0434] Determining a second motion state of the detector at a second time based on the first motion state and the second integration step, using the second motion state as a new first motion state, the second integration step as a new first integration step, and the second time as a new first time, wherein the second time is determined based on the first time and the second integration step;

[0435] When the second time does not reach the latest time of the set time period, the iteration is triggered again, otherwise the iteration ends;

[0436] Based on the initial motion state and the second motion state in each iteration, a trajectory of the detector within a set time period is generated.

[0437] The controller 152 implements the steps of the above-mentioned trajectory optimization method by running the executable instructions, and repeated parts are not repeated here.

[0438] Based on the same inventive concept, the embodiment of the present application provides an electronic device that can realize the functions of the trajectory optimization device discussed above. Please refer to Figure 16 , the device includes a processor 161 and a memory 162, wherein the memory 162 is used to store program instructions;

[0439] The processor 161 calls the program instructions stored in the memory and executes the program instructions to implement the following steps:

[0440] Obtaining a co-state initial value, setting an integration step and a duty cycle reference time, and determining an initial motion state of the detector at the earliest time in a set time period based on the co-state initial value;

[0441] Set the earliest time as the first time, set the integration step as the first integration step, and the initial motion state as the first motion state, and trigger the following iteration:

[0442] determining a first duty cycle function value based on the first time and the duty cycle reference time, and determining a second duty cycle function value based on the first time, the first integration step, and the duty cycle reference time;

[0443] comparing a product of the first duty cycle function value and the second duty cycle function value with a set threshold to obtain a first comparison result, determining whether a thrust switching point exists for the detector within a target time interval based on the first comparison result, and determining a second integration step size based on the presence of the thrust switching point, wherein the target time interval is determined based on the first time and the first integration step size;

[0444] Determining a second motion state of the detector at a second time based on the first motion state and the second integration step, using the second motion state as a new first motion state, the second integration step as a new first integration step, and the second time as a new first time, wherein the second time is determined based on the first time and the second integration step;

[0445] When the second time does not reach the latest time of the set time period, the iteration is triggered again, otherwise the iteration ends;

[0446] Based on the initial motion state and the second motion state in each iteration, a trajectory of the detector within a set time period is generated.

[0447] The processor 161 implements the steps of the above-mentioned trajectory optimization method by running the executable instructions, and repeated parts are not repeated here.

[0448] Based on the same inventive concept, embodiments of the present application provide a computer-readable storage medium, a computer program product comprising: computer program code, which, when executed on a computer, causes the computer to execute any of the trajectory optimization methods discussed above. Because the principles underlying the problem solved by the computer-readable storage medium are similar to those of the trajectory optimization method, the implementation of the computer-readable storage medium can be referred to as the implementation of the method, and any repetitions will not be repeated here.

[0449] Based on the same inventive concept, embodiments of the present application further provide a computer program product, comprising: computer program code, which, when executed on a computer, causes the computer to execute any of the trajectory optimization methods discussed above. Because the principles underlying the problems solved by the computer program products are similar to those of the trajectory optimization methods, the implementation of the computer program products can be referenced to the implementation of the methods, and any repetitions will not be repeated.

[0450] Those skilled in the art will appreciate that the embodiments of the present application can be provided as methods, systems, or computer program products. Therefore, the present application can adopt the form of a complete hardware embodiment, a complete software embodiment, or an embodiment in combination with software and hardware. Moreover, the present application can adopt the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to magnetic disk storage, CD-ROM, optical storage, etc.) that contain computer-usable program code.

[0451] The present application is described with reference to the flowcharts and / or block diagrams of the methods, devices (systems), and computer program products according to the present application. It should be understood that each process and / or block in the flowchart and / or block diagram, as well as the combination of processes and / or blocks in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowchart and / or block diagram. Figure 1 a process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0452] These computer program instructions may also be stored in a computer readable memory that can direct a computer or other programmable data processing device to work in a specific manner, so that the instructions stored in the computer readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 a process or multiple processes and / or boxes Figure 1 The function specified in one or more boxes.

[0453] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of user-operated steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing instructions for executing on the computer or other programmable device to implement the process. Figure 1 a process or multiple processes and / or boxes Figure 1 A step that specifies a function in one or more boxes.

[0454] Obviously, those skilled in the art may make various changes and modifications to this application without departing from the spirit and scope of this application. Thus, if these modifications and variations of this application fall within the scope of the claims of this application and their equivalents, this application is intended to include these modifications and variations.

Claims

1. A trajectory optimization method, characterized in that: The method includes: Obtaining a co-state initial value, setting an integration step and a duty cycle reference time, and determining an initial motion state of the detector at the earliest time in a set time period based on the co-state initial value; Set the earliest time as the first time, set the integration step as the first integration step, and the initial motion state as the first motion state, and trigger the following iteration: determining a first duty cycle function value based on the first time and the duty cycle reference time, and determining a second duty cycle function value based on the first time, the first integration step, and the duty cycle reference time; comparing a product of the first duty cycle function value and the second duty cycle function value with a set threshold to obtain a first comparison result, determining whether a thrust switching point exists for the detector within a target time interval based on the first comparison result, and determining a second integration step size based on the presence of the thrust switching point, wherein the target time interval is determined based on the first time and the first integration step size; Determining a second motion state of the detector at a second time based on the first motion state and the second integration step, using the second motion state as a new first motion state, the second integration step as a new first integration step, and the second time as a new first time, wherein the second time is determined based on the first time and the second integration step; If the second time does not reach the latest time of the set time period, the iteration is triggered again, otherwise the iteration ends; Based on the initial motion state and the second motion state in each iteration, a trajectory of the detector within a set time period is generated.

2. The method according to claim 1, wherein The determining, based on the first comparison result, whether the detector has a thrust switching point within a target time interval includes: If the first comparison result indicates that the product value is less than or equal to the set threshold, determining that the detector is at the thrust switching point within the target time interval; If the first comparison result indicates that the product value is greater than the set threshold, it is determined that the thrust switching point does not exist in the detector within the target time interval.

3. The method according to claim 2, wherein The determining of the second integration step size based on the existence of the thrust switching point comprises: If it is determined that the detector is at the thrust switching point within the target time interval, the first duty cycle function value is compared with the set threshold, and the second duty cycle function value is compared with the set threshold; if the first duty cycle function value is greater than or equal to the set threshold, and the second duty cycle function value is less than or equal to the set threshold, the second integration step size is determined based on the first duty cycle function value and the set thrust duration; if the first duty cycle function value is less than or equal to the set threshold, and the second duty cycle function value is greater than or equal to the set threshold, the second integration step size is determined based on the first duty cycle function value, the set thrust duration, and the set forced coasting duration; If it is determined that the thrust switching point does not exist in the detector within the target time interval, the first integration step is used as the second integration step.

4. The method according to claim 1, wherein Before determining a first duty cycle function value based on the first time and the duty cycle reference time, the method further includes: determining a first cosine value of a Sun-Earth-Probe SEP angle / a Sun-Probe-Earth SPE angle at the first time based on the first motion state; comparing the first cosine value with a set cosine threshold to obtain a second comparison result, and determining whether a thrust switching point exists for the detector within the target time interval based on the second comparison result; Based on the existence of the thrust switching point, it is determined whether to update the first integration step size.

5. The method according to claim 4, wherein The setting of the cosine threshold includes a first cosine threshold and a second cosine threshold, the first cosine threshold is greater than the second cosine threshold, and the determining, based on the second comparison result, whether the detector has a thrust switching point within the target time interval includes: If the second comparison result indicates that the first cosine value is greater than or equal to the first cosine threshold, predicting a second cosine value of the SEP angle / SPE angle at a third time based on the first motion state, the first cosine value, and the first integration step, and comparing the second cosine value with the first cosine threshold; if the second cosine value is less than the first cosine threshold, determining that the probe has the thrust switching point within the target time interval; and if the second cosine value is greater than or equal to the first cosine threshold, determining that the probe does not have the thrust switching point within the target time interval, wherein the third time is determined based on the first time and the first integration step; If the second comparison result indicates that the first cosine value is greater than the second cosine threshold, and the first cosine value is less than the first cosine threshold, determining whether the probe has the thrust switching point within the target time interval based on the first time and the set solar transit shutdown constraint interval; If the first comparison result indicates that the first cosine value is less than or equal to the second cosine threshold, it is determined that the thrust switching point does not exist in the detector within the target time interval.

6. The method according to claim 5, wherein The determining, based on the first time and the set solar eclipse shutdown constraint interval, whether the probe has the thrust switching point within the target time interval includes: If the first time is within the set solar transit shutdown constraint interval, determining that the probe does not have the thrust switching point within the target time interval; If the first time is not within the set solar transit shutdown constraint interval, a third motion state of the probe at a fourth time is determined based on the first motion state, the set solar transit early shutdown duration, and the first time. Based on the third motion state and the first integration step, a third cosine value of the SEP angle / SPE angle at a fifth time is predicted, and the third cosine value is compared with the first cosine threshold. If the third cosine value is greater than or equal to the first cosine threshold, it is determined that the probe has the thrust switching point within the target time interval. If the third cosine value is less than the first cosine threshold, it is determined that the probe does not have the thrust switching point within the target time interval. The fourth time is determined based on the set solar transit early shutdown duration and the first time, and the fifth time is determined based on the fourth time and the first integration step.

7. The method according to claim 6, wherein The determining whether to update the first integration step size based on the existence of the thrust switching point includes: If it is determined that the detector is at the thrust switching point within the target time interval, determining the thrust switching point using a set root method, and updating the first integration step size based on the thrust switching point and the first time; If it is determined that the detector does not have the thrust switching point within the target time interval, the first integration step size is not updated.

8. The method according to claim 7, wherein After determining the thrust switching point by using the set root method, the method further includes: If the thrust switching point is the time of sunrise, the duty cycle reference time is updated based on the thrust switching point and the set sunrise delay power-on time.

9. The method according to claim 1, wherein Before determining the first duty cycle function value based on the first time and the duty cycle reference time, the method further includes: determining a first available power of the detector at the first time based on the first motion state and a set detector system operating power; predicting a second available power of the detector at a third time based on the first motion state, the first available power, and the first integration step; comparing the first available power with a set power threshold to obtain a third comparison result, comparing the second available power with the set power threshold to obtain a fourth comparison result, and determining whether the thrust switching point exists within the target time interval based on the third and fourth comparison results; Based on the existence of the thrust switching point, it is determined whether to update the first integration step size.

10. The method according to claim 9, wherein The setting power threshold includes a first power threshold and a second power threshold, the first power threshold is greater than the second power threshold, and the determining, based on the third comparison result and the fourth comparison result, whether the detector has the thrust switching point within the target time interval includes: If the third comparison result indicates that the first available power is greater than or equal to the first power threshold, and the fourth comparison result indicates that the second available power is greater than or equal to the first power threshold, determining that the probe does not have the thrust switching point within the target time interval; If the third comparison result indicates that the first available power is greater than or equal to the first power threshold, and the fourth comparison result indicates that the second available power is less than the first power threshold, determining that the detector is at the thrust switching point within the target time interval; If the third comparison result indicates that the first available power is less than or equal to the second power threshold, and the fourth comparison result indicates that the second available power is less than or equal to the second power threshold, determining that the probe does not have the thrust switching point within the target time interval; If the third comparison result indicates that the first available power is less than or equal to the second power threshold, and the fourth comparison result indicates that the second available power is greater than the second power threshold, determining that the detector is at the thrust switching point within the target time interval; If the third comparison result indicates that the first available power is greater than the second power threshold and less than the first power threshold, the second available power is compared with the first available power and the third available power respectively to obtain a fifth comparison result. Based on the fifth comparison result, it is determined whether the probe has a thrust switching point within the target time interval, and the third available power is the available power of the probe at a sixth time, and the sixth time is earlier than the first time.

11. The method according to claim 10, wherein The determining, based on the fifth comparison result, whether the detector has a thrust switching point within the target time interval includes: If the fifth comparison result indicates that the second available power is greater than or equal to the first available power, and the second available power is less than or equal to the third available power, determining that the probe does not have the thrust switching point within the target time interval; If the fifth comparison result indicates that the second available power is greater than the third available power, determining that the detector is at the thrust switching point within the target time interval; If the fifth comparison result indicates that the second available power is less than the first available power, it is determined that the probe is at the thrust switching point within the target time interval.

12. The method according to any one of claims 10-11, characterized in that: The determining whether to update the first integration step size based on the existence of the thrust switching point includes: If it is determined that the detector is at the thrust switching point within the target time interval, determining the thrust switching point using a set root method, and updating the first integration step size based on the thrust switching point and the first time; If it is determined that the detector does not have the thrust switching point within the target time interval, the first integration step size is not updated.

13. The method according to claim 12, wherein: After updating the first integration step size based on the thrust switching point and the first time, the method further includes: determining a fourth motion state of the probe at the thrust switching point based on the first motion state and the first integration step size, and determining a fourth available power of the probe at the thrust switching point based on the fourth motion state and a set probe system operating power; determining a second co-state initial value of the detector at the thrust switching point based on a first co-state initial value of the detector at the first time and the first integration step; determining, based on a correspondence between each available power interval and a thruster specific impulse, a first thruster specific impulse corresponding to the available power interval to which the fourth available power belongs, and a second thruster specific impulse corresponding to the available power interval to which the first available power belongs; determining a first switching function value based on the fourth motion state, the second co-state initial value, and the first thruster specific impulse; determining a second switching function value based on the first motion state, the first co-state initial value, and the second thruster specific impulse; If the second switching function value is greater than a set homotopy parameter and the first switching function value is less than the set homotopy parameter, the duty cycle reference time is updated based on the thrust switching point.

14. The method according to claim 1, wherein The initial value of the co-state is the initial value of the co-state in the i-th iteration in the process of solving the detector fuel optimization problem using the homotopy algorithm, 0≤i≤N-1, N is the set number of homotopy steps, and the iteration is triggered again when the second time does not reach the latest time of the set time period. Otherwise, after the iteration ends, the method further includes: Based on the second motion state in the last iteration and the set motion state of the detector at the latest time in the set time period, the co-state initial value is updated, and the updated co-state initial value is used as the co-state initial value of the (i+1)th iteration.

15. A detector, characterized in that: The thruster comprises a thruster and a controller, wherein the controller is configured to execute the steps of the method according to any one of claims 1 to 14.

16. An electronic device, characterized in that: include: a memory for storing program instructions; A processor is configured to call program instructions stored in the memory and execute the steps of the method according to any one of claims 1 to 14 according to the obtained program instructions.

17. A computer-readable storage medium, characterized in that The computer-readable storage medium stores a computer program, wherein the computer program includes program instructions. When the program instructions are executed by a computer, the computer is caused to perform the method according to any one of claims 1 to 14.

Images