CN106020231A - Hypersonic air vehicle reentry trajectory optimization method based on reentry point parameter - Google Patents

Abstract

The invention provides a hypersonic air vehicle reentry trajectory optimization method based on a reentry point parameter. A series of problems such as long optimization time of a reentry trajectory optimization process, separate design of a vertical trajectory and a horizontal trajectory, inability of guaranteeing global optimization, fast optimization requiring model simplification are solved. An accurate dynamical model considering earth flattening, convected acceleration, Coriolis acceleration items is established, and various complicated constraint conditions are analyzed, the hypersonic air vehicle reentry trajectory optimization method focuses on researching a main factor of a reentry trajectory optimization result-reentry point parameter uncertainty, and in addition, by analyzing the uncertainty expansion problem of the reentry point parameter, the mapping relation between the reentry optimization trajectory and the reentry point parameter is acquired, and therefore a reentry trajectory is optimized quickly according to the reentry point parameter. The hypersonic air vehicle reentry trajectory optimization method is suitable for being used in an air vehicle trajectory optimization field, and is advantageous in that calculation efficiency is high, and a strong engineering application value is provided.

Description

Hypersonic aircraft reentry trajectory optimization method based on reentry point parameter

Technical field

The present invention relates to track optimizing technical field, be specifically related to a kind of hypersonic flight based on reentry point parameter Device reentry trajectory optimization method.

Background technology

In recent years, hypersonic aircraft gradually becomes and implements whole world rapid strike and effective work of holding air superiority Tool, is paid close attention to by countries in the world.The height that the U.S. started by advanced person's research project office of Ministry of National Defence is leading in June, 2003 is ultrasonic Speed aircraft project (CAV), this aircraft project obtains preliminary progress, have developed Lockheed Martin Corporation and open (CAV-L rises for the lifting body profile (CAV-H, lift-drag ratio scope is about 3.5～5.0) sent out and the improvement bipyramid appearance of Boeing Resistance is about 2.0～2.5 than scope).Russia's lightning science produces association and takes the lead in starting to design iron hammer in 2012 (Hammer) hypersonic aircraft, European Space Agency, Japan, India and China the most all step up to carry out hypersonic aircraft afterwards Research work.But, hypersonic aircraft flight speed is very fast (flight Mach number is generally higher than 5), and flight environment of vehicle exists Bigger uncertainty, on the other hand the flight time of re-entry flight track is long, and shared thru-flight time scale is high, flight bar Part is severe, and end strike effect is had highly important impact.In order to ensure that aircraft efficient stable flies, optimize design one The flight path of flight path, especially reentry stage is particularly important.

The purpose of track optimizing is to meet Dynamic Constraints, boundary condition, process about to determine in flight course simultaneously The optimum control amount of the constraints such as bundle, essence is a kind of Optimal Control Problem.Track optimizing problem substantially can be divided into two classes: indirectly Method and direct method.Indirect method is by control variable being expressed as state variable and the function of association's state variable, then solving 2 points Boundary value problem.Direct method is by being solved along flight path discretization, the direct method using parameter optimization by the equation of motion, To target function direct searching optimization.Indirect method solving precision is higher, but it is generally required to accurately provide the initial value of association's state variable, and be somebody's turn to do The proposition of initial value is the most often difficult to ensure that its accuracy.Therefore direct method occupies the leading of track optimizing field gradually.Pseudo-spectrometry, makees For most widely used a kind of direct method, joining simultaneously discrete state variable and a control variable by some and solve, solving precision is relatively Height, applies more convenient.

But, above-mentioned optimization method is often the longest, it is impossible to realizes completing on line to calculate, needs off-line to complete with reference to rail The optimization of mark, often need to bind entrance aircraft before transmission.The method optimized on some lines also tends to need to be to handled environment Simplify, lack the research to three-dimensional optimized track.

The finite time track such as disclosing a kind of hypersonic aircraft in CN201410216389.3 quickly generates Method, the method, by track optimizing problem is converted into convex optimization problem, reaches the purpose of rapid solving.But, the party Method needs special CVXGEN software to be compiled optimization problem, and coding and conversion process are complex, and optimize process In do not consider the compression of the Earth, the aceleration of transportation and the impact of Coriolis acceleration item and no-fly zone and way point constraint, only test Demonstrate,prove the motion of fore-and-aft plane, sidestep maneuver track optimizing has been lacked necessary discussion.

CN201510051589.2 discloses a kind of hypersonic aircraft reentry trajectory based on goal programming online Optimization method, the method uses speed-elevation plane method to calculate reentry corridor, the angle of attack is set to piecewise linear function, from And realize reentry trajectory on-line optimization.But the method needs the change curve of the most prefabricated angle of attack, longitudinal track and transverse path Need separately design, use sequential quadratic programming algorithm or interior point method cannot ensure global optimum's performance of result.

Summary of the invention

It is an object of the invention to provide a kind of hypersonic aircraft reentry trajectory optimization side based on reentry point parameter Method, this invention solves and optimizes overlong time, longitudinal track and transverse path need during reentry trajectory optimizes in prior art Separately to design, cannot ensure that global optimum maybe must carry out the technical problem of model simplification ability rapid Optimum.

The present invention provides the basic ideas of method: being conceived to research affects the principal element of reentry trajectory optimum results Reentry point parameter uncertainty, by being analyzed the uncertain expansion problem of reentry point parameter, obtains reentering optimization track With the mapping relations of reentry point parameter, thus realize one reentry trajectory of rapid Optimum according to reentry point parameter.The method is suitable for Optimal design of trajectory in hypersonic aircraft reentry stage.

See Fig. 1, the invention provides a kind of hypersonic aircraft reentry trajectory optimization side based on reentry point parameter Method, comprises the following steps:

Step S100: consider the compression of the Earth, the aceleration of transportation and the impact of Coriolis acceleration item, set up hypersonic flight The kinetic model of device re-entry；

Step S200: Various Complex constraint is analyzed, sets up Non-linear Optimal Model；Described Complex Constraints can be Process constraints, end conswtraint, no-fly zone and way point constraint, be certainly not limited to this.

Step S300: use Gauss puppet spectrometry to solve the optimization problem under name re-entry mode parameter, obtain name optimization Track and corresponding state variable, be normalized discrete time point, obtains being available for the reference time benchmark of follow-up use；

Step S400: according to point collocation solve the polynomial principle of GENERALIZED CHAOTIC set up meet probability distribution reenter point-like State parameter joins point sampling space, calculates corresponding weight and orthogonal polynomial；

Step S500: each in sample space assembles a little for optimizing initial value in step S400, uses Gauss puppet spectrometry Solve, obtain series of optimum trajectory parameters, then solve the output corresponding with step S300 Plays discrete time point Variable；

Step S600: calculate the polynomial coefficient of GENERALIZED CHAOTIC；

Step S700: according to the polynomial coefficient of GENERALIZED CHAOTIC obtained in step S600, solve and fly for hypersonic The optimization track of row device true reentry point parameter；

(1) consider the compression of the Earth, the aceleration of transportation and the impact of Coriolis acceleration item, set up hypersonic aircraft and reenter The kinetic model of inflight phase is as follows:

$\begin{array}{c}\stackrel{\·}{r}=V\sin \mathrm{\θ}\\ \stackrel{\·}{\mathrm{\λ}}=\frac{V\cos \mathrm{\θ}\sin \mathrm{\σ}}{r\cos \mathrm{\φ}}\\ \stackrel{\·}{\mathrm{\φ}}=\frac{V\cos \mathrm{\θ}\cos \mathrm{\σ}}{r\cos \mathrm{\φ}}\\ \stackrel{\·}{V}=-D+g_r\sin \mathrm{\θ}+C_V+{\tilde{C}}_V\\ \stackrel{\·}{\mathrm{\θ}}=\frac{1}{V}\[L\cos \mathrm{\μ}+(\frac{V^2}{r}+g_r)\cos \mathrm{\θ}\]+C_{\mathrm{\θ}}+{\hat{C}}_{\mathrm{\θ}}+{\tilde{C}}_{\mathrm{\θ}}\\ \stackrel{\·}{\mathrm{\σ}}=\frac{L\sin \mathrm{\μ}}{V\cos \mathrm{\θ}}+\frac{V}{r}\cos \mathrm{\θ}\sin \mathrm{\σ}\tan \mathrm{\φ}+C_{\mathrm{\σ}}+{\hat{C}}_{\mathrm{\σ}}+{\tilde{C}}_{\mathrm{\σ}}\end{array}---\left(1\right)$

Wherein, r be the earth's core away from, λ be longitude, φ be dimension, V be speed, θ be flight path angle, σ be flight path yaw angle, α Being angle of heel for the angle of attack, μ, wherein r, λ, φ, V, θ and σ are state variable, α and μ is control variable, and D is that aerodynamic drag is accelerated Degree, L are lift acceleration, another g_rFor acceleration of gravity along the component in arrow direction, the earth's core, g_ωFor the component in arrow direction, vertical the earth's core, g_rAnd g_ωCan be expressed as:

$\begin{array}{c}{g}_r=-\frac{{\mathrm{\μ}}_E}{r^2}\[1+\frac{3}{2}J{\left(\frac{a_e}{r}\right)}^2(1-5{\sin }^2\mathrm{\φ})\]\\ g_{\mathrm{\ω}}=-\frac{3{\mathrm{\μ}}_E}{r^2}J{\left(\frac{a_e}{r}\right)}^2\sin \mathrm{\φ}\end{array}---\left(2\right)$

Wherein, μ_EIt is to consider the Section 2 zonal harmonic coefficient of the compression of the Earth, a for Gravitational coefficient of the Earth, J_eFor earth semi-major axis. Aerodynamic drag and lift can be expressed as:

$\begin{array}{c}D=\frac{C_D{\mathrm{\ρV}}^2{S}_{ref}}{2m}\\ L=\frac{C_L{\mathrm{\ρV}}^2{S}_{ref}}{2m}\end{array}---\left(3\right)$

Wherein m is quality, S_refIt is atmospheric density for aircraft area of reference, ρ, can be expressed asWherein ρ₀=1.225kg/m³、h₀=7100m, h are height.C_DAnd C_LIt is respectively resistance and lift coefficient.The earth in accurate model is certainly Turn item (C_V、C_θ、C_σ), Coriolis acceleration itemWith aceleration of transportation itemCan be expressed as:

$\begin{array}{c}{C}_V=g_{\mathrm{\ω}}\·(\cos \mathrm{\σ}\cos \mathrm{\θ}\cos \mathrm{\φ}+\sin \mathrm{\θ}\sin \mathrm{\φ})\\ C_{\mathrm{\θ}}=(g_{\mathrm{\ω}}/V)\·(-\cos \mathrm{\σ}\sin \mathrm{\θ}\cos \mathrm{\φ}+\cos \mathrm{\θ}\sin \mathrm{\φ})\\ C_{\mathrm{\σ}}=-(g_{\mathrm{\ω}}/V\cos \mathrm{\θ})\·\cos \mathrm{\φ}\sin \mathrm{\φ}\\ {\hat{C}}_{\mathrm{\θ}}=2{\mathrm{\ω}}_E\cos \mathrm{\φ}\sin \mathrm{\σ}\\ {\hat{C}}_{\mathrm{\σ}}=2{\mathrm{\ω}}_E\·(\sin \mathrm{\φ}-\tan \mathrm{\θ}\cos \mathrm{\φ}\sin \mathrm{\σ})\\ {\tilde{C}}_V={\mathrm{\ω}}_E^{2}r\cos \mathrm{\φ}\·(\sin \mathrm{\θ}\cos \mathrm{\φ}-\cos \mathrm{\θ}\sin \mathrm{\φ}\cos \mathrm{\σ})\\ {\tilde{C}}_{\mathrm{\θ}}=({\mathrm{\ω}}_E^{2}r/V)\·\cos \mathrm{\φ}\·(\cos \mathrm{\θ}\cos \mathrm{\φ}+\sin \mathrm{\θ}\sin \mathrm{\φ}\cos \mathrm{\σ})\\ {\tilde{C}}_{\mathrm{\σ}}=({\mathrm{\ω}}_E^{2}r/V\cos \mathrm{\θ})\·\cos \mathrm{\σ}\sin \mathrm{\φ}\cos \mathrm{\φ}\end{array}---\left(4\right)$

Wherein ω_EIt it is rotational-angular velocity of the earth.

In order to improve computational efficiency and computational accuracy, the present invention carries out nondimensionalization process when kinetics equation model. The earth's core is away from using earth radius R with height etc.₀=6378.14km carries out change of scale, speed and time and is respectively adopted Convert；

(2) Various Complex constraint (including process constraints, end conswtraint, no-fly zone and way point constraint etc.) is carried out point Analysis, sets up Non-linear Optimal Model

First, process constraints (stationary point hot-fluid, dynamic pressure, overload) is retrained:

$\begin{array}{c}\stackrel{\·}{Q}=c_{\stackrel{\·}{Q}}{\mathrm{\ρ}}^{0.5}{V}^{3.15}\≤{\stackrel{\·}{Q}}_{\max }\\ q=\frac{1}{2}{\mathrm{\ρV}}^2\≤q_{\max }\\ n=\frac{\sqrt{L^2+D^2}}{{\mathrm{mg}}_0}\≤n_{\max }\end{array}---\left(5\right)$

Wherein,For stationary point heat flow modulus,q_maxAnd n_maxIt is respectively the maximum allowable of stationary point hot-fluid, dynamic pressure and overload Value；

Secondly, controlling process variable needs to retrain, and meets the requirement of actual control system:

$\begin{array}{c}{\mathrm{\α}}_{\min }\≤\mathrm{\α}\≤{\mathrm{\α}}_{\max }\\ {\mathrm{\μ}}_{\min }\≤\mathrm{\μ}\≤{\mathrm{\μ}}_{\max }\end{array}---\left(6\right)$

It addition, the SOT state of termination the most also has certain constraint:

$\begin{array}{c}{V}_f\≥V_{f\min }\\ h_f\≥h_{f\min }\\ \left|{\mathrm{\θ}}_f\right|\≤{\mathrm{\θ}}_{f\max }\\ \left|{\mathrm{\φ}}_f\right|\≤{\mathrm{\φ}}_{f\max }\end{array}---\left(7\right)$

Wherein, subscript " f " represents the respective SOT state of termination；

Finally, no-fly zone is studied, the present invention illustrates as a example by cylinder no-fly zone and ellipsoid no-fly zone.Circle The statement of post no-fly zone needs two parameter (major semiaxisAnd semi-minor axis).Ellipsoid no-fly zone also needs to height parameter State.Aircraft can be defined as (φ to the distance at center, no-fly zone in launching coordinate system_d,λ_d,h_d), the most no-fly District's constraint can be expressed as

The present invention can also be for the reentry trajectory of additional way point constraint (position that track has to pass through) and is optimized, its The way point coordinate of middle horizontal path point constraint can be expressed as

After Complex Constraints is analyzed, binding kinetics equation, Non-linear Optimal Model can be expressed as

$\begin{array}{cc}find& \mathrm{\μ}\left(t\right)and\mathrm{\α}\left(t\right)\\ \min & J_m=-L(x_{0,m},\mathrm{\μ}(t),\mathrm{\α}(t\left)\right)or{J}_m=t_f\\ s.t.& \left\{\begin{array}{c}\stackrel{\·}{x}=f(x_{0,m},\mathrm{\μ}(t),\mathrm{\α}(t\left)\right)\\ g(x,\mathrm{\μ}(t),\mathrm{\α}(t),x_f)\≤0\end{array}\right.\end{array}---\left(9\right)$

Wherein, J_mFor optimization object function, it can be range or flight time.

(3) use Gauss puppet spectrometry solve name re-entry mode parameter under optimization problem, obtain name optimize track and Corresponding output variable Y=[r, λ, φ, V, θ, σ, μ, α, t_f,J_m], the discrete time point during optimization can be expressed asWherein N_tRepresent total number of discrete time point.Discrete time point can carry out change of scale to [0,1] district Between:Discrete time point after normalized can be expressed as

(4) solve the polynomial principle of GENERALIZED CHAOTIC according to point collocation and set up reentry point parameter uncertainty δ=[δ₁,δ₂, δ₃,δ₄,δ₅,δ₆] join point sampling space, the multi-dimensional orthogonal multinomial Ψ that jth output variable is corresponding_j(δ) can be by monotropic Measure polynomial tensor product to obtain:

${\mathrm{\Ψ}}_j\left(\mathrm{\δ}\right)=\underset{i=1}{\overset{6}{\Π}}{\mathrm{\ψ}}_i^{\left(l_i\right)}\left({\mathrm{\δ}}_i\right),j=1,2,...,N_Y---\left(10\right)$

WhereinRepresent l_iRank univariate polynomials, N_YRepresent the dimension of output variable.Many for each multi-dimensional orthogonal Item formula Ψ_j(δ), the combination of univariate polynomials exponent number has uniqueness.

(5) assemble a little for optimizing initial value with each in sample space, use Gauss puppet spectrometry to solve, can obtain Series of optimum track output variable Y^(m), and interpolation solves(wherein) corresponding Output variable.

(6) the polynomial coefficient of GENERALIZED CHAOTIC is calculated

${\hat{x}}_j=\underset{m=1}{\overset{Q}{\Σ}}Y\left({\mathrm{\δ}}_m\right){\mathrm{\Ψ}}_j\left({\mathrm{\δ}}_m\right){\mathrm{\τ}}_m,j=1,2,\mathrm{...},P---\left(11\right)$

Wherein Q is for joining a sum, δ_m=[δ_1,m,δ_2,m,δ_3,m,δ_4,m,δ_5,m,δ_6,m] represent that m assembles corresponding reentering Point parameter uncertainty, τ_mFor corresponding weight, Y (δ_m) for using m to assemble a conduct optimization initial value, count according to step S500 Output variable Y calculated^(m)。

(7) forecast by track or sensor obtains reentry point parameter accurately, thus the output variable of reentry trajectory Can be expressed as

$Y_i\left(\mathrm{\δ}\right)=\underset{j=1}{\overset{P}{\Σ}}{\hat{x}}_{ij}{\mathrm{\Ψ}}_j\left(\mathrm{\δ}\right),i=1,2,\mathrm{...},N_Y---\left(12\right)$

Wherein Y_i(δ) represent optimize track i-th output variable,For corresponding GENERALIZED CHAOTIC multinomial coefficient, P table Show the polynomial item number of GENERALIZED CHAOTIC.

During use, step S100～S600 can be completed by calculated off line, then by many for the GENERALIZED CHAOTIC of gained Binomial coefficient loads in the computer of hypersonic aircraft, and line completes step S700, and step S700 obtains really at aircraft Cutting and reenter after parameter according to the polynomial coefficient of GENERALIZED CHAOTIC, solve for hypersonic aircraft true reentry point parameter is excellent Change track, it is thus achieved that optimum results.

The technique effect of the present invention:

1, the hypersonic aircraft reentry trajectory optimization method based on reentry point parameter that the present invention provides, establishes and examines Considering the compression of the Earth, the aceleration of transportation and the precise kinetic model of Coriolis acceleration item, computational accuracy is higher, flies closer to actual Row situation.

2, the hypersonic aircraft reentry trajectory optimization method based on reentry point parameter that the present invention provides, by many Plant Complex Constraints (including process constraints, end conswtraint, no-fly zone and way point constraint etc.) to be analyzed, preferably plan optimization Path so that optimum results can meet Complex Battlefield Environments requirement.

3, the hypersonic aircraft reentry trajectory optimization method based on reentry point parameter that the present invention provides, by using Gauss puppet spectrometry, thus realize quick obtaining globally optimal solution.

4, the hypersonic aircraft reentry trajectory optimization method based on reentry point parameter that the present invention provides, by setting up Optimize track and the mapping relations of reentry point parameter, thus realize getting final product rapid Optimum only in accordance with a certain reentry point state determined Go out a reentry trajectory.The process employs off-line and the algorithm combined online, optimization process more for elapsed time is handed over Completed by off-line algorithm, saved the cost in line computation, improve computational efficiency.

Specifically refer to the hypersonic aircraft reentry trajectory optimization method based on reentry point parameter according to the present invention Various embodiments described below, by apparent for the above and other aspect making the present invention proposed.

Accompanying drawing explanation

The hypersonic aircraft reentry trajectory optimization method flow process based on reentry point parameter that Fig. 1 provides for the present invention is shown It is intended to；

Fig. 2 provides method to reenter the speed of initial value Optimization Solution example one gained according to 1 times of standard deviation for using the present invention The nominal speed altitude curve comparison diagram of degree altitude curve and Gauss puppet spectrometry calculated off line is (with maximum range as target Function)；

Fig. 3 provides method to reenter the speed of initial value Optimization Solution example one gained according to 2 times of standard deviations for using the present invention The nominal speed altitude curve comparison diagram of degree altitude curve and Gauss puppet spectrometry calculated off line is (with maximum range as target Function)；

Fig. 4 provides method to reenter the speed of initial value Optimization Solution example one gained according to 3 times of standard deviations for using the present invention The nominal speed altitude curve comparison diagram of degree altitude curve and Gauss puppet spectrometry calculated off line is (with maximum range as target Function)；

Fig. 5 provides method to reenter initial value Optimization Solution example two integration gained according to 3 times of standard deviations for using the present invention The nominal altitude time graph comparison diagram of high temporal curve and Gauss puppet spectrometry calculated off line is (with maximum range as target Function)；

Fig. 6 provides method to reenter initial value Optimization Solution example two integration gained according to 3 times of standard deviations for using the present invention The nominal latitude, longitude curve comparison figure of latitude, longitude curve and Gauss puppet spectrometry calculated off line is (with maximum range as target Function)；

Fig. 7 provides method to reenter initial value Optimization Solution example two integration gained according to 3 times of standard deviations for using the present invention The nominal speed altitude curve comparison diagram of speed-altitude curve and Gauss puppet spectrometry calculated off line is (with the shortest flight time For object function, increase way point constraint)；

Fig. 8 provides method to reenter initial value Optimization Solution example two integration gained according to 3 times of standard deviations for using the present invention Latitude, longitude curve with the nominal latitude, longitude curve comparison figure of Gauss puppet spectrometry calculated off line (with the shortest flight time is Object function, increases way point constraint).

Detailed description of the invention

The accompanying drawing of the part constituting the application is used for providing a further understanding of the present invention, and the present invention's is schematic real Execute example and illustrate for explaining the present invention, being not intended that inappropriate limitation of the present invention.

Two kind application example checking beneficial effects of the present invention are given below.

In each example, name initial value is both configured to h₀=80km, V₀=6500m/s, θ₀=0deg, λ₀=0deg, φ₀= 0deg、σ₀=90deg, chooses and affects 3 maximum initial values of reentry trajectory optimum results (highly, speed, flight path angle), false Fixed its meets Gauss distribution, and standard deviation is expressed asThen reality is again Enter dotted state can be expressed asWherein k_iExpression standard deviation multiple (| k_i| The probability of≤1 is 68.3%, | k_i| the probability of≤2 is 95.5%, | k_i| the probability of≤3 is 99.7%), according to 3 times of standard deviations Principle, the biggest, the actual probability occurred is the least.Cylinder no-fly zone it is, in this cylinder no-fly zone handled by each example Heart position isMajor semiaxis and semi-minor axis are respectivelyEllipsoid Center, no-fly zone isMajor semiaxis, semi-minor axis and height are respectivelyOther constrained parameters see table 1.

Table 1 constrained parameters table

Example one according to the output variable of reentry point state parameter calculation optimization track for guidance system, i.e. in step Output variable Y=[r, λ, φ, V, θ, σ, t in S700_f,J_m]。

The track optimizing result of 3 kinds of uncertain levels of this case study, | k_i|=1,2,3.In example, " reference value " represents Using the Gauss puppet spectrometry optimum results with actual reentry point state parameter as initial value, " optimal value " represents the base using the present invention In the optimization track that reentry point parameter uncertainty expansion method obtains with actual reentry point state parameter for initial value, " nominal value " Represent that using Gauss puppet spectrometry calculated off line in the name of reentry point state parameter is the optimum results of initial value.

The method that Fig. 2～4 compared for using the present invention to provide optimizes track with traditional Gauss puppet spectrometry calculated off line name The speed-altitude curve obtained.It can be seen that along with probabilistic increase, in the name of reentry point state parameter is initial value Deviation between optimum results and reference locus is increasing, and uses the optimal value of the present invention to be substantially unaffected.In table 2 right Compare the optimum results with maximum range as object function, it can be seen that the optimum results of the present invention (| k_i| when=2) and reference It is worth the most identical.Table 3 compared for the elapsed time of different calculation methods, it can be seen that it is excellent that the present invention provides needed for method The change time (| k_i| when=2) it is far smaller than Gauss puppet spectrometry, illustrate that the present invention provides method can meet aircraft computer Calculating requirement.

The table 2 optimum results contrast table with maximum range as object function

Level of uncertainty	Reference value	Optimal value (error)	Nominal value (error)
				|ki|=1	15073.86km	15093.63km (0.13%)	14630.36km (2.94%)
|ki|=2	15607.23km	15609.94km (0.02%)	14630.36km (6.26%)
				|ki|=3	16126.66km	16158.12km (0.20%)	14630.36km (9.28%)

Table 3 calculation consumption time contrast table (condition: notebook computer 2.0GB RAM, 3GHz CPU)

Example two according to the control variable of reentry point state parameter calculation optimization track for control system, i.e. in step Output variable Y=[θ, σ, t in S700_f,J_m]。

In example two study the worst situation of optimum results (uncertain level | k_i|=3), use the present invention to provide Method calculation optimization control variable, then adds in kinetics equation by optimal control variable, and integration obtains kinestate, and right State parameter and control parameter contrast.

Acquired results, as shown in Fig. 5～6, compared for height curve in time and the horizontal movement of Different Optimization method respectively Trail change rule.Can be seen that the optimization method based on reentry point parameter that the present invention uses, optimum results is close to reference value Effect is preferable, and name integration track has with reference value and separates largely, is not suitable for the control system of hypersonic aircraft System uses.

In order to further illustrate beneficial effects of the present invention, the shortest flight time is used to be simultaneously introduced boat as optimizing index Waypoint retrains.Acquired results, as shown in Fig. 7～8, compared for speed-altitude curve and the horizontal movement rail of integration gained respectively The Changing Pattern of mark.It can be seen that different optimizing index functions all can be used by the method for the present invention, optimal value and reference value Registration is higher, and optimum results is preferable.

In sum, the present invention expands as starting point with reentry point parameter uncertainty, has used generalized polynomial chaos The hybrid algorithm that method and traditional optimization combine, it is proposed that off-line prepares and the optimisation strategy of online planning, can efficient solution Certainly band Complex Constraints and the hypersonic aircraft reentry trajectory optimization problem of reentry point parameter uncertainty.The method is not only The calculating time is short, and implementation is simple, additionally it is possible to make optimization track meet Various Complex constraints, it is ensured that feasibility, simultaneously The method that the present invention proposes is without carrying out model simplification, and optimum results meets global optimum, has the strongest engineering application and is worth.

Although the present invention is disclosed above with preferred embodiment, so it is not limited to the present invention, and any this area is general Logical technical staff, without departing from the spirit and scope of the present invention, when making various change and retouching, the therefore protection of the present invention Scope is when defining in the range of standard depending on claims.

Those skilled in the art will understand that the scope of the present invention is not restricted to example discussed above, it is possible to carries out it Some changes and amendment, the scope of the present invention limited without deviating from appended claims.Although oneself is through in accompanying drawing and explanation Book illustrates and describes the present invention in detail, but such explanation and description are only explanations or schematic, and nonrestrictive. The present invention is not limited to the disclosed embodiments.

By to accompanying drawing, the research of specification and claims, when implementing the present invention, those skilled in the art are permissible Understand and realize the deformation of the disclosed embodiments.In detail in the claims, term " includes " being not excluded for other steps or element, And indefinite article " " or " a kind of " are not excluded for multiple.Some measure quoted in mutually different dependent claims The fact does not means that the combination of these measures can not be advantageously used.It is right that any reference marker in claims is not constituted The restriction of the scope of the present invention.

Claims (4)

1. a hypersonic aircraft reentry trajectory optimization method based on reentry point parameter, it is characterised in that include following Step:

Step S100: consider the compression of the Earth, the aceleration of transportation and the impact of Coriolis acceleration item, set up hypersonic aircraft again Enter the kinetic model of inflight phase；

Step S200: Various Complex constraint is analyzed, sets up Non-linear Optimal Model；

Step S300: use Gauss puppet spectrometry to solve the optimization problem under name re-entry mode parameter, obtain name and optimize track With corresponding output variable Y=[r, λ, φ, V, θ, σ, μ, α, t_f,J_m], discrete time point is normalized and is expressed asDiscrete time point, obtain reference time benchmark,

Described discrete time point is expressed asWherein N_tRepresent total number of discrete time point；

Described discrete time point is expressed as to [0,1] interval through change of scale

Step S400: solve the foundation of GENERALIZED CHAOTIC polynomial principle according to point collocation and meet the reentry point parameter of probability distribution not Definitiveness δ=[δ₁,δ₂,δ₃,δ₄,δ₅,δ₆] join point sampling space, calculate corresponding weight and orthogonal polynomial,

The multi-dimensional orthogonal multinomial Ψ that jth output variable is corresponding_j(δ) obtained by the tensor product of univariate polynomials:

${\mathrm{\Ψ}}_j\left(\mathrm{\δ}\right)=\underset{i=1}{\overset{6}{\Π}}{\mathrm{\ψ}}_i^{\left(l_i\right)}\left({\mathrm{\δ}}_i\right),j=1,2,...,N_Y---\left(10\right)$

WhereinRepresent l_iRank univariate polynomials, N_YRepresent the dimension of output variable；

Step S500: each in sample space assembles a little for optimizing initial value in step S400, uses Gauss puppet spectrometry to carry out Solve, obtain output variable Y that series of optimum track is corresponding^(m), and interpolation solves and step S300 Plays discrete time Point correspondenceThe output variable in moment；

Step S600: calculate the polynomial coefficient of GENERALIZED CHAOTIC by formula (11)

${\hat{x}}_j=\underset{m=1}{\overset{Q}{\Σ}}Y\left({\mathrm{\δ}}_m\right){\mathrm{\Ψ}}_j\left({\mathrm{\δ}}_m\right){\mathrm{\τ}}_m,j=1,2,...,P---\left(11\right)$

Wherein Q is for joining a sum, δ_m=[δ_1,m,δ_2,m,δ_3,m,δ_4,m,δ_5,m,δ_6,m] represent that m assembles a corresponding reentry point parameter Uncertain, τ_mFor corresponding weight, Y (δ_m) for using m to assemble a conduct optimization initial value；

Step S700: according to the polynomial coefficient of calculated GENERALIZED CHAOTIC in step S600, solve and fly for hypersonic The optimization track of row device true reentry point parameter.

Hypersonic aircraft reentry trajectory optimization method based on reentry point parameter the most according to claim 1, it is special Levying and be, forecast by track or sensor obtains real reentry point parameter in described step S700, described reentry point is joined The output variable optimizing track that number is corresponding is:

$Y_i\left(\mathrm{\δ}\right)=\underset{j=1}{\overset{P}{\Σ}}{\hat{x}}_{ij}{\mathrm{\Ψ}}_j\left(\mathrm{\δ}\right),i=1,2,...,N_Y---\left(12\right)$

Wherein, Y_i(δ) represent optimize track i-th output variable,Represent wide for corresponding GENERALIZED CHAOTIC multinomial coefficient, P The justice polynomial item number of chaos.

Hypersonic aircraft reentry trajectory optimization method based on reentry point parameter the most according to claim 1, it is special Levying and be, the kinetic model setting up hypersonic aircraft re-entry in described step S100 is:

$\begin{array}{c}\stackrel{\·}{r}=V\sin \mathrm{\θ}\\ \stackrel{\·}{\mathrm{\λ}}=\frac{V\cos \mathrm{\θ}\sin \mathrm{\σ}}{r\cos \mathrm{\φ}}\\ \stackrel{\·}{\mathrm{\φ}}=\frac{V\cos \mathrm{\θ}\cos \mathrm{\σ}}{r\cos \mathrm{\φ}}\\ \stackrel{\·}{V}=-D+g_r\sin \mathrm{\θ}+C_V+{\tilde{C}}_V\\ \stackrel{\·}{\mathrm{\θ}}=\frac{1}{V}\[L\cos \mathrm{\μ}+\left(\frac{V^2}{r}+g_r\right)\cos \mathrm{\θ}\]+C_{\mathrm{\θ}}+{\hat{C}}_{\mathrm{\θ}}+{\tilde{C}}_{\mathrm{\θ}}\\ \stackrel{\·}{\mathrm{\σ}}=\frac{L\sin \mathrm{\μ}}{V\cos \mathrm{\θ}}+\frac{V}{r}\cos \mathrm{\θ}\sin \mathrm{\σ}\tan \mathrm{\φ}+C_{\mathrm{\σ}}+{\hat{C}}_{\mathrm{\σ}}+{\tilde{C}}_{\mathrm{\σ}}\end{array}---\left(1\right)$

Wherein, r be the earth's core away from, λ be longitude, φ be dimension, V be speed, θ be flight path angle, σ be that flight path yaw angle, α are for attacking Angle, μ are angle of heel, and wherein r, λ, φ, V, θ and σ are state variable, α and μ is control variable, and D is that aerodynamic drag acceleration, L are Lift acceleration, another g_rFor acceleration of gravity along the component in arrow direction, the earth's core, g_ωFor the component in arrow direction, vertical the earth's core, g_rAnd g_ω It is expressed as:

$\begin{array}{c}{g}_r=-\frac{{\mathrm{\μ}}_E}{r^2}\[1+\frac{3}{2}J{\left(\frac{a_e}{r}\right)}^2\left(1-5{\sin }^2\mathrm{\φ}\right)\]\\ g_{\mathrm{\ω}}=-\frac{3{\mathrm{\μ}}_E}{r^2}J{\left(\frac{a_e}{r}\right)}^2\sin \mathrm{\φ}\end{array}---\left(2\right)$

Wherein, μ_EIt is to consider the Section 2 zonal harmonic coefficient of the compression of the Earth, a for Gravitational coefficient of the Earth, J_eFor earth semi-major axis；

Aerodynamic drag and lift are expressed as:

$\begin{array}{c}D=\frac{C_D{\mathrm{\ρV}}^2{S}_{ref}}{2m}\\ L=\frac{C_L{\mathrm{\ρV}}^2{S}_{ref}}{2m}\end{array}---\left(3\right)$

Wherein m is quality, S_refIt is atmospheric density for aircraft area of reference, ρ, is expressed asWherein ρ₀= 1.225kg/m³、h₀=7100m, h are height；

C_DAnd C_LIt is respectively resistance and lift coefficient, the earth rotation terms (C in accurate model_V、C_θ、C_σ), Coriolis acceleration itemWith aceleration of transportation itemIt is expressed as:

$\begin{array}{c}{C}_V=g_{\mathrm{\ω}}\·\left(\cos \mathrm{\σ}\cos \mathrm{\θ}\cos \mathrm{\φ}+\sin \mathrm{\θ}\sin \mathrm{\φ}\right)\\ C_{\mathrm{\θ}}=\left(g_{\mathrm{\ω}}/V\right)\·\left(-\cos \mathrm{\σ}\sin \mathrm{\θ}\cos \mathrm{\φ}+\cos \mathrm{\θ}\sin \mathrm{\φ}\right)\\ C_{\mathrm{\σ}}=-\left(g_{\mathrm{\ω}}/V\cos \mathrm{\θ}\right)\·\cos \mathrm{\φ}\sin \mathrm{\σ}\\ {\hat{C}}_{\mathrm{\θ}}=2{\mathrm{\ω}}_E\cos \mathrm{\φ}\sin \mathrm{\σ}\\ {\hat{C}}_{\mathrm{\σ}}=2{\mathrm{\ω}}_E\·\left(\sin \mathrm{\φ}-\tan \mathrm{\θ}\cos \mathrm{\φ}\sin \mathrm{\σ}\right)\\ {\tilde{C}}_V={\mathrm{\ω}}_E^{2}r\cos \mathrm{\φ}\·\left(\sin \mathrm{\θ}\cos \mathrm{\φ}-\cos \mathrm{\θ}\sin \mathrm{\φ}\cos \mathrm{\σ}\right)\\ {\tilde{C}}_{\mathrm{\θ}}=\left({\mathrm{\ω}}_E^{2}r/V\right)\·\cos \mathrm{\φ}\·\left(\cos \mathrm{\θ}\cos \mathrm{\φ}+\sin \mathrm{\θ}\sin \mathrm{\φ}\cos \mathrm{\σ}\right)\\ {\tilde{C}}_{\mathrm{\σ}}=\left({\mathrm{\ω}}_E^{2}r/V\cos \mathrm{\θ}\right)\·\cos \mathrm{\σ}\sin \mathrm{\φ}\cos \mathrm{\φ}\end{array}---\left(4\right)$

Wherein ω_EIt it is rotational-angular velocity of the earth.

Hypersonic aircraft reentry trajectory optimization method based on reentry point parameter the most according to claim 1, it is special Levy and be, described step S200 sets up Non-linear Optimal Model and comprises the following steps:

Step S210, retrains process constraints:

$\begin{array}{c}\stackrel{\·}{Q}=c_{\stackrel{\·}{Q}}{\mathrm{\ρ}}^{0.5}{V}^{3.15}\≤{\stackrel{\·}{Q}}_{\max }\\ q=\frac{1}{2}{\mathrm{\ρV}}^2\≤q_{\max }\\ n=\frac{\sqrt{L^2+D^2}}{{\mathrm{mg}}_0}\≤n_{\max }\end{array}---\left(5\right)$

WhereinFor stationary point heat flow modulus,q_maxAnd n_maxIt is respectively stationary point hot-fluid, dynamic pressure and the maximum permissible value of overload；

Step S220, retrains controlling process variable:

$\begin{array}{c}{\mathrm{\α}}_{\min }\≤\mathrm{\α}\≤{\mathrm{\α}}_{max}\\ {\mathrm{\μ}}_{\min }\≤\mathrm{\μ}\≤{\mathrm{\μ}}_{max}\end{array}---\left(6\right)$

Step S230, the SOT state of termination the most also has certain constraint:

$\begin{array}{c}{V}_f\≥V_{f\min }\\ h_f\≥h_{f\min }\\ \left|{\mathrm{\θ}}_f\right|\≤{\mathrm{\θ}}_{f\max }\\ \left|{\mathrm{\φ}}_f\right|\≤{\mathrm{\φ}}_{f\max }\end{array}---\left(7\right)$

Wherein, subscript " f " represents the SOT state of termination of corresponding state variable；

After Complex Constraints is analyzed, binding kinetics equation, Non-linear Optimal Model is expressed as:

find μ(t)and α(t)

min J_m=-L (x_0,m,μ(t),α(t))or J_m=t_f (9)

$\begin{array}{cc}s.t.& \left\{\begin{array}{c}\stackrel{\·}{x}=f\left(x_{0,m},\mathrm{\μ}\left(t\right),\mathrm{\α}\left(t\right)\right)\\ g\left(x,\mathrm{\μ}\left(t\right),\mathrm{\α}\left(t\right),x_f\right)\≤0\end{array}\right.\end{array}$

Wherein, J_mFor optimization object function, it is set to range or flight time.