CN114872934B - Lunar surface aircraft track optimization method based on lossless convexity and Chebyshev orthogonal mating point method - Google Patents
Lunar surface aircraft track optimization method based on lossless convexity and Chebyshev orthogonal mating point method Download PDFInfo
- Publication number
- CN114872934B CN114872934B CN202210568718.5A CN202210568718A CN114872934B CN 114872934 B CN114872934 B CN 114872934B CN 202210568718 A CN202210568718 A CN 202210568718A CN 114872934 B CN114872934 B CN 114872934B
- Authority
- CN
- China
- Prior art keywords
- aircraft
- track
- lunar
- representing
- optimization
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 71
- 238000005457 optimization Methods 0.000 title claims abstract description 70
- 230000013011 mating Effects 0.000 title claims abstract description 21
- 230000001174 ascending effect Effects 0.000 claims abstract description 22
- 239000000446 fuel Substances 0.000 claims abstract description 14
- 230000000630 rising effect Effects 0.000 claims abstract description 12
- 230000001131 transforming effect Effects 0.000 claims abstract description 7
- PEDCQBHIVMGVHV-UHFFFAOYSA-N Glycerine Chemical compound OCC(O)CO PEDCQBHIVMGVHV-UHFFFAOYSA-N 0.000 claims description 28
- 239000011159 matrix material Substances 0.000 claims description 14
- 230000008569 process Effects 0.000 claims description 10
- 238000010606 normalization Methods 0.000 claims description 9
- 230000001133 acceleration Effects 0.000 claims description 8
- 230000009466 transformation Effects 0.000 claims description 6
- 238000004422 calculation algorithm Methods 0.000 description 18
- 238000009826 distribution Methods 0.000 description 5
- 230000008901 benefit Effects 0.000 description 4
- 230000008859 change Effects 0.000 description 4
- 238000006467 substitution reaction Methods 0.000 description 4
- 238000012546 transfer Methods 0.000 description 4
- 238000006243 chemical reaction Methods 0.000 description 3
- 230000000295 complement effect Effects 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 239000006185 dispersion Substances 0.000 description 3
- 230000009977 dual effect Effects 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 238000004088 simulation Methods 0.000 description 3
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 3
- 238000011217 control strategy Methods 0.000 description 2
- 238000013461 design Methods 0.000 description 2
- 238000002360 preparation method Methods 0.000 description 2
- 206010048669 Terminal state Diseases 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 238000009411 base construction Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 238000012938 design process Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 201000010099 disease Diseases 0.000 description 1
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 1
- 238000011010 flushing procedure Methods 0.000 description 1
- 230000007774 longterm Effects 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 230000035772 mutation Effects 0.000 description 1
- 238000009828 non-uniform distribution Methods 0.000 description 1
- 230000010355 oscillation Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000011946 reduction process Methods 0.000 description 1
- 230000035945 sensitivity Effects 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 238000003860 storage Methods 0.000 description 1
- 230000008685 targeting Effects 0.000 description 1
Classifications
-
- B—PERFORMING OPERATIONS; TRANSPORTING
- B64—AIRCRAFT; AVIATION; COSMONAUTICS
- B64G—COSMONAUTICS; VEHICLES OR EQUIPMENT THEREFOR
- B64G1/00—Cosmonautic vehicles
- B64G1/22—Parts of, or equipment specially adapted for fitting in or to, cosmonautic vehicles
- B64G1/24—Guiding or controlling apparatus, e.g. for attitude control
- B64G1/242—Orbits and trajectories
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T10/00—Road transport of goods or passengers
- Y02T10/10—Internal combustion engine [ICE] based vehicles
- Y02T10/40—Engine management systems
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Chemical & Material Sciences (AREA)
- Combustion & Propulsion (AREA)
- Radar, Positioning & Navigation (AREA)
- Aviation & Aerospace Engineering (AREA)
- Feedback Control In General (AREA)
Abstract
A lunar surface large-scale aircraft track optimization method based on lossless convexity and Chebyshev orthogonal mating point method relates to the technical field of lunar surface aircraft flight track optimization and is used for solving the problem of power rising and power falling two-stage track optimization of lunar surface aircraft. The technical key points of the invention include: track optimization is carried out on the ascending track and the descending track of the aircraft in the lunar surface flight stage: constructing a lunar surface manned aircraft track optimization model which comprises a dynamics model equation, constraint conditions and an objective function; transforming the trajectory optimization model, including transforming constraint conditions and dynamic model equations; further discretizing a kinetic model equation by using a Chebyshev orthogonal mating point method; and carrying out iterative solution on the transformed track optimization model to obtain the state quantity and the control quantity of the ascending track or the descending track. The invention realizes the track planning of the ascending section and the descending section of the lunar surface aircraft aiming at the optimal fuel, and has universality.
Description
Technical Field
The invention relates to the technical field of lunar surface aircraft flight trajectory optimization, in particular to a lunar surface large-range aircraft trajectory optimization method based on a lossless convexity and chebyshev orthogonal mating point method.
Background
In order to better develop and utilize resources of the moon distributed in various areas, the moon base may be built in different areas of the moon, and the transportation of large-range and long-distance personnel and materials between the moon is inevitably required, so that the existing conventional transportation modes of lunar vehicles and the like cannot meet the requirements. In addition, the rescue work of personnel, the detection work in a large range and the like all need to be used for completing the carrier which takes off on the moon surface, flies according to a planned route and lands autonomously, and the transportation cost can be greatly reduced due to repeated utilization. However, because of the special environments such as overweight and no atmosphere in the moon, the fuel consumed by the aircraft in the two flight stages of power rising and power braking can be greatly increased, so that the transportation cost is greatly increased, and a trajectory optimization algorithm which is suitable for the moon environment and takes the fuel optimization as an index is required to reasonably plan the two flight stages of the aircraft.
In order to solve the above problems, a lunar surface aircraft designed based on task demands of water ice exploitation and transfer is proposed in document [1], and the aircraft is analyzed and designed from subsystems such as a loading mechanism, a manned system, an avionics system, a power system and the like, and main flight tasks and schemes thereof are proposed and planned. According to the designed improved trajectory, the aircraft can transport two astronauts and 170 kg of goods for a large-scale transportation task within 6 minutes, and the voyage can reach 57km. However, it is disadvantageous in that neither the flight time nor the flight distance can meet the requirement for wide range transfer. The literature [2] proposes an on-line guidance algorithm of the air suction hypersonic aircraft with process constraint and terminal state constraint based on a pseudo-spectrum method, which can remarkably reduce the on-line calculation amount and improve the algorithm performance. The method does not need the design process of an inner loop feedback controller, and can continuously generate open loop suboptimal control on line. However, the control time of the smaller flight mission range is shorter, and the control strategy can be obtained online, but the control strategy cannot be directly applied to the track optimization of the lunar surface aircraft due to different constraint conditions. A idea of solving the problem of Mars landing by means of convex optimization is given in document [3 ]. The essence of this idea is to highlight the whole optimization problem. Mapping the original strong nonlinear constraint condition into convex constraint by introducing a relaxation variable method, and proving that the original property of the constraint condition is not changed; the index function and the constraint condition are linearized through linearization methods such as Taylor expansion, discrete and the like, so that the overall optimization problem is converted into a convex optimization problem, and the problem can be solved through a faster convex optimization algorithm such as an interior point method and the like. However, it is less versatile and requires higher requirements for model simplification and convexity.
In summary, it is an important issue to design a fast trajectory optimization algorithm that enables a lunar surface aircraft to perform safe power up and power down and is premised on energy optimization.
Disclosure of Invention
In view of the problems, the invention provides a lunar surface large-scale aircraft track optimization method based on a lossless convexity and Chebyshev orthogonal mating point method, which is used for solving the track optimization problem of a lunar surface aircraft in two stages of power rising and power falling.
A lunar surface large-scale aircraft track optimization method based on lossless convexity and Chebyshev orthogonal mating point method comprises the following steps of track optimization for an ascending section track and a descending section track of an aircraft in a lunar surface flight stage:
step one, constructing a lunar surface manned aircraft track optimization model; the method specifically comprises the following steps:
establishing a kinetic model equation of the lunar surface manned aircraft under a lunar inertial spherical coordinate system;
establishing constraint conditions and an objective function of track optimization by taking fuel optimization as a target, wherein the constraint conditions comprise an initial moment state, a terminal moment state, an engine thrust amplitude, a thrust angle amplitude and a safety distance;
step two, transforming the track optimization model to obtain a transformed track optimization model; the method specifically comprises the following steps:
based on a lossless convexity theory, introducing a relaxation variable to transform constraint conditions in a track optimization model;
introducing a time control quantity to transform the dynamics model equation;
discretizing the transformed dynamic model equation by using a Chebyshev orthogonal mating point method;
and thirdly, carrying out iterative solution on the transformed track optimization model to obtain a position sequence, a speed sequence, an engine thrust amplitude sequence and a pitch angle and yaw angle sequence of the aircraft of all position points in the ascending section track or the descending section track.
Further, in the first step, the lunar inertial spherical coordinate system is a spherical coordinate system established by taking a lunar centroid as an origin of coordinates, and points in the spherical coordinate system are represented by (r, θ, phi), wherein r represents a distance from the point to the lunar centroid, and θ represents an included angle between a position vector of the point and a lunar north pole, namely a polar angle; phi represents the angle between the projection of the position vector of the point on the equatorial plane of the moon and the initial meridian plane of the moon, namely the azimuth angle.
Further, in the first step, the kinetic model equation is normalized, and the normalized kinetic model equation is established as follows:
1) Ascending section:
2) A descending section:
in the formula, v r 、v θ 、v φ Respectively representing the radial, tangential and normal speeds of the aircraft; μ represents a lunar gravitational constant; t represents the engine thrust amplitude; alpha and beta respectively represent the included angles between the thrust direction and the radial direction and the normal direction of the aircraft; m represents the aircraft mass; v e Representing engine specific impulse;respectively represent r, theta, phi and v r 、v θ 、v φ Derivative of m with respect to time t;
wherein the rising sectionm 0 Representing the aircraft mass at the initial moment, R F Representing the distance from the end point of the rising segment to the center of mass of the moon, v F Representing the aircraft speed corresponding to the end point of the ascending section, T max Representing a maximum thrust amplitude of the engine;
descending sectionR 0 Representing the distance from the beginning of the descent segment to the center of mass of the moon, v 0 Representing the aircraft speed corresponding to the start point of the descent segment;
the constraint conditions are as follows:
defining an initial time t 0 Status: x (t) 0 )＝[r 0 ；θ 0 ；φ 0 ；v r0 ；v θ0 ；v φ0 ；m 0 ]Wherein r is 0 、θ 0 、φ 0 、v r0 、v θ0 、v φ0 、m 0 Respectively represent r, theta, phi and v r 、v θ 、v φ A value at an initial time;
defining a terminal time t F Status: x (t) F )＝[r F ；θ F ；φ F ；v rF ；v θF ；v φF ；m F ]Wherein r is F 、θ F 、φ F 、v rF 、v θF 、v φF 、m F Respectively represent r, theta, phi and v r 、v θ 、v φ A value at the terminal time;
setting the range of the thrust amplitude of the engine, namely meeting the following conditions:T min representing a minimum thrust amplitude of the engine;
setting the thrust angle amplitude range of the aircraft, namely, meeting the following conditions: alpha is more than or equal to 0 and less than or equal to pi; beta is more than or equal to 0 and less than or equal to pi;
setting a safe distance range: r is more than or equal to (R+H)/R 0 Wherein R represents the average radius of the moon, and H represents a preset safety height threshold;
the objective function is: j=m 0 -m F The method comprises the steps of carrying out a first treatment on the surface of the Wherein m is F Representing the aircraft mass at the terminal moment.
Further, constraint conditions after the relaxation variable is introduced to be transformed in the second step are as follows:
wherein z is 0 ＝lnm 0 ，z F ＝lnm F Z= lnm (τ), τ representing the normalization time; σ represents a relaxation variable; u (u) r 、u θ 、u φ The radial, tangential and normal accelerations of the aircraft in the spherical coordinate system are indicated, respectively.
Further, the kinetic model equation after the time control quantity is introduced and transformed in the second step is as follows:
wherein u is t Representing the control quantity of the terminal time, which is the derivative of the time t to the normalized time tau, namelyRepresenting the derivative of z with respect to the normalized time τ; in the rising section, c=c 1 In the descending segment, c=c 2 ；
For the transformed kinetic model equation, the state quantity is expressed as:
x(τ)＝[r；θ；φ；v r ；v θ ；v φ ；z；t]
the control amount is expressed as follows:
p(τ)＝[u r ；u θ ；u φ ；σ；u t ]
the state quantity and the control quantity are uniformly defined as a variable matrix, expressed as:
y(τ)＝[x(τ)；p(τ)]
expressing the transformed kinetic model equation as the derivative of x (τ)The function for y (τ), expressed as:
at the initial time tau 0 The state quantity and the control quantity are subjected to first-order Taylor expansion, and the finally obtained transformed kinetic model equation is as follows:
wherein y (τ) 0 )＝[x(τ 0 )；p(τ 0 )]Representing an initial variable matrix.
Further, in the second step, the discrete process of the transformed dynamic model equation by using the Chebyshev orthogonal mating point method comprises time domain transformation and mating point discrete; the discretized kinetic model equation is:
wherein X (τ) i ) Representing the discrete state quantity; y (τ) k ) Representing a discrete variable matrix consisting of state quantities and control quantities; d (D) k,i Representing the derivative value of the Lagrangian basis function at the discrete points after the discretization; i and k each represent the location of an interpolation node.
And further, in the third step, the transformed track optimization model is subjected to iterative solution by using an original-dual interior point method.
The beneficial technical effects of the invention are as follows: the invention realizes the track planning of the ascending section and the descending section of the lunar surface aircraft aiming at the fuel optimization, improves the precision of the track optimization algorithm and the smoothness of the track curve, expands the application range of the convex optimization algorithm frame and enhances the universality of the algorithm.
Drawings
The invention may be better understood by reference to the following description taken in conjunction with the accompanying drawings, which are included to provide a further illustration of the preferred embodiments of the invention and to explain the principles and advantages of the invention, together with the detailed description below.
FIG. 1 is a flow chart of a method for optimizing lunar surface aircraft trajectories based on lossless projection and Chebyshev's orthogonal mating point method in accordance with an embodiment of the present invention;
FIG. 2 is a schematic diagram of a lunar inertial Cartesian coordinate system and a spherical coordinate system according to an embodiment of the present invention;
FIG. 3 is a schematic layout of a lunar surface manned aircraft in an embodiment of the present invention;
FIG. 4 is a graph showing an example of a change in position in a lunar inertial Cartesian coordinate system according to an embodiment of the present invention;
FIG. 5 is a graph illustrating an example of velocity change in a lunar inertial Cartesian coordinate system according to an embodiment of the present invention;
FIG. 6 is an exemplary graph of thrust sequence in a lunar inertial Cartesian coordinate system according to an embodiment of the present invention.
Detailed Description
In order that those skilled in the art will better understand the present invention, exemplary embodiments or examples of the present invention will be described below with reference to the accompanying drawings. It is apparent that the described embodiments or examples are only implementations or examples of a part of the invention, not all. All other embodiments or examples, which may be made by one of ordinary skill in the art without undue burden, are intended to be within the scope of the present invention based on the embodiments or examples herein.
In order to meet the requirements of large-scale personnel and cargo transfer tasks in the later stage of moon base construction, the lunar surface aircraft has the characteristics of reusability and low transportation cost. Since the moon is free of the atmosphere, it is not possible to provide lift for the aircraft when ascending and braking force when descending, and these two phases of the mission consume a large amount of fuel. Therefore, the invention provides a lunar surface aircraft track optimization method, which optimizes the track of the ascending section and the track of the descending section of the aircraft on the premise of optimizing energy, and aims to realize the ascending and the descending of the aircraft on the premise of consuming less fuel as possible.
The invention combines the lossless salifying idea with the orthogonal point by using the Chebyshev orthogonal polynomial, and the mixed algorithm idea is to perform proper salifying treatment before the orthogonal point, so that the convergence of the original algorithm can be improved to a certain extent, the algorithm solving speed can be accelerated, the application range can be expanded, and the algorithm universality can be enhanced.
The method of the invention optimizes the ascending section track and the descending section track of the aircraft in the flying stage of the moon surface according to the following steps:
step one, constructing a lunar surface manned aircraft track optimization model; the method specifically comprises the following steps: establishing a kinetic model equation of the lunar surface manned aircraft under a lunar inertial spherical coordinate system; establishing constraint conditions and an objective function of track optimization by taking fuel optimization as a target, wherein the constraint conditions comprise an initial moment state, a terminal moment state, an engine thrust amplitude, a thrust angle amplitude and a safety distance;
step two, transforming the track optimization model to obtain a transformed track optimization model; the method specifically comprises the following steps: based on a lossless convexity theory, introducing a relaxation variable to transform constraint conditions in a track optimization model; introducing a time control quantity to transform the dynamics model equation; discretizing the transformed dynamic model equation by using a Chebyshev orthogonal mating point method;
and thirdly, carrying out iterative solution on the transformed track optimization model to obtain a position sequence, a speed sequence, an engine thrust amplitude sequence and a pitch angle and yaw angle sequence of the aircraft of all position points in the ascending section track or the descending section track.
The embodiment of the invention provides a lunar surface large-scale aircraft track optimization method based on a lossless convexity and Chebyshev orthogonal mating point method, which comprises the following steps as shown in figure 1:
step S1: constructing an initial optimization problem of the lunar surface manned aircraft; namely, establishing a dynamic model, constraint conditions and index functions;
in order to facilitate description of movement of a lunar surface aircraft in a wide range of transfers, according to embodiments of the present invention, related concepts are defined and a lunar inertial Cartesian coordinate system and a spherical coordinate system are established as follows:
(1) Lunar north: the moon rotation axis intersects with the moon along the rotation direction.
(2) The lunar equator: a great circle intersecting the moon through the center of mass of the moon and perpendicular to the plane of the rotation axis.
(3) Lunar origin: in long-term observation, the average position of the projected point of the centroid of the earth on the lunar surface.
(4) Lunar primary meridian plane: a large circle where a plane formed by the lunar origin, the lunar centroid and the lunar north intersects the moon.
(5) Lunar inertial cartesian coordinate system: a Cartesian coordinate system is established by taking a lunar centroid as a coordinate origin, a lunar equatorial plane is taken as a basic plane, a Z axis points to a lunar north pole perpendicular to the lunar equatorial plane, an X axis points to a lunar longitude origin direction in the lunar equatorial plane, and X, Y, Z three axes accord with a right-handed screw rule to form a right-handed system. The coordinate diagram is shown in fig. 2.
(6) Lunar inertial sphere coordinate system: and establishing a spherical coordinate system by taking the lunar centroid as a coordinate origin. Any point in the lunar inertial Cartesian coordinate system can be represented by three ordered real numbers (r, θ, φ). Wherein r is the distance from the point to the center of mass of the moon; θ is the angle between the bit vector of the point and the north pole of the moon, i.e. the polar angle; phi is the angle between the projection of the position vector of the point on the lunar equatorial plane and the lunar primary meridian plane, namely the azimuth angle. The coordinate diagram is shown in fig. 2.
In order to facilitate the description of the angular variation of the aircraft during the flight, a kinetic model of the lunar surface manned aircraft is built under a lunar inertial spherical coordinate system:
wherein r, theta and phi respectively represent the distance from the lunar center, the polar angle and the azimuth angle of the aircraft under a lunar inertial spherical coordinate system; v r 、v θ 、v φ Respectively representing the radial, tangential and normal speeds of the aircraft; alpha and beta respectively represent the included angles between the thrust direction and the radial direction and the normal direction of the aircraft; mu is the lunar gravitational constant; t is the thrust amplitude of the engine; m is the mass of the aircraft; v e The engine is specific flushing;representing the derivative of each variable with respect to time, respectively.
In addition, based on the kinetic equation, other variables and constants are defined as follows:
(1) Defining a state quantity function matrix: x (t) = [ r (t); θ (t); phi (t); v r (t)；v θ (t)；v φ (t)；m(t)]。
(2) Defining a control quantity function matrix: u (T) = [ T (T); alpha (t); beta (t) ].
(3) Defining an initial time state: x (t) 0 )＝[r 0 ；θ 0 ；φ 0 ；v r0 ；v θ0 ；v φ0 ；m 0 ]Wherein t is 0 Indicating the initial moment of flight, r 0 、θ 0 、φ 0 、v r0 、v θ0 、v φ0 、m 0 Respectively represent each variable at the initial time t 0 Values at that time.
(4) Defining a terminal moment state: x (t) F )＝[r F ；θ F ；φ F ；v rF ；v θF ；v φF ；m F ]Wherein t is F Indicating the end time, i.e. the end time, r of the flight F 、θ F 、φ F 、v rF 、v θF 、v φF 、m F Respectively represent the initial variables at the terminal time t F Values at that time.
In order to improve the convergence property of the solution, the condition of occurrence of a disease state during the solution is avoided, a proper constant is selected as a standard (called a normalization factor), and each variable in the problem is normalized and unified into a dimensionless variable.
In the normalization process of the ascending section dynamics equation, selecting the maximum thrust of the engineForce amplitude T max Distance R from the end point of the rising section to the center of the moon F Initial mass m of aircraft 0 Aircraft speed v corresponding to the end of the rise F As a normalization factor; in the normalization process of the descending-section dynamics equation, selecting the maximum thrust T of the engine max Distance R from the beginning of descent segment to moon center 0 Initial mass m of aircraft 0 Aircraft speed v corresponding to the start of descent 0 As a normalization factor. The equation is obtained after normalization as follows:
the constant factor c is defined here and is not actually physically significant, but is merely convenient for further deduction below. Wherein the rising section is optimized by takingWhen optimizing descending segment, taking +.>
After the dynamics equation is normalized, the relevant constraint conditions of the trajectory optimization problem can be obtained as follows:
(1) Variable engine thrust amplitude:T min representing a minimum thrust amplitude of the engine;
(2) Thrust angle amplitude: alpha is more than or equal to 0 and less than or equal to pi; beta is more than or equal to 0 and less than or equal to pi;
(3) Safety distance constraint conditions: r is more than or equal to (R+H)/R 0 Wherein R and H represent the average radius of the moon and the set safety height, respectively. The safe distance is the height of the aircraft from the moon surface, i.e. the flight path is required to be at least not lower than a height, which may be typically set to 2km, to avoid collision with mountains.
Targeting fuel optimisation, establishing an objective function or indexA function. Using the total mass m of the aircraft at the initial moment 0 With the total mass m of the aircraft at the moment of termination F The difference in (2) represents the mass of fuel consumed, expressed in mathematical form: j=m 0 -m F 。
In view of the above description, the original trajectory optimization problem can be initially established: and under the constraint conditions of satisfying an aircraft dynamics equation, an initial moment state, a terminal moment state, an engine thrust amplitude, a thrust angle amplitude and a safety distance, solving a state quantity and control quantity function matrix, so that the fuel is optimal, namely the consumed fuel is minimum. Can be expressed in the mathematical form as follows:
Minimize:J＝m 0 -m F
step S2: variable substitution and relaxation: introducing a relaxation variable and transforming constraint conditions and a kinetic equation;
according to the embodiment of the invention, due to the existence of complex nonlinear constraint, a great amount of operations in the optimizing process of the track optimizing problem consume iteration at a local extreme point. To enhance the reliability and convergence of the algorithm, the nonlinearity of the algorithm is reduced to a certain extent by variable substitution. The introduced variable substitution relationship is as follows:
where z has no special physical meaning, for convenience of derivation only, z= lnm (τ), τ representing the normalization time; u is the total acceleration generated by the thrust of the engine, u r 、u θ 、u φ The magnitudes of the acceleration of the aircraft in radial, tangential and normal directions in the spherical coordinate system are shown, respectively.
The kinetic equation at this time is obtained by variable substitution:
it can be seen that at this time there is no multiplication or division relation between the control quantity and the state quantity, but a new constraint is introduced: the module of the acceleration vector formed by the accelerations in three directions should be equal to the total acceleration generated by the thrust, and the mathematical expression is as follows:
in order to reduce the strong non-convexity that it brings as much as possible, while increasing the algorithm convergence, a relaxation variable σ is introduced here according to the lossless convexity theory. The aim is to relax the total acceleration u, changing its constraint into a range instead of a certain value.
Definition of the definitionMeanwhile, the thrust amplitude constraint at this time becomes:
T min e -z ≤σ≤T max e -z
the second-order taylor expansion can be obtained at the initial moment:
the constraint conditions at this time can be obtained after finishing and are expressed as follows:
wherein z is 0 ＝lnm 0 、z F ＝lnm F The value of the variable z at the initial time and the final time are indicated respectively.
Step S3: introducing a time control amount; namely, transforming the dynamic model;
according to the embodiment of the invention, the problem of terminal time uncertainty is solved for the two-stage track optimization problem of the power rising section and the power falling section, and although one terminal time can be given manually according to experience, the autonomy of the algorithm is greatly reduced. Therefore, in the embodiment of the invention, the terminal time is also used as one of the control quantities, and the time is also used as a state variable and is listed in a dynamics equation.
Defining normalized time tau and terminal time control quantity u t And defining the time relationship as follows:
t＝u t τ,τ∈[0,1]
the derivative of the original left side of the kinetic equation on the time t is converted into the derivative of the normalized time tau, and a new state variable equation is added to obtain the kinetic equation as follows:
the defined state quantity at this time is:
x(τ)＝[r；θ；φ；v r ；v θ ；v φ ；z；t]
the control amount is defined as follows:
p(τ)＝[u r ；u θ ；u φ ；σ；u t ]
the state quantity and the control quantity are uniformly defined as a variable matrix:
y(τ)＝[x(τ)；p(τ)]
the initial variable matrix may be expressed as:
y(τ 0 )＝[x(τ 0 )；p(τ 0 )]
in summary, the kinetic equation can be expressed asWith respect to the function of y (τ), the mathematical form is as follows:
in order to reduce the nonlinearity, the dynamics equation is subjected to first-order taylor expansion at the initial state and the initial control quantity, and the following form is obtained:
wherein,,the transpose of the gradient matrix at the initial time of each variable representing the kinetic equation. />
Step S4: performing dispersion of a kinetic equation by using a Chebyshev orthogonal mating point method;
according to the embodiment of the invention, the core idea of the Chebyshev orthogonal fitting point method is to use the zero point of the Chebyshev polynomial as an interpolation node, use the Lagrange polynomial as an interpolation basis function, discrete a continuous kinetic equation and convert the original problem of optimizing continuous variables into the problem of optimizing discrete variables on the fitting points.
The advantage of using chebyshev polynomial zero interpolation over other node interpolation methods is: the method can comprise two points of an initial state and a termination state, discrete results are not required to be amplified, and each variable cannot generate larger mutation near two ends; because the distribution of the Chebyshev polynomial zero points on [ -1,1] is not uniform, the distribution is more dense at two ends and more sparse in the middle, the distribution of the interpolation nodes can avoid the Dragon phenomenon to a great extent; the extremum of the Lagrangian interpolation remainder can be minimized, so that the interpolated result is closest to the original function value.
The Chebyshev orthogonal mating point method mainly comprises two steps of time domain transformation and mating point dispersion respectively:
(1) Time domain transform
Since chebyshev polynomials satisfy the orthogonal relationship in the time domain [ -1,1], and there is a recurrence relationship as follows:
T 0 (τ)＝1,T 1 (τ)＝τ,
T n+1 (τ)＝2τT n (τ)-T n-1 (τ)
the polynomial function T (τ) here is a fixed expression of the chebyshev polynomial of the first type, independent of the thrust function described above, and is therefore described in order to avoid ambiguity.
In order to take advantage of the non-uniform distribution of zero points on [ -1,1] of the chebyshev polynomial, the time domain of the original problem needs to be replaced by [ -1,1] as follows:
(2) Distribution points and dispersion
The fitting point is the zero point on the Chebyshev polynomial [ -1,1 ]. When the fitting point is infinite, the final solution will fully converge to the original optimization problem, but interpolation will cause larger errors due to the Dragon phenomenon. Even though chebyshev polynomials are less likely than other methods, it is still not preferable to choose a particularly large number of nodes, which would not only severely affect the speed but also produce severe oscillations making the result infeasible. Therefore, the proper number of distribution points should be selected, and the occurrence of the Dragon phenomenon is avoided as much as possible. The suitability of the number of the preparation points is related to the problem of actual processing, and the number of the preparation points can be simply adjusted according to the precision effect. General deduction is made here, no specific assignment is made to the number of points, the number of dots is denoted as n.
First, the first class of Chebyshev polynomials is found in [ -1,1]N zeros above, denoted as tau i I=1, 2, …, n, while satisfying T (τ i )＝0,i＝1,2,...,n。
Then, the state variable and the control variable are discretized by utilizing Lagrange interpolation, and the expression is as follows:
wherein,,
note that the kinetic equationTo the left of (a) is to derive the state variable from the time factor τ, and therefore also to interpolate the derivative of the state variable, expressed as follows: />
The final solution is the state and control variables at the discrete points, and all that is required is the constraint of the kinetic equation at the discrete points, so although the interpolation function is continuous, only the function values at the discrete points on the interpolation function are considered, expressed mathematically as follows:
wherein k=2,..n; i and k each represent the location of an interpolation node.
The discretized total variable matrix may also be defined based on discretized control variables and state variables:
Y(τ k )＝[X(τ k )；P(τ k )]
additionally provided with differential matrixThe expression meaning is: after deriving the Lagrange interpolation basis function, the derivative value of the Lagrange basis function at the discrete point is obtained, and the mathematical expression is as follows:
where i=1,..n, k=2,..n,
in summary, the kinetic constraints after discretization can be obtained as follows:
step S5: the final optimization problem is obtained and iterative solution is carried out;
according to the embodiment of the invention, the optimization problem is finally obtained through the corresponding processes of the above conversion and the conversion as follows:
Minimize：J＝m 0 -m F
Subject t0：
the constraint condition and the index function form conform to a more typical second-order cone planning problem in the convex optimization problem, and the second-order cone planning problem can be solved iteratively by using an original-dual interior point method in a traditional convex optimization algorithm.
The basic idea of the original-dual interior point method is to iteratively represent the inequality constraint as an equality constraint problem. Firstly, converting an original problem into a dual problem, and then solving an original dual gap to obtain a complementary condition. The problem is converted into an optimal solution for solving the complementary condition, the solution can be carried out by a Newton method or other numerical optimization methods, and finally the original variable and the dual variable are continuously updated in the process of solving the complementary condition optimal solution in an iteration mode until the iteration stop condition is met.
Further, taking a power reduction process of a conceptual lunar surface aircraft as an example, the technical effect of the invention is verified in a simulation way.
The aircraft generally comprises a nose, a manned cabin and a power cabin, and the full load weight is designed to be 12t. 4.7 tons of storage tanks of LH2 and LO2 with diameters of 1.2 meters and 0.6 meters are placed around the manned cabin, an auxiliary engine is arranged around the manned cabin for fixed-point landing, and a posture-adjusting thruster is placed around the manned cabin. The main engine is arranged at the bottom of the power cabin and mainly provides power for the ascending section and also provides decelerating thrust for the descending section. The general layout is schematically shown in fig. 3.
The aircraft transport mission is designed as: personnel return to the vicinity of the lunar base from the area where the lunar water ice is exploited, the vicinity of the lunar base rich in water ice (89.9 degrees S,180 degrees W) is selected as a starting point, and a lunar base candidate address (2.5 degrees N,86.5 degrees E) is selected as a target point. The whole process is divided into a vertical ascending section, a power ascending section, a track section, a power descending section and a vertical descending section. The highest height of the whole process is planned to be not more than 15km, and the safety distance is set to be 2km in height.
The invention is applicable to two parts of a power ascending section and a power descending section, so the power descending section is taken as an example for simulation explanation, and relevant parameters are selected as shown in the following table 1:
table 1 simulation parameter table
| Initial azimuth angle phi 0 | 70.32° |
| Initial polar angle θ 0 | 25° |
| Target azimuth angle phi F | 86.5° |
| Target polar angle theta F | 2.5° |
| Initial moon center distance r 0 /km | 1758 |
| Target moon center distance r F /km | 1740 |
| Initial velocity vector v 0 /m·s -1 | [-762.8,869.4,-1204.6] |
| Target velocity vector v F /m·s -1 | [0,0,0] |
| Initial mass m 0 /kg | 12000 |
The number of discrete points is selected to be n=200, modeling, transformation and solving are carried out according to the optimization flow described by the embodiment of the invention, and finally, a position sequence, a speed sequence, a thrust amplitude sequence and a pitch angle and yaw angle sequence of the aircraft on 200 discrete points are obtained. Through the conversion, the position change, the speed change and the thrust sequence of the aircraft in three directions under the lunar inertial Cartesian coordinate system can be obtained, as shown in figures 4-6.
It should be noted that, the position curve, the velocity curve and the thrust sequence obtained by the method are all pre-trajectory designs performed by the aircraft before the actual flight mission, and can be used as a nominal trajectory curve or a nominal thrust sequence in actual control. In specific use, different nominal values may be selected according to the type of the controller, such as: if position control is actually used, the space coordinates of the aircraft are monitored in real time and correspondingly adjusted by taking the position curve obtained by solving as a reference or standard; if the engine thrust and the attitude of the aircraft are controlled in practice, the main engine can be controlled to match the solved thrust sequence, and the attitude-adjusting engine is controlled to match the solved pitch angle and yaw angle sequences. The result of the embodiment of the invention is used as a standard to control, so that the flight track with optimal fuel in theory can be obtained.
While the invention has been described with respect to a limited number of embodiments, those skilled in the art, having benefit of the above description, will appreciate that other embodiments are contemplated within the scope of the invention as described herein. The disclosure of the present invention is intended to be illustrative, but not limiting, of the scope of the invention, which is defined by the appended claims.
The documents cited in the present invention are as follows:
[1]D Akin."Project ASHLAIN:A Lunar Flying Vehicle for Rapid Universal Surface Access."Aiaa Space Conference&Exposition,2013.
[2]Da Zhang，Lei Liu，Yongji Wang.On-line Ascent Phase Trajectory Optimal Guidance Algorithm based on Pseudo-spectral Method and Sensitivity Updates.Journal of Navigation,68(6),pp.1056-1074，2015.
[3]Acikmese B,Ploen S R.Convex programming approach to powered descent guidance for mars landing.Journal of Guidance,Control and Dynamics,30(5),1353-1366，2007.
Claims (5)
1. a lunar surface large-scale aircraft track optimization method based on lossless convexity and Chebyshev orthogonal mating point method is characterized in that the track optimization is carried out on the ascending section track and the descending section track of an aircraft in the flying stage of the lunar surface according to the following steps:
step one, constructing a lunar surface manned aircraft track optimization model; the method specifically comprises the following steps:
establishing a kinetic model equation of the lunar surface manned aircraft under a lunar inertial spherical coordinate system; the lunar inertial spherical coordinate system is a spherical coordinate system established by taking a lunar centroid as a coordinate origin, points in the spherical coordinate system are represented by (r, theta, phi), wherein r represents the distance from the points to the lunar centroid, and theta represents the included angle between the position vector of the points and the north pole of the lunar, namely the polar angle; phi represents the included angle between the projection of the position vector of the point on the lunar equatorial plane and the lunar primary meridian plane, namely the azimuth angle; normalizing the kinetic model equation, and establishing the normalized kinetic model equation as follows:
1) Ascending section:
2) A descending section:
in the formula, v r 、v θ 、v φ Respectively representing the radial, tangential and normal speeds of the aircraft; μ represents a lunar gravitational constant; t represents the engine thrust amplitude; alpha and beta respectively represent the included angles between the thrust direction and the radial direction and the normal direction of the aircraft; m represents the aircraft mass; v e Representing engine specific impulse;respectively represent r, theta, phi and v r 、v θ 、v φ Derivative of m with respect to time t;
wherein the rising sectionm 0 Representing the aircraft mass at the initial moment, R F Representing the distance from the end point of the rising segment to the center of mass of the moon, v F Representing the aircraft speed corresponding to the end point of the ascending section, T max Representing a maximum thrust amplitude of the engine;
descending sectionR 0 Representing the distance from the beginning of the descent segment to the center of mass of the moon, v 0 Representing the aircraft speed corresponding to the start point of the descent segment;
establishing constraint conditions and an objective function of track optimization by taking fuel optimization as a target, wherein the constraint conditions comprise an initial moment state, a terminal moment state, an engine thrust amplitude, a thrust angle amplitude and a safety distance; the constraint conditions are as follows:
defining an initial time t 0 Status: x (t) 0 )＝[r 0 ；θ 0 ；φ 0 ；v r0 ；v θ0 ；v φ0 ；m 0 ]Wherein r is 0 、θ 0 、φ 0 、v r0 、v θ0 、v φ0 、m 0 Respectively represent r, theta, phi and v r 、v θ 、v φ A value of m at the initial time;
defining a terminal time t F Status: x (t) F )＝[r F ；θ F ；φ F ；v rF ；v θF ；v φF ；m F ]Wherein r is F 、θ F 、φ F 、v rF 、v θF 、v φF 、m F Respectively represent r, theta, phi and v r 、v θ 、v φ The value of m at the terminal moment;
setting the range of the thrust amplitude of the engine, namely meeting the following conditions:T min representing a minimum thrust amplitude of the engine;
setting the thrust angle amplitude range of the aircraft, namely, meeting the following conditions: alpha is more than or equal to 0 and less than or equal to pi; beta is more than or equal to 0 and less than or equal to pi;
setting a safe distance range: r is more than or equal to (R+H)/R 0 Wherein R represents the average radius of the moon, and H represents a preset safety height threshold;
the objective function is: j=m 0 -m F The method comprises the steps of carrying out a first treatment on the surface of the Wherein m is F Representing the aircraft mass at the moment of the terminal;
step two, transforming the track optimization model to obtain a transformed track optimization model; the method specifically comprises the following steps:
based on a lossless convexity theory, introducing a relaxation variable to transform constraint conditions in a track optimization model;
introducing a time control quantity to transform the dynamics model equation;
discretizing the transformed dynamic model equation by using a Chebyshev orthogonal mating point method;
and thirdly, carrying out iterative solution on the transformed track optimization model to obtain a position sequence, a speed sequence, an engine thrust amplitude sequence and a pitch angle and yaw angle sequence of the aircraft of all position points in the ascending section track or the descending section track.
2. The lunar surface large-scale aircraft trajectory optimization method based on lossless protrusion and chebyshev orthogonal mating point method according to claim 1, wherein the constraint condition after the relaxation variable is introduced for transformation in the second step is as follows:
wherein z is 0 ＝lnm 0 ，z F ＝lnm F Z= lnm (τ), τ representing the normalization time; σ represents a relaxation variable; u (u) r 、u θ 、u φ The radial, tangential and normal accelerations of the aircraft in the spherical coordinate system are indicated, respectively.
3. The lunar surface large-scale aircraft trajectory optimization method based on lossless protrusion and chebyshev orthogonal mating point method according to claim 2, wherein a dynamics model equation after time control quantity is introduced for transformation in the second step is as follows:
wherein u is t Representing the control quantity of the terminal time, which is the derivative of the time t to the normalized time tau, namely Representing the derivative of z with respect to the normalized time τ; in the rising section, c=c 1 In the descending segment, c=c 2 ；
For the transformed kinetic model equation, the state quantity is expressed as:
x(τ)＝[r；θ；φ；v r ；v θ ；v φ ；z；t]
the control amount is expressed as follows:
p(τ)＝[u r ；u θ ；u φ ；σ；u t ]
the state quantity and the control quantity are uniformly defined as a variable matrix, expressed as:
y(τ)＝[x(τ)；p(τ)]
expressing the transformed kinetic model equation as the derivative of x (τ)The function for y (τ), expressed as:
at the initial time tau 0 The state quantity and the control quantity are subjected to first-order Taylor expansion, and the finally obtained transformed kinetic model equation is as follows:
wherein y (τ) 0 )＝[x(τ 0 )；p(τ 0 )]Representing an initial variable matrix.
4. The lunar surface large-scale aircraft trajectory optimization method based on lossless protrusion and chebyshev orthogonal fitting point method according to claim 3, wherein the process of discretizing the transformed kinetic model equation by using chebyshev orthogonal fitting point method in the second step comprises time domain transformation and fitting point discretization; the discretized kinetic model equation is:
wherein X (τ) i ) Representing the discrete state quantity; y (τ) k ) Representing a discrete variable matrix consisting of state quantities and control quantities; d (D) k,i Representing the derivative value of the Lagrangian basis function at the discrete points after the discretization; i and k each represent the location of an interpolation node.
5. The lunar surface large-scale aircraft trajectory optimization method based on lossless protrusion and chebyshev orthogonal mating point method according to claim 4, wherein in the third step, the transformed trajectory optimization model is iteratively solved by using an original-dual interior point method.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210568718.5A CN114872934B (en) | 2022-05-24 | 2022-05-24 | Lunar surface aircraft track optimization method based on lossless convexity and Chebyshev orthogonal mating point method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210568718.5A CN114872934B (en) | 2022-05-24 | 2022-05-24 | Lunar surface aircraft track optimization method based on lossless convexity and Chebyshev orthogonal mating point method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114872934A CN114872934A (en) | 2022-08-09 |
| CN114872934B true CN114872934B (en) | 2023-10-13 |
Family
ID=82676787
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210568718.5A Active CN114872934B (en) | 2022-05-24 | 2022-05-24 | Lunar surface aircraft track optimization method based on lossless convexity and Chebyshev orthogonal mating point method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114872934B (en) |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109976154A (en) * | 2019-03-04 | 2019-07-05 | 北京理工大学 | A kind of aerial vehicle trajectory optimization method based on chaos multinomial and the convex optimization of sequence |
| CN111931131A (en) * | 2020-08-12 | 2020-11-13 | 北京航空航天大学 | Online trajectory planning method and system for power landing segment of planetary probe |
| CN112051854A (en) * | 2020-09-23 | 2020-12-08 | 北京理工大学 | Rapid planning method for optimal trajectory of lunar soft landing considering complex constraints |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8489260B2 (en) * | 2008-12-16 | 2013-07-16 | California Institute Of Technology | Method and apparatus for powered descent guidance |
| US20120265374A1 (en) * | 2011-04-15 | 2012-10-18 | Thomas Edward Yochum | Aircraft vertical trajectory optimization method |
| US8655589B2 (en) * | 2012-01-25 | 2014-02-18 | Mitsubishi Electric Research Laboratories, Inc. | System and method for controlling motion of spacecrafts |
-
2022
- 2022-05-24 CN CN202210568718.5A patent/CN114872934B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109976154A (en) * | 2019-03-04 | 2019-07-05 | 北京理工大学 | A kind of aerial vehicle trajectory optimization method based on chaos multinomial and the convex optimization of sequence |
| CN111931131A (en) * | 2020-08-12 | 2020-11-13 | 北京航空航天大学 | Online trajectory planning method and system for power landing segment of planetary probe |
| CN112051854A (en) * | 2020-09-23 | 2020-12-08 | 北京理工大学 | Rapid planning method for optimal trajectory of lunar soft landing considering complex constraints |
Non-Patent Citations (1)
| Title |
|---|
| 基于凸优化和切比雪夫伪谱法的再入轨迹优化;陈嘉澍;;电子技术与软件工程(16);第71-75页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114872934A (en) | 2022-08-09 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Wang et al. | Improved sequential convex programming algorithms for entry trajectory optimization | |
| Morgan et al. | Model predictive control of swarms of spacecraft using sequential convex programming | |
| CN106054604B (en) | Reentry vehicle robust optimal method of guidance based on Model Predictive Control Theory | |
| Zhu et al. | Geometric optimal control and applications to aerospace | |
| Jiang et al. | Integrated guidance for Mars entry and powered descent using reinforcement learning and pseudospectral method | |
| CN112051854B (en) | Rapid planning method for optimal trajectory of lunar soft landing considering complex constraints | |
| CN110750850A (en) | Three-dimensional profile optimization design method, system and medium under strong constraint complex task condition | |
| Wang et al. | Predictor-corrector entry guidance with waypoint and no-fly zone constraints | |
| Polsgrove et al. | Human Mars Entry, Descent, and Landing Architecture Study: Rigid Decelerators | |
| Wang et al. | Optimal trajectory-tracking guidance for reusable launch vehicle based on adaptive dynamic programming | |
| Bae et al. | Convex optimization-based entry guidance for spaceplane | |
| Wang et al. | Autonomous rendezvous guidance via deep reinforcement learning | |
| CN114872934B (en) | Lunar surface aircraft track optimization method based on lossless convexity and Chebyshev orthogonal mating point method | |
| Chen et al. | Steady Glide Dynamics and Guidance of Hypersonic Vehicle | |
| Benedikter et al. | Stochastic control of launch vehicle upper stage with minimum-variance splash-down | |
| Restrepo et al. | Structured adaptive model inversion controller for Mars atmospheric flight | |
| Haghparast et al. | A cubature Kalman filter for parameter identification and output-feedback attitude control of liquid-propellant satellites considering fuel sloshing effects | |
| Wang | Maximum-normal-load entry trajectory optimization for hypersonic glide vehicles | |
| Zhou et al. | Ascent trajectory optimization for air‐breathing vehicles in consideration of launch window | |
| Dilão et al. | Dynamic trajectory control of gliders | |
| CN112325711A (en) | Carrier rocket orbit height keeping control method based on balanced flight theory | |
| Bestaoui et al. | Geometry of Translational Trajectories for an Autonomous Vehicle With Wind Effect | |
| Radhakrishnan et al. | 6D trajectory, guidance and control development for air-breathing phase of reusable launch vehicle | |
| Li et al. | Improved guidance algorithm considering terminal attitude constraints for spacecrafts via optimal control theory | |
| Ashraf et al. | Iterative guidance scheme for satellite launch vehicles |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |