CN114872934A - Lunar surface large-range aircraft trajectory optimization method based on lossless convexity and Chebyshev orthogonal fitting method - Google Patents
Lunar surface large-range aircraft trajectory optimization method based on lossless convexity and Chebyshev orthogonal fitting method Download PDFInfo
- Publication number
- CN114872934A CN114872934A CN202210568718.5A CN202210568718A CN114872934A CN 114872934 A CN114872934 A CN 114872934A CN 202210568718 A CN202210568718 A CN 202210568718A CN 114872934 A CN114872934 A CN 114872934A
- Authority
- CN
- China
- Prior art keywords
- aircraft
- lunar
- optimization
- representing
- trajectory
- 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.)
- Granted
Links
Images
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-range aircraft trajectory optimization method based on lossless convexity and a Chebyshev orthogonal fitting method relates to the technical field of lunar surface aircraft flight trajectory optimization and is used for solving the two-stage trajectory optimization problems of power rise and power fall of a lunar surface aircraft. The technical points of the invention comprise: performing trajectory optimization on the ascending section trajectory and the descending section trajectory of the aircraft in the lunar surface flight phase: constructing a lunar manned vehicle trajectory optimization model which comprises a dynamic model equation, constraint conditions and an objective function; transforming the track optimization model, including transforming constraint conditions and a dynamic model equation; further discretizing a dynamic model equation by utilizing a Chebyshev orthogonal fitting method; and carrying out iterative solution on the transformed track optimization model to obtain the state quantity and the control quantity of the ascending section track or the descending section track. The method realizes the trajectory planning of the ascending section and the descending section of the lunar surface aircraft aiming at fuel optimization, and has universality.
Description
Technical Field
The invention relates to the technical field of optimization of flight trajectories of lunar surface aircrafts, in particular to a lunar surface large-range aircraft trajectory optimization method based on lossless convexity and Chebyshev orthogonal fitting method.
Background
In order to better develop and utilize resources of the moon distributed in various regions, a moon base may be built in different regions of the moon, and large-range and long-distance personnel and material transportation between the months is inevitable, at which time, the conventional carrying modes of the moon rover and the like cannot meet the requirements. In addition, rescue work of personnel, large-scale detection work and the like all need to use a carrier which can finish lunar takeoff, fly according to a planned route and independently land, and the transportation cost can be greatly reduced by recycling. However, due to the heavy mass of the aircraft and the special environment such as no atmosphere in the moon, the fuel consumed in the two flight phases of power-up and power-braking is greatly increased, which leads to a great increase in transportation cost.
In order to solve the above problems, document [1] proposes a lunar aircraft designed based on the mission requirements of water ice mining and transfer, which is analyzed and designed by subsystems such as a load mechanism, a manned system, an avionic system and a power system, and proposes a main flight mission and scheme thereof. According to the designed improved ballistic trajectory, the aircraft can transport two astronauts and 170 kilograms of goods in 6 minutes for a large-scale transportation task, and the range can reach 57 km. However, it is disadvantageous in that neither the time of flight nor the distance of flight can meet the requirements for large-scale transfers. The document [2] provides an air-breathing hypersonic aircraft online guidance algorithm with process constraint and terminal state constraint based on a pseudo-spectrum method, which can obviously reduce online calculated amount and improve 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 flight mission range is small, the control time is short, and although the control strategy can be obtained online, the control strategy cannot be directly applied to the trajectory optimization of the lunar aircraft due to the difference of the constraint conditions. The idea of solving the Mars landing problem by using a convex optimization method is given in the document [3 ]. The essence of this idea is to highlight the entire optimization problem. The original strong nonlinear constraint condition is mapped into convex constraint by introducing a relaxation variable method, and the original property of the convex constraint condition is not changed by proving; the index function and the constraint condition are linearized through Taylor expansion, dispersion and other linearization methods, so that the overall optimization problem is converted into a convex optimization problem, and the solution can be performed through a faster convex optimization algorithm such as an interior point method. However, the universality is low, and the requirements on model simplification and convexity are high.
In conclusion, it is an important problem to design a fast trajectory optimization algorithm for safely ascending and descending power of a lunar aircraft and based on the premise of optimal energy.
Disclosure of Invention
In view of the above problems, the invention provides a lunar large-scale aircraft trajectory optimization method based on lossless convexity and Chebyshev orthogonal fitting method, which is used for solving the trajectory optimization problem of two stages of power rising and power falling of a lunar aircraft.
A lunar surface large-range aircraft trajectory optimization method based on lossless convexity and Chebyshev orthogonal fitting method is characterized in that trajectory optimization is carried out on an ascent trajectory and a descent trajectory of an aircraft in a lunar surface flight phase according to the following steps:
step one, constructing a lunar manned aircraft trajectory optimization model; the method specifically comprises the following steps:
establishing a dynamic model equation of the lunar manned aircraft under a lunar inertial spherical coordinate system;
establishing a constraint condition and an objective function of track optimization by taking fuel optimization as a target, wherein the constraint condition comprises 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 convex theory, introducing a relaxation variable to transform a constraint condition in the track optimization model;
introducing a time control quantity to transform the kinetic model equation;
discretizing the transformed dynamic model equation by utilizing a Chebyshev orthogonal fitting 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 at 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 a coordinate origin, and points in the spherical coordinate system are represented by (r, θ, φ), wherein r represents a distance from a point to the lunar centroid, and θ represents an included angle between a vector of a point and a north pole of the moon, namely a polar angle; phi represents the included angle between the projection of the point position vector on the equatorial plane of the moon and the meridian plane of the moon.
Further, the dynamic model equation is normalized in the first step, and the normalized dynamic 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; mu denotes the gravitational force of the moonA constant value; t represents the thrust amplitude of the engine; alpha and beta respectively represent included angles between the thrust direction of the aircraft and the radial direction and the normal direction; m represents the aircraft mass; v. of e Representing the engine specific impulse;
respectively represent r, theta, phi and v r 、v θ 、v φ The derivative of m with respect to time t;
wherein the rising section
m 0 Representing the aircraft mass at the initial moment, R F Representing the distance, v, from the end of the rise to the center of mass of the moon F Indicating the aircraft speed, T, corresponding to the end of the approach max Representing a maximum thrust magnitude of the engine;
descending section
R 0 Representing the distance, v, from the beginning of the descent segment to the center of mass of the moon 0 Representing the speed of the aircraft corresponding to the starting point of the descending section;
the constraint conditions are as follows:
defining an initial time t 0 The state is as follows: 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 The state is as follows: 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 φ At the time of terminalA value;
setting the range of the thrust amplitude of the engine, namely satisfying the following conditions:
T min representing a minimum thrust magnitude of the engine;
and setting the amplitude range of the thrust angle of the aircraft, namely meeting the following requirements: 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 not less than (R + H)/R 0 Wherein R represents the average radius of the moon, and H represents a preset safe height threshold value;
the objective function is: j ═ m 0 -m F (ii) a Wherein m is F Representing the aircraft mass at the terminal moment.
Further, the constraint conditions after the relaxation variables are introduced in the step two and transformed are as follows:
in the formula, z 0 ＝lnm 0 ，z F ＝lnm F Lnm (τ), τ denotes normalized time; σ represents a relaxation variable; u. of r 、u θ 、u φ Respectively representing the radial, tangential and normal accelerations of the aircraft in a spherical coordinate system.
Further, the dynamic model equation after the time control quantity is introduced and transformed in the step two is as follows:
in the formula u t Representing the terminal time control quantity, is the derivative of time t over normalized time τ, i.e.
To representThe derivative of z with respect to normalized time τ; in the ascending section, c ═ c 1 In the descending section, 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 quantity is expressed as:
p(τ)＝[u r ；u θ ；u φ ；σ；u t ]
the state quantity and the control quantity are uniformly defined as a variable matrix, and are expressed as follows:
y(τ)＝[x(τ)；p(τ)]
expressing the transformed kinetic model equation as the derivative of x (τ)
The function for y (τ) is expressed as:
at an initial time τ 0 The state quantity and the control quantity are subjected to first-order Taylor expansion, and the finally obtained transformed dynamic model equation is in the form as follows:
in the formula, y (τ) 0 )＝[x(τ 0 )；p(τ 0 )]And represents an initial variable matrix.
Further, the discretization process of the transformed dynamic model equation by utilizing the Chebyshev orthogonal fitting method in the second step comprises time domain transformation and fitting point discretization; the discretized kinetic model equation is as follows:
in the formula, X (τ) i ) Representing the state quantity after dispersion; y (τ) k ) Representing a discrete variable matrix consisting of state quantity and control quantity; d k,i Representing derivative values of Lagrangian basis functions at the discrete points after the dispersion; i and k each represent the position of an interpolation node.
Further, in the third step, the transformed trajectory optimization model is iteratively solved by using an original-dual interior point method.
The beneficial technical effects of the invention are as follows: the method realizes the trajectory planning of the ascending section and the descending section of the lunar aircraft aiming at fuel optimization, improves the precision of the trajectory optimization algorithm and the smoothness of a trajectory curve, expands the application range of a convex optimization algorithm frame, and enhances the universality of the algorithm.
Drawings
The present invention may be better understood by reference to the following description taken in conjunction with the accompanying drawings, which are incorporated in and form a part of this specification, and which are used to further illustrate preferred embodiments of the present invention and to explain the principles and advantages of the present invention.
FIG. 1 is a flow chart of a lunar aircraft trajectory optimization method based on lossless convexity and Chebyshev orthogonal fitting methods in an embodiment of the present invention;
FIG. 2 is a schematic diagram of a lunar inertial Cartesian coordinate system and a spherical coordinate system in accordance with an embodiment of the present invention;
FIG. 3 is a schematic layout of a lunar manned aircraft in an embodiment of the invention;
FIG. 4 is an exemplary diagram of a change in position in a lunar inertial Cartesian coordinate system in accordance with an embodiment of the present invention;
FIG. 5 is an exemplary graph of velocity variations in a lunar inertial Cartesian coordinate system in accordance with an embodiment of the present invention;
FIG. 6 is an exemplary thrust sequence under a lunar inertial Cartesian coordinate system in an embodiment of the present invention.
Detailed Description
In order that those skilled in the art will better understand the disclosure, exemplary embodiments or examples of the disclosure are described below with reference to the accompanying drawings. It is obvious that the described embodiments or examples are only some, but not all embodiments or examples of the invention. All other embodiments or examples obtained by a person of ordinary skill in the art based on the embodiments or examples of the present invention without any creative effort shall fall within the protection scope of the present invention.
In order to meet the requirements of large-scale personnel and goods transfer tasks in the later construction stage of a lunar base, the lunar surface aircraft has the characteristics of repeated use and low transportation cost. Since the moon is not free of atmosphere and cannot provide lift for the aircraft when ascending or provide braking when descending, these two phases of the mission will consume a significant amount of fuel. Therefore, the invention provides a lunar aircraft trajectory optimization method, which optimizes the trajectories of the ascending section and the descending section of an aircraft on the premise of energy optimization, and aims to realize the ascending and the descending of the aircraft on the premise of consuming as little fuel as possible.
The invention combines the idea of lossless convexity with the positive mating point by utilizing the Chebyshev orthogonal polynomial for use, and the mixed algorithm idea is to perform proper convexity treatment before the orthogonal mating point, so that the convergence of the original algorithm can be improved to a certain extent, the algorithm solving speed is accelerated, the application range is expanded, and the algorithm universality is enhanced.
The method of the invention optimizes the trajectories of the ascending section and the descending section of the aircraft in the lunar surface flight phase according to the following steps:
step one, constructing a lunar manned aircraft trajectory optimization model; the method specifically comprises the following steps: establishing a dynamic model equation of the lunar manned aircraft under a lunar inertial spherical coordinate system; establishing a constraint condition and an objective function of track optimization by taking fuel optimization as a target, wherein the constraint condition comprises 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 convex theory, introducing a relaxation variable to transform a constraint condition in the track optimization model; introducing a time control quantity to transform the kinetic model equation; discretizing the transformed dynamic model equation by utilizing a Chebyshev orthogonal fitting 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-range aircraft trajectory optimization method based on lossless convexity and Chebyshev orthogonal fitting method, as shown in figure 1, the optimization method comprises the following steps:
step S1: constructing an initial optimization problem of the lunar manned aircraft; namely, establishing a dynamic model, a constraint condition and an index function;
according to the embodiment of the invention, in order to facilitate the description of the movement of the lunar vehicle in the large-scale transfer, the relevant concepts are defined and the lunar inertial cartesian coordinate system and the spherical coordinate system are respectively established as follows:
(1) lunar north pole: the intersection point of the moon rotation axis and the moon along the rotation direction.
(2) Equator of the moon: and a great circle which passes through the center of mass of the moon and is perpendicular to the plane of the rotation shaft and intersected with the moon.
(3) Origin of lunar longitude: in long term observations, the mean position of the projected point of the earth's centroid on the lunar surface.
(4) Lunar book junior noon surface: and a great circle which is formed by a plane formed by the longitude origin of the moon, the center of mass of the moon and the north pole of the moon and intersects the moon.
(5) Lunar inertial cartesian coordinate system: a Cartesian coordinate system is established by taking the center of mass of the moon as a coordinate origin, the equatorial plane of the moon is a basic plane, the Z axis is perpendicular to the equatorial plane of the moon and points to the north pole of the moon, the X axis points to the longitude origin direction of the moon in the equatorial plane of the moon, and the X, Y, Z three axes accord with a right-hand screw rule to form a right-hand system. The schematic diagram of coordinates is shown in fig. 2.
(6) Lunar inertial sphere coordinate system: and establishing a spherical coordinate system by taking the center of mass of the moon as the origin of coordinates. Any point in the lunar inertial cartesian coordinate system can be represented by the three ordered real numbers (r, θ, φ). Wherein r is the distance from the point to the center of mass of the moon; theta is an included angle between the potential vector of the point and the north pole of the moon, namely a polar angle; phi is the included angle between the projection of the position vector of the point on the equatorial plane of the moon and the meridian plane of the moon. The schematic diagram of coordinates is shown in fig. 2.
In order to facilitate the description of the angle change of the aircraft in the flight process, a dynamic model of the lunar manned aircraft is established under a lunar inertial spherical coordinate system:
in the formula, r, theta and phi respectively represent the distance from the moon center, the polar angle and the azimuth angle of the aircraft under a lunar inertial spherical coordinate system; v. of r 、v θ 、v φ Respectively representing the radial, tangential and normal speeds of the aircraft; alpha and beta respectively represent included angles between the thrust direction of the aircraft and the radial direction and the normal direction; mu is a constant value of lunar gravity; t is the thrust amplitude of the engine; m is the mass of the aircraft; v. of e Is the engine specific impulse;
respectively, the derivative of each variable with respect to time.
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. of r (t)；v θ (t)；v φ (t)；m(t)]。
(2) Defining a control quantity function matrix: u (t) ═ t (t); α (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 representing the variables at an initial time t 0 The value of time.
(4) Defining a terminal time state: x (t) F )＝[r F ；θ F ；φ F ；v rF ；v θF ；v φF ；m F ]Wherein, t F Representing the terminal, i.e. end, time of flight, r F 、θ F 、φ F 、v rF 、v θF 、v φF 、m F Respectively representing the initial variables at terminal time t F The value of time.
In order to improve the convergence property of the solution and avoid the occurrence of a pathological condition during solution, a proper constant is selected as a standard (called a normalization factor), and all variables in the problem are normalized and unified into dimensionless variables.
Selecting the maximum thrust amplitude T of the engine in the normalization process of the ascending section dynamic equation max Distance R between the end point of the ascending section and the center of the moon F Initial mass m of the aircraft 0 And the speed v of the aircraft corresponding to the end of the ascent F As a normalization factor; in the normalization process of the descending section dynamic equation, the maximum thrust T of the engine is selected max Distance R between the starting point of the descending section and the moon center 0 Initial mass m of the aircraft 0 Aircraft speed v corresponding to the beginning of descent segment 0 As a normalization factor. The following equation is obtained after normalization:
the constant factor c is defined here and has no practical physical significance, but only to facilitate further derivation below. Wherein when the rising section is optimized, take
When optimizing the descending section, take
After the kinetic 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 magnitude 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 not less than (R + H)/R 0 Wherein R and H represent the moon average radius and the set safety height, respectively. The safety distance is the height of the aircraft from the surface of the moon, i.e. the flight trajectory is required to be at least not lower than an altitude to avoid collision with a mountain, which may be set to 2km in general.
An objective function, or referred to as an index function, is established with the fuel optimization as the target. Using total mass m of aircraft at initial moment 0 Total mass m of aircraft at terminal moment F The difference in (b) represents the fuel mass consumed, expressed in mathematical form: j ═ m 0 -m F 。
Combining the above description, the original trajectory optimization problem can be initially established: under the constraint conditions of meeting the aircraft dynamics equation, the initial moment state, the terminal moment state, the engine thrust amplitude, the thrust angle amplitude and the safety distance, the state quantity and control quantity function matrix is solved, so that the fuel is optimal, namely the consumed fuel is minimum. Can be expressed in 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 the complex nonlinear constraint, a large amount of operations are consumed in the optimization process of the track optimization problem in the iteration of the local extreme point. In order to enhance the reliability and convergence of the algorithm, the nonlinearity is firstly reduced to a certain extent by variable substitution. Introducing a variable substitution relationship as follows:
where z has no particular physical meaning, simply for convenience of derivation, and is lnm (τ), τ denotes normalized time; u is the magnitude of the total acceleration generated by the thrust of the engine, u r 、u θ 、u φ The radial, tangential and normal acceleration of the aircraft in a spherical coordinate system are respectively shown.
Through variable substitution, the kinetic equation at the moment is obtained as follows:
it can be seen that at this time, there is no multiplication or division relationship between the control quantity and the state quantity, but a new constraint condition is introduced: the modulus of an acceleration vector formed by the accelerations in three directions is equal to the magnitude of the total acceleration generated by the thrust, and the mathematical expression is as follows:
in order to reduce the strong non-convexity brought by the method as much as possible and improve the convergence of the algorithm, a relaxation variable sigma is introduced according to a lossless convexity theory. The objective is to relax the total acceleration u, constraining it to a range rather than a deterministic value.
Definition of
At the same time, the thrust amplitude constraint at this time will become:
T min e -z ≤σ≤T max e -z
performing a second order taylor expansion at the initial time may result in:
the constraints at this point in time after finishing are expressed as follows:
wherein z is 0 ＝lnm 0 、z F ＝lnm F Respectively, as the value of the variable z at the initial and terminal times.
Step S3: introducing a time control quantity; namely, the dynamic model is transformed;
according to the embodiment of the invention, the two-stage track optimization problem of the power ascending section and the power descending section belongs to the problem of uncertain terminal time, and although one terminal time can be manually given according to experience, the autonomy of the algorithm is greatly reduced by doing so. Therefore, in the embodiment of the invention, the terminal time is also taken as one of the control quantities, and the time is also taken as a state variable and is listed in the kinetic equation.
Defining a normalized time tau and a terminal time control u t And defining the time relationship as follows:
t＝u t τ,τ∈[0,1]
the original derivation of the left side of the kinetic equation to the time t is converted into the derivation of the normalized time tau, and a new state variable equation is added to obtain the following kinetic equation:
the state quantities are defined here as:
x(τ)＝[r；θ；φ；v r ；v θ ；v φ ；z；t]
defining the control quantity as:
p(τ)＝[u r ；u θ ；u φ ；σ；u t ]
the state quantity and the control quantity are uniformly defined as a variable matrix as follows:
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 as
With respect to the function of y (τ), the mathematical form is as follows:
to reduce the nonlinearity, the kinetic equation is further subjected to a first-order Taylor expansion at the initial state and the initial controlled variable, and the result is given by the following form:
wherein the content of the first and second substances,
the transpose of the gradient matrix of the variables representing the kinetic equation at the initial instant.
Step S4: performing dispersion of a kinetic equation by utilizing a Chebyshev orthogonal fitting method;
according to the embodiment of the invention, the core idea of the Chebyshev orthogonal fitting method is that a zero point of a Chebyshev polynomial is used as an interpolation node, a Lagrangian polynomial is used as an interpolation basis function, a continuous kinetic equation is dispersed, and the original problem of optimizing continuous variables is converted into the problem of optimizing discrete variables on fitting points.
Compared with the interpolation method of other nodes, the zero-point interpolation method adopting the Chebyshev polynomial has the advantages that: the method can comprise two points of an initial state and a termination state, the discrete result does not need to be amplified, and each variable cannot have larger mutation near the two ends; because the distribution of Chebyshev polynomial zero points on [ -1,1] is not uniform, the Chebyshev polynomial zero points are distributed more densely at two ends and sparsely in the middle, the distribution of the interpolation nodes can avoid the occurrence of the grid phenomenon to a great extent; the extreme value of the Lagrange interpolation remainder can be minimized, so that the interpolated result is closest to the original function value.
The Chebyshev orthogonal fitting method mainly comprises two steps, namely time domain transformation and fitting point dispersion:
(1) time domain transformation
Since the Chebyshev polynomial satisfies 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 (τ) is a fixed expression of the first chebyshev polynomial, and is not related to the thrust function, so as to avoid ambiguity.
In order to utilize the characteristic of uneven distribution of zero point of Chebyshev polynomial on [ -1,1], the time domain of the original problem needs to be replaced by the following equation of [ -1,1 ]:
(2) coordination and dispersion
And the point matching is to solve the zero point on the Chebyshev polynomial [ -1,1 ]. When the collocation point is infinite, the final solution is completely converged to the original optimization problem, but the interpolation causes larger errors due to the Runge phenomenon. Even though chebyshev polynomials are less prone to this than other methods, it is not desirable to choose a particularly large number of nodes, otherwise not only the speed is severely affected but also severe oscillations are produced making the result unfeasible. Therefore, the proper number of distribution points is selected to avoid the occurrence of the dragon lattice phenomenon as much as possible. Whether the number of the configuration points is proper or not is related to the problem of actual processing, and simple adjustment can be performed according to the precision effect. Here, general derivation is performed, and no specific assignment is made to the number of nodes, so the number of nodes is represented as n.
Firstly, a first Chebyshev polynomial is required to be in [ -1,1 [ -1 [ ]]N zeros, denoted as τ i I is 1,2, …, n, and satisfies T (tau) i )＝0,i＝1,2,...,n。
And then, dispersing the state variables and the control variables by utilizing Lagrange interpolation, wherein the expression is as follows:
wherein the content of the first and second substances,
note that the kinetic equation
To the left of (1) is to differentiate the state variable by a time factor τ, and therefore also to interpolate and discretize the derivative of the state variable, the expression is as follows:
the final solution is the state variables and the control variables at discrete points, and all that is needed is the constraints of the kinetic equations at discrete points, so although the interpolation function is continuous, only the function values at discrete points on the interpolation function are considered, and the mathematical expression is as follows:
wherein, k is 2.., n; i and k each represent the position of an interpolation node.
The discretized total variable matrix can also be defined according to the discretized control variables and state variables:
Y(τ k )＝[X(τ k )；P(τ k )]
additionally setting a differential matrix
The expression means: after derivation of the Lagrange interpolation basis function, the derived value of the Lagrange basis function at the discrete point is obtained, and the mathematical expression is as follows:
wherein, i is 1, n, k is 2, n,
in summary, the following dynamic constraints can be obtained:
step S5: obtaining a final optimization problem and carrying out iterative solution;
according to the embodiment of the invention, through the corresponding processes of transformation and conversion, the optimization problem is finally obtained as follows:
Minimize:J＝m 0 -m F
the form of the constraint condition and the index function accords with a typical second-order cone planning problem in a 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 primitive-dual interior point method is to iteratively represent inequality constraints as an equality constraint problem. Firstly, an original problem is converted into a dual problem, and then an original dual gap is solved to obtain a complementary condition. The problem is converted into an optimal solution for solving complementary conditions, the solution can be carried out through a Newton method or other numerical optimization methods, and finally the original variables and the dual variables are continuously updated in the process of iteratively solving the optimal solution for the complementary conditions until the iteration stop conditions are met.
Further, by taking the power descending process of the conceptual lunar aircraft as an example, the technical effect of the invention is simulated and verified.
The aircraft is composed of a nose, a manned cabin and a power cabin, and the full load weight is designed to be 12 t. 4.7 tons of storage tanks of LH2 and LO2 with the diameters of 1.2 meters and 0.6 meter are arranged around the manned cabin, meanwhile, an auxiliary engine is arranged around the manned cabin for fixed-point landing, and a posture adjusting thruster is arranged around the auxiliary engine. The main engine is arranged at the bottom of the power cabin and mainly provides power for the ascending section and also provides deceleration thrust for the descending section. The general layout is schematically shown in fig. 3.
The aircraft transportation mission is designed as follows: the method realizes that personnel return to the vicinity of the lunar base from a region where the lunar water ice is exploited, the vicinity (89.9 degrees S and 180 degrees W) of the lunar polar region where the water ice is rich is selected as a starting point, and the candidate address (2.5 degrees N and 86.5 degrees E) of the lunar base 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 maximum height of the whole process is planned to be not more than 15km, and the safety distance is set to be 2 km.
The invention is suitable for two parts of a power ascending section and a power descending section, so that the power descending section is taken as an example for simulation description, and relevant parameters are selected as shown in the following table 1:
TABLE 1 simulation parameters Table
| Initial azimuth angle phi 0 | 70.32° |
| Initial polar angle θ 0 | 25° |
| Target azimuth angle phi F | 86.5° |
| Target polar angle θ 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 the discrete points is selected to be n, which is 200, modeling, conversion and solving are carried out according to the optimization process of 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 the 200 discrete points are obtained. Through the transformation, the position change, the speed change and the thrust sequence in three directions of the aircraft under the lunar inertia Cartesian coordinate system can be obtained, as shown in FIGS. 4-6.
It should be noted that the position curve, the speed curve and the thrust sequence obtained by the present invention are designed for the pre-trajectory of the aircraft before the actual flight mission, and these curves can be used as the nominal trajectory curve or the nominal thrust sequence in the actual control. When the controller is used specifically, different nominal values can be selected according to different types of controllers, such as: if position control is actually used, real-time monitoring is carried out on the space coordinate of the aircraft and corresponding adjustment is carried out by taking the position curve obtained by solving as a reference or standard; if the thrust of the engine and the attitude of the aircraft are actually controlled, the main engine can be controlled to match and solve the obtained thrust sequence, and the attitude adjusting engine is controlled to match with the obtained pitch angle and yaw angle sequence. The result of the embodiment of the invention is used as a standard for control, and the theoretically optimal flight path of the fuel 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 this description, will appreciate that other embodiments can be devised which do not depart from the scope of the invention as described herein. The present invention has been disclosed in an illustrative rather than a restrictive sense, and the scope of the present invention 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 (7)
1. a lunar surface large-range aircraft trajectory optimization method based on lossless convexity and Chebyshev orthogonal fitting method is characterized in that trajectory optimization is carried out on an ascending section trajectory and a descending section trajectory of an aircraft in a lunar surface flying stage according to the following steps:
step one, constructing a lunar manned aircraft trajectory optimization model; the method specifically comprises the following steps:
establishing a dynamic model equation of the lunar manned aircraft under a lunar inertial spherical coordinate system;
establishing a constraint condition and an objective function of track optimization by taking fuel optimization as a target, wherein the constraint condition comprises 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 convex theory, introducing a relaxation variable to transform a constraint condition in the track optimization model;
introducing a time control quantity to transform the kinetic model equation;
discretizing the transformed dynamic model equation by utilizing a Chebyshev orthogonal fitting 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 method for optimizing the trajectories of the lunar surface wide-range aircraft based on the lossless convexity and the Chebyshev orthogonal fitting method as claimed in claim 1, wherein 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, θ, φ), wherein r represents a distance from a point to the lunar centroid, and θ represents an included angle between a position vector of a point and a north pole of the moon, namely a polar angle; phi represents the included angle between the projection of the point position vector on the equatorial plane of the moon and the meridian plane of the moon.
3. The method for optimizing the lunar surface large-range aircraft trajectory based on the lossless convexity and Chebyshev orthogonal fitting method according to claim 2, wherein in the first step, the dynamic model equation is normalized, and the normalized dynamic 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 moon gravity constant; t represents the thrust amplitude of the engine; alpha and beta respectively represent included angles between the thrust direction of the aircraft and the radial direction and the normal direction; m represents the aircraft mass; v. of e Representing the engine specific impulse;
respectively represent r, theta, phi and v r 、v θ 、v φ The derivative of m with respect to time t;
wherein the rising section
m 0 Representing the aircraft mass at the initial moment, R F Representing the distance, v, from the end of the rise to the center of mass of the moon F Indicating aircraft speed corresponding to terminal point of ascentDegree, T max Representing a maximum thrust magnitude of the engine;
descending section
R 0 Representing the distance, v, from the beginning of the descent segment to the center of mass of the moon 0 Representing the speed of the aircraft corresponding to the starting point of the descending section;
the constraint conditions are as follows:
defining an initial time t 0 The state is as follows: 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 The state is as follows: 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 φ Value at terminal time;
setting the range of the thrust amplitude of the engine, namely satisfying the following conditions:
T min representing a minimum thrust magnitude of the engine;
and setting the amplitude range of the thrust angle of the aircraft, namely meeting the following requirements: 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 not less than (R + H)/R 0 Wherein R represents the average radius of the moon, and H represents a preset safe height threshold value;
the objective function is: j ═ m 0 -m F (ii) a Wherein m is F Representing the aircraft mass at the terminal moment.
4. The lunar surface large-scale aircraft trajectory optimization method based on lossless convexity and Chebyshev orthogonal fitting method according to claim 3, wherein the constraint conditions after the relaxation variables are introduced in the second step and transformed are as follows:
in the formula, z 0 ＝ln m 0 ，z F ＝ln m F Z ═ ln m (τ), τ denotes normalized time; σ represents a relaxation variable; u. of r 、u θ 、u φ Respectively representing the radial, tangential and normal accelerations of the aircraft in a spherical coordinate system.
5. The lunar surface large-range aircraft trajectory optimization method based on the lossless convexity and Chebyshev orthogonal fitting method according to claim 4, wherein the dynamic model equation after the time control quantity is introduced and transformed in the second step is as follows:
in the formula u t Representing the terminal time control quantity, is the derivative of time t over normalized time τ, i.e.
Represents the derivative of z with respect to normalized time τ; in the ascending section, c ═ c 1 In the descending section, 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 quantity is expressed as:
p(τ)＝[u r ；u θ ；u φ ；σ；u t ]
the state quantity and the control quantity are uniformly defined as a variable matrix, and are expressed as follows:
y(τ)＝[x(τ)；p(τ)]
expressing the transformed kinetic model equation as the derivative of x (τ)
The function for y (τ) is expressed as:
at an initial time τ 0 The state quantity and the control quantity are subjected to first-order Taylor expansion, and the finally obtained transformed dynamic model equation is in the form as follows:
in the formula, y (τ) 0 )＝[x(τ 0 )；p(τ 0 )]And represents an initial variable matrix.
6. The lunar surface large-range aircraft trajectory optimization method based on the lossless convexity and the Chebyshev orthogonal fitting method as claimed in claim 5, wherein the discretization process of the transformed kinetic model equation by the Chebyshev orthogonal fitting method in the second step comprises time domain transformation and fitting point discretization; the discretized kinetic model equation is as follows:
in the formula, X (τ) i ) Representing the state quantity after dispersion; y (τ) k ) Representing a discrete variable matrix consisting of state quantity and control quantity; d k,i Representing derivative values of Lagrange basis functions at the discrete points after the dispersion; i and k each represent the position of an interpolation node.
7. The lunar surface wide-range aircraft trajectory optimization method based on the lossless convexity and Chebyshev orthogonal fitting method as claimed in claim 6, 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 true CN114872934A (en) | 2022-08-09 |
| CN114872934B 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 (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20100228409A1 (en) * | 2008-12-16 | 2010-09-09 | 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 |
| US20130187008A1 (en) * | 2012-01-25 | 2013-07-25 | Piyush Grover | System and Method for Controlling Motion of Spacecrafts |
| 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 |
-
2022
- 2022-05-24 CN CN202210568718.5A patent/CN114872934B/en active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20100228409A1 (en) * | 2008-12-16 | 2010-09-09 | 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 |
| US20130187008A1 (en) * | 2012-01-25 | 2013-07-25 | Piyush Grover | System and Method for Controlling Motion of Spacecrafts |
| 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 |
|---|
| 陈嘉澍;: "基于凸优化和切比雪夫伪谱法的再入轨迹优化", 电子技术与软件工程, no. 16, pages 71 - 75 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114872934B (en) | 2023-10-13 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Sun et al. | Adaptive backstepping control of spacecraft rendezvous and proximity operations with input saturation and full-state constraint | |
| CN103488814B (en) | Closed loop simulation system suitable for controlling attitude of reentry vehicle | |
| Zheng et al. | Planar path following control for stratospheric airship | |
| Huang et al. | Saturated backstepping control of underactuated spacecraft hovering for formation flights | |
| Wang et al. | A robust predictor–corrector entry guidance | |
| Boyi et al. | Performance limitations in trajectory tracking control for air-breathing hypersonic vehicles | |
| Wang et al. | Predictor-corrector entry guidance with waypoint and no-fly zone constraints | |
| Wang et al. | Optimal trajectory-tracking guidance for reusable launch vehicle based on adaptive dynamic programming | |
| Gong et al. | Mars entry guidance for mid-lift-to-drag ratio vehicle with control constraints | |
| Liu et al. | Mars entry trajectory planning with range discretization and successive convexification | |
| Zhang et al. | Aerocapture trajectory planning using hierarchical differential dynamic programming | |
| Chen et al. | Steady Glide Dynamics and Guidance of Hypersonic Vehicle | |
| CN114872934A (en) | Lunar surface large-range aircraft trajectory optimization method based on lossless convexity and Chebyshev orthogonal fitting method | |
| Long et al. | Vector trajectory method for obstacle avoidance constrained planetary landing trajectory optimization | |
| Liu et al. | Predictor-corrector guidance for entry with terminal altitude constraint | |
| De Oliveira et al. | Reusable launchers re-entry controlled dynamics simulator | |
| Wang et al. | Optimal 3-dimension trajectory-tracking guidance for reusable launch vehicle based on back-stepping adaptive dynamic programming | |
| He et al. | Surrogate-based entire trajectory optimization for full space mission from launch to reentry | |
| Li et al. | Re-entry guidance method based on decoupling control variables and waypoint | |
| Wang et al. | Event-triggered trajectory-tracking guidance for reusable launch vehicle based on neural adaptive dynamic programming | |
| 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 | |
| Jiang | Robust optimization of Mars entry trajectory under uncertainty | |
| Kazmierczak et al. | Improvements to modeling and trajectory simulation of mars aerocapture and aerobraking |
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 |