CN109508030B - Collaborative analysis reentry guidance method considering multi-forbidden flight area constraint - Google Patents

Abstract

The invention discloses a collaborative analysis reentry guidance method considering multiple flight control zone constraints, which comprises the following steps: 1. then, entering guidance problem description; 2. a time-of-flight analytic solution based on the rotating earth model; 3. considering a multi-forbidden flight area and an analysis reentry guidance method of arrival time constraint; 4. a multi-aircraft arrival time coordination scheme. The invention has the advantages that: (1) the hypersonic glide trajectory flight time accurate analytic solution based on the rotating earth model is provided, the prediction error is kept within 3%, and the hypersonic glide trajectory flight time accurate analytic solution is suitable for the condition of a section with a very-high longitudinal lift-drag ratio. (2) The cooperative flight problem that a plurality of hypersonic glide aircrafts are guided to arrive at a target simultaneously in a multi-flight-forbidden region environment can be solved based on a three-dimensional reentry trajectory analytic solution and a flight time analytic solution. (3) A multi-objective numerical iterative programming scheme based on online ballistic simulation considering terminal time, speed and height requirements is designed. The quasi-Newton method is improved by using the directional derivative, and the trajectory simulation times are reduced.

Description

Collaborative analysis reentry guidance method considering multi-forbidden flight area constraint

Technical Field

The invention relates to a collaborative analysis reentry guidance method considering multi-forbidden flight zone constraint, and belongs to the fields of aerospace technology, weapon technology and guidance control.

Background

The cooperative formation flight technology has important significance for improving system combat efficiency, has wide application prospect, and has been widely and deeply researched in the field of unmanned aerial vehicles at present. However, for hypersonic gliding aircraft, formation flight is difficult. Because such aircraft are high in initial speed, but unpowered, it is necessary to precisely manage flight energy through appropriate lateral maneuvers. However, for aircrafts with different skid conditions and different interference degrees, the speed attenuation curve and the transverse maneuvering amplitude are different, so that the formation of a flying formation with relatively stable relative position is difficult.

In the aspect of cooperative formation flight technology, common formation forms include wedge formation, echelon formation, transverse formation, longitudinal formation, V-shaped formation and the like, and complex tasks such as cooperative reconnaissance, defense, attack and the like can be realized. In addition, by utilizing the leader to generate the updraft in the vortex, the formation flight can effectively save fuel and increase the range. However, under the constraint of time-varying characteristics of speed and transverse maneuvering amplitude, no theory or method capable of effectively dealing with the problem of cooperative formation flight of hypersonic glide aircrafts exists at present.

Yu w., Chen w., Analytical entry for no-fly-zone avoidance (analytic reentry guidance law for no-fly zone avoidance) [ J ] Aerospace Science and Technology, 2018, 72: 426, 442, an earth rotation compensation model with high precision is proposed, a three-dimensional reentry trajectory analytic solution with better performance is obtained, and then an incremental roll reversal sequence planning scheme is designed based on the analytic solution, so that the aim of coping with a large number of flight forbidden zones at a low roll reversal cost is fulfilled. In the following description of the present invention, it is referred to as "document one".

Gaxley C, Atmospheric reentry [ M ] The RAND Corporation, 1960. time of flight analytical solutions under constant longitudinal lift-drag ratio conditions are given for a fixed earth model. In the following description of the present invention, it is abbreviated as "document two".

Disclosure of Invention

The invention aims to solve the problem that the hypersonic glide vehicle collaborative formation flying cannot be effectively dealt with at present, provides a collaborative analysis reentry guidance method considering multi-no-fly zone constraint, and solves the collaborative flying problem that a plurality of hypersonic glide vehicles are guided to simultaneously reach a target in a multi-no-fly zone environment.

The main content of the invention is divided into three aspects: (1) a high-precision flight time analytic solution is deduced for the rotating earth model; (2) aiming at a single aircraft, an analytic reentry guidance method considering multiple no-fly zones and arrival time constraints is designed; (3) aiming at a plurality of aircrafts, a cooperative flight scheme for guiding a plurality of hypersonic gliding aircrafts to simultaneously reach a target in a multi-no-fly zone environment is designed. The longitudinal reference profile planning method capable of considering both the energy management requirement and the arrival time requirement is designed by utilizing a longitudinal analytic solution and a time analytic solution, and the problem can be effectively solved.

The invention discloses a collaborative analysis reentry guidance method considering multiple flight control zone constraints, which comprises the following steps:

step 1: reentry guidance problem description

Here, a three-dimensional reentrant particle dynamics model of the spherical rotating earth is used, where the longitude λ, latitude φ and altitude H are used to describe the aircraft position information, while the velocity V, ballistic inclination γ and heading ψ are used to describe the velocity information.

The extremely high flying speeds result in a very harsh thermal environment. Therefore, in order to ensure that each subsystem of the aircraft is working properly, the reentry flight trajectory needs to meet the heat flux density

And the incoming flow pressure q and the overload n are restrained. Considering flight control systems with limited capabilities, it is desirable to limit the range and rate of change of the angle of attack and roll. In addition to conventional constraints, a number of circular no-fly zone constraints are further considered herein.

The horizontal distance from the aircraft to the target in the reentry section is S_TAEMAnd then terminates. The desired terminal height at this time is H_TAEMTerminal velocity is V_TAEMTerminal heading error | Δ ψ_TAEMAnd terminal roll angle σ_TAEMAnd | satisfies the constraint condition. In addition, the method can be used for producing a composite materialHere, the desired terminal time is also designated as t_TAEM。

Step 2: time-of-flight analytic solution based on rotating earth model

To cope with the time-of-arrival constraints, a time-of-flight resolution is derived here. The energy is defined as follows

Wherein V is velocity, H is height, R_eIs the earth's mean radius and μ is the gravitational constant. The derivative of the energy with respect to time t is

Wherein g is the acceleration of gravity. In the context of the rotating earth, there are complex non-linear equations

Where γ is the ballistic inclination, D is the drag, m is the mass, ω is_eIs the angular velocity of the earth's rotation, phi is the latitude, and psi is the heading angle. Right item 3 in formula (4)

Item 4

Is an inertial force component caused by the rotation of the earth. The inertial force component in the direction of velocity is considered herein as an additional resistance, as follows

Defining an equivalent resistance of

Substituting equations (3) to (5) into equation (2) and taking the reciprocal thereof to obtain

Defining a longitudinal lift-to-drag ratio of

L therein₁L cos (σ) is the component of lift in the longitudinal plane, plan

Is composed of

Wherein a is₀，a₁，a₂Are coefficients of a quadratic polynomial.

Now make use of

The profile estimates the equivalent drag (note as steady glide) under steady glide conditions

). Time derivative of trajectory inclination as follows

Suppose that

The steady glide longitudinal lift can be estimated as follows

Wherein, Delta L₁Is the inertia force along the vertical directionIs regarded as an additional longitudinal lift, the formula is as follows

At the same time, Δ D can also be reduced to be

Then

Can be estimated by

Substituting equation (12) into equation (6), and substituting E for V can be obtained

Because H < R_eTherefore, H can be set to be the middle value H of glide height^*Further let R^*＝R_e+H^*. Although the inertial force is a small quantity, when the velocity approaches the first cosmic velocity, the inertial force may cause the denominator of the above equation to be zero, and thus singularity occurs. To avoid the above singularity, the inertial force is considered as a small quantity and a first order Taylor expansion is performed as follows

Wherein the simplification symbol h_z1、h_mIs defined as follows

h_m＝2E+μ/R^*(16)

According to Δ L₁And Δ D, h can be_z1Is divided into two parts h_PAF1And h_PAF2As follows

h_z1＝h_PAF1+h_PAF2(17)

Wherein the content of the first and second substances,

h_PAF1＝-2R^*Vω_ecos(φ)sin(ψ) (18)

due to h_PAF2Is much less than h_PAF1Thus, a linear function pair h can be roughly employed_PAF2Carry out an approximation as follows

h_PAF2(E)≈K_PAF2(1)V+k_PAF2(0)(20)

Wherein the coefficient k_PAF2(1)And k_PAF2(0)Determined by two endpoints, as follows

Wherein h is_PAF2(E) Calculated by substituting the current state into equation (19), and h_PAF2(E_TAEM) Is calculated by substituting the desired terminal state into equation (19).

Actually, the hypersonic gliding aircraft flies to the target around one great circle marked as the generalized equator, the aircraft transversely maneuvers around the generalized equator, and the azimuth angle of the generalized equator

Can be considered as some weighted average of the aircraft heading angle, so that the invention utilizes the azimuth of the generalized equator in the event that the subsequent heading angle curve is not yet known

Substitution of h_PAF2(E) And h_PAF2(E_TAEM) Heading angle in the formula. Up to this point, in the formula (14), the other parameters except the argument E are constant values, and the formula (14) may be integrated. Through a large amount of derivation and arrangement, an expression of a flight time analytic solution can be obtained as follows

Wherein the coefficient k_t(1)、k_t(2)、k_t(3)、k_t(4)、k_t(5)、k_t(6)、k_t(7)Is expressed as follows

And step 3: analytic re-entry guidance method considering multiple no-fly zones and arrival time constraints

Aiming at a single aircraft, a reentry guidance method considering multiple no-fly zones and arrival time constraints is designed. According to the ballistic characteristics, the reentry process is generally divided into three stages: a descent phase, a steady glide phase and a height adjustment phase.

S31: descending section

In order to avoid excessive heat flux caused by falling into dense atmospheric layers, the descent segment aircraft glides at the maximum available angle of attack, zero roll angle. When the lift force is enough to support the aircraft to glide smoothly, the attack angle is smoothly transited to the reference attack angle, and then the stable gliding stage is entered.

S32: steady glide phase

The steady glide phase is the longest, most important and most complex reentry flight phase, and the guidance scheme of the phase rapidly plans a reference trajectory meeting the constraint of a no-fly area and arrival time on line by using a three-dimensional reentry trajectory analytic solution and a flight time analytic solution. The guidance process in the balanced gliding stage is designed as follows:

s321, designing a reference attack angle section, and determining a corresponding reference lift-drag ratio section;

s322, reasonably designing parameterized longitudinal lift-drag ratio by considering energy management requirements and arrival time constraints

Profile while allowing determination of a parameterized transverse lift-to-drag ratio

A section;

s323, in order to compensate the effect of the earth rotation

And

on the basis of the profile, the equivalent vertical and horizontal lift-drag ratios (

And

) A section;

s324, solving longitudinal profile parameters by utilizing a longitudinal analytical solution in the literature I and the flight time analytical solution obtained in the step 2 according to energy management and arrival time requirements;

and S325, updating a roll reversal sequence once according to the no-fly zone and the terminal position constraint by using a longitudinal analysis solution and a transverse analysis solution in the first document every 60S.

S326, adjusting a reference roll angle, tracking an equivalent longitudinal lift-drag ratio section, and reversing according to the planned roll;

and S327, in order to inhibit ballistic oscillation, introducing ballistic damping feedback into a reference attack angle and a roll angle so as to obtain a guidance instruction.

S328, limiting the roll angle instruction according to process constraints;

and S329, repeating the process from the step S322 until the smooth gliding phase is finished.

S33: height adjustment stage

At the last reversal point (note that the energy of the aircraft corresponding to this node is E_BR(nTR)) Thereafter, the aircraft enters the altitude adjustment phase, but at some node before the last reversal point (denoted as the corresponding aircraft energy E)_ST) The aircraft starts the preparation for the altitude adjustment phase. The invention designs a multi-target numerical iterative planning scheme based on online trajectory simulation. In the scheme, a quasi-Newton method iteration algorithm is improved by utilizing the directional derivative, so that the trajectory simulation times are reduced, and the iteration convergence speed can be greatly improved. The process of the multi-target numerical value iterative programming scheme is as follows:

s331, predicting the terminal time and the terminal state by using trajectory simulation;

s332, judging whether the terminal speed or the terminal time meets the requirement; if yes, skipping to step S333, otherwise skipping to step S334;

s333, adjusting a longitudinal lift-drag ratio profile and a last inversion point by using a quasi-Newton iteration method according to the terminal speed and the terminal time requirement, and jumping to the step S331;

s334, judging whether the height of the terminal meets the requirement; if yes, jumping to step S336, otherwise jumping to step S335;

s335, correcting the reference attack angle parameter according to the terminal height error, and skipping to the step S331;

and S336, finishing the iterative planning.

When E ═ E_STIn the guidance method, the multi-target numerical value iterative programming scheme is utilized to plan a longitudinal lift-drag ratio, a residual flight distance and a reference profile of relative energy of flight time, and follow-up tracks are finely adjusted to meet the requirements of terminal time, speed and height. Since only a few simulations of the subsequent small segment of trajectory are needed, no significant computational burden is incurred. When E is_st＞E＞E_BR(nTR)When the guidance method tracks in E ═ E_STThe reference profile of the longitudinal lift-drag ratio, the residual flight distance and the flight time relative energy obtained by planning the multi-target numerical iteration scheme at the moment E_BR(nTR)And when the angle is more than E, determining the reference roll angle by using a proportional guidance law to eliminate the course error, and tracking the flight time and the residual flight distance reference profile by finely adjusting the attack angle to ensure that the requirements of terminal time and speed are met.

And 4, step 4: multi-aircraft arrival time coordination scheme

In general, the time for an aircraft to reach a destination is affected by two important factors: the time control during the launch and the reentry flight, the former for a wide range of time adjustments and the latter for interference avoidance and fine time adjustments. However, since the present invention does not consider the flight procedure of the boost section and the push-down section, the start time of the reentry section is used instead of the transmission time, and the end time of the reentry section is used instead of the arrival time. In practical situations, the three flight processes need to be comprehensively considered. Although uncertainty factors during the boosting process may cause deviations in the starting time and starting state of the reentry segment, these deviation factors can be effectively overcome during the reentry flight using the guidance method described above.

Common temporal coordination problems are: synchronous to the same target, asynchronous to the same target, synchronous to different targets, asynchronous to different targets, and combinations thereof. Since the coping strategies for these problems are similar, the following describes the coordination scheme by taking the problem of synchronous arrival at the same target as an example.

Suppose there is n_HGVAn aircraft, denoted as HGV_i，i＝1，2，…，n_HGV. After the launch point and the target point are given, the initial state of the reentry section of the aircraft can be determined according to the performance of the booster rocket and the selected active section guidance scheme. And then, determining corresponding longitudinal lift-drag ratio profile parameters by using a longitudinal analytical solution according to the energy management requirement. Furthermore, the off-line trajectory simulation may be used to predict the time of flight t for all aircraft_EF(HGVi)，i＝1，2，…，n_HGVSelecting a flight time t with the longest time_EF(max)And reserve a transmission preparation time t_preAnd boost time of flight t_boostThe desired time of arrival t may be determined_TAEM＝t_pre+t_boost+t_EF(max)And further use t_TAEM-t_EF(HGVi)The re-entry period start time for each aircraft may be determined. Then, according to the selected guidance method for the boost phase, the corresponding launch time can be determined.

The invention has the advantages that:

(1) the hypersonic glide trajectory flight time accurate analytic solution based on the rotating earth model is invented, the prediction error is kept within 3%, and the hypersonic glide trajectory flight time accurate analytic solution is suitable for the situation of a section with a very-high longitudinal lift-drag ratio. The traditional prediction method does not consider the influence of earth rotation, is only suitable for the condition of constant longitudinal lift-drag ratio, and has larger error;

(2) the method can rapidly plan the multi-constraint reference trajectory meeting the no-fly zone, the arrival time, the terminal energy and the like on line based on the three-dimensional reentry trajectory analytic solution and the flight time analytic solution, and solves the problem of cooperative flight for guiding a plurality of hypersonic glide aircrafts to arrive at a target simultaneously in the multi-no-fly zone environment.

(3) A multi-objective numerical iterative programming scheme based on online ballistic simulation considering terminal time, speed and height requirements is designed. Meanwhile, in order to reduce the calculation burden, a quasi-Newton method is improved by using the directional derivative, and the ballistic simulation times are reduced.

Drawings

FIG. 1 is a schematic flow chart of the present invention;

FIG. 2 is a reference lift-to-drag ratio profile;

FIG. 3 is a flow chart of a multi-objective numerical iterative planning scheme during the height adjustment phase;

FIG. 4 is a ground projection of an aircraft trajectory;

FIG. 5 is an altitude-velocity curve for an aircraft;

FIG. 6 is an angle of attack versus time curve for an aircraft;

FIG. 7 is a roll angle versus time curve for an aircraft;

FIG. 8 is a plot of remaining flight distance versus energy for an aircraft;

FIG. 9 is a time-of-flight versus energy curve for an aircraft.

FIGS. 10a, b are Monte Carlo simulated aircraft reentry trajectories;

FIGS. 11a, b are aircraft altitude-speed curves for Monte Carlo simulations;

FIGS. 12a, b are aircraft angle of attack curves for Monte Carlo simulations;

FIGS. 13a, b are aircraft roll angle curves from Monte Carlo simulations;

FIGS. 14a and b are the terminal speed and terminal heading error distributions of the Monte Carlo simulation of the aircraft;

FIGS. 15a, b are Monte Carlo simulated aircraft arrival time profiles.

Detailed Description

The present invention will be described in further detail with reference to the accompanying drawings and examples.

The main content of the invention is divided into three aspects: (1) a high-precision flight time analytic solution is deduced for the rotating earth model; (2) aiming at a single aircraft, an analytic reentry guidance method considering multiple no-fly zones and arrival time constraints is designed; (3) aiming at a plurality of aircrafts, a cooperative flight scheme for guiding a plurality of hypersonic gliding aircrafts to simultaneously reach a target in a multi-no-fly zone environment is designed. The work flow diagram of the invention is shown in fig. 1, wherein, when the analytic reentry guidance method is designed, the control terminal speed and the arrival time are both controlled by the longitudinal lift-drag ratio, so that the control is an under-actuated problem.

The whole process comprises the following steps:

step 1: reentry guidance problem description

Here, a three-dimensional reentrant particle dynamics model of the spherical rotating earth is used, where the longitude λ, latitude φ and altitude H are used to describe the aircraft position information, while the velocity V, ballistic inclination γ and heading ψ are used to describe the velocity information.

The extremely high flying speeds result in a very harsh thermal environment. Therefore, in order to ensure that each subsystem of the aircraft is working properly, the reentry flight trajectory needs to meet the heat flux density

And the incoming flow pressure q and the overload n are restrained. Considering flight control systems with limited capabilities, it is desirable to limit the range and rate of change of the angle of attack and roll. In addition to conventional constraints, a number of circular no-fly zone constraints are further considered herein.

The horizontal distance from the aircraft to the target in the reentry section is S_TAEMAnd then terminates. The desired terminal height at this time is H_TAEMTerminal velocity is V_TAEMTerminal heading error | Δ ψ_TAEMAnd terminal roll angle σ_TAEMAnd | satisfies the constraint condition. In addition, the desired terminal time is also designated as t_TAEM。

Step 2: time-of-flight analytic solution based on rotating earth model

To cope with the time-of-arrival constraints, a time-of-flight resolution is derived here. The energy is defined as follows

Wherein V is velocity, H is height, R_eIs the earth's mean radius and μ is the gravitational constant. The derivative of the energy with respect to time t is

Wherein g is the acceleration of gravity. In the context of the rotating earth, there are complex non-linear equations

Where γ is the ballistic inclination, D is the drag, m is the mass, ω is_eIs the angular velocity of the earth's rotation, phi is the latitude, and psi is the heading angle. Right item 3 in formula (4)

Item 4

Is an inertial force component caused by the rotation of the earth. The inertial force component in the direction of velocity is considered herein as an additional resistance, as follows

Defining an equivalent resistance of

By substituting equations (33) to (35) into equation (32) and taking the reciprocal thereof, the expression

Defining a longitudinal lift-to-drag ratio of

L therein₁L cos (σ) is the component of lift in the longitudinal plane, and the longitudinal lift-drag ratio is formulated as follows,

wherein a is₀，a₁，a₂Are coefficients of a quadratic polynomial.

Now make use of

The profile estimates the equivalent drag (note as steady glide) under steady glide conditions

). Time derivative of trajectory inclination as follows

Suppose γ ≈ 0 and

the smooth gliding longitudinal lift L can be estimated_1(SG)As follows

Wherein, Delta L₁Is the component of the inertia force along the vertical direction and is regarded as the additional longitudinal lift force, and the formula is as follows

At the same time, Δ D can also be reduced to be

Then the equivalent resistance in the steady glide state

Can be estimated by

Substituting equation (42) into equation (36) and replacing V with E can result in

Because H < R_eTherefore, H can be set to be the middle value H of glide height^*Further, let the average earth-center distance R^*＝R_e+H^*. Although the inertial force is a small quantity, when the velocity approaches the first cosmic velocity, the inertial force may cause the denominator of the above equation to be zero, and thus singularity occurs. To avoid the above singularity, the inertial force is considered as a small quantity and a first order Taylor expansion is performed as follows

Wherein the simplification symbol h_z1、h_mIs defined as follows

h_m＝2E+μ/R^*(46)

According to Δ L₁And Δ D, h can be_z1Is divided into two parts h_PAF1And h_PAF2As follows

h_z1＝h_PAF1+h_PAF2(47)

Wherein the content of the first and second substances,

h_PAF1＝-2R^*Vω_ecos(φ)sin(ψ) (48)

due to h_PAF2Is much less than h_PAF1Thus, a linear function pair h can be roughly employed_PAF2Carry out an approximation as follows

h_PAF2(E)≈k_PAF2(1)V+k_PAF2(0)(50)

Wherein the coefficient k_PAF2(1)And k_PAF2(0)Determined by two endpoints, as follows

Wherein h is_PAF2(E) Calculated by substituting the current state into equation (49), and h_PAF2(E_TAEM) The desired terminal state is calculated by substituting the equation (49).

Actually, the hypersonic gliding aircraft flies to the target around one great circle marked as the generalized equator, the aircraft transversely maneuvers around the generalized equator, and the azimuth angle of the generalized equator

Can be considered as some weighted average of the aircraft heading angle, so that the invention utilizes the azimuth of the generalized equator in the event that the subsequent heading angle curve is not yet known

Substitution of h_PAF2(E) And h_PAF2(E_TAEM) Heading angle in the formula. Up to this point, in the formula (44), the other parameters except the argument E are constant values, and the formula (44) can be integrated. The expression of the time-of-flight analytical solution derived is as follows

Wherein the coefficient k_t(1)、k_t(2)、k_t(3)、k_t(4)、k_t(5)、k_t(6)、k_t(7)Is expressed as follows

The analysis of the time-of-flight analytical solution of the present invention compared with the results of the conventional method and the numerical ballistic simulation is shown in the first embodiment.

And step 3: analytic re-entry guidance method considering multiple no-fly zones and arrival time constraints

Aiming at a single aircraft, a reentry guidance method considering multiple no-fly zones and arrival time constraints is designed. According to the ballistic characteristics, the reentry process is generally divided into three stages: a descent phase, a steady glide phase and a height adjustment phase. The initial height of the descending section is high, the atmosphere is thin, and the aircraft quickly falls to the height. In the stage of stable gliding, the atmospheric density is moderate, the lifting force is enough to balance the gravity, and the gliding height is gently reduced. Thanks to the high lift and high speed, the gliding distance of the aircraft at this stage can reach tens of thousands of kilometers. When the aircraft is close enough to the target, the altitude adjustment phase is entered. The aircraft then obtains the desired flight height by appropriately adjusting the angle of attack. Compared with the previous guidance method (document one), there are two main improvements: (1) according to the energy management requirement and the arrival time requirement, a planning method of the reference profile is improved by using a flight time analytic solution and a longitudinal analytic solution; (2) in order to accurately fine-tune the track in the height adjustment stage, a numerical iterative planning scheme and a corresponding tracking strategy are designed according to terminal time, speed and height constraints. The following steps 4-6 describe the guidance schemes of the three phases in detail.

And 4, step 4: descending section guidance method design

In order to avoid excessive heat flux caused by falling into dense atmospheric layers, the descent segment aircraft glides at the maximum available angle of attack, zero roll angle. When the lift force is enough to support the aircraft to glide smoothly, the attack angle is smoothly transited to the reference attack angle, and then the stable gliding stage is entered.

And 5: design of guidance method in stable gliding stage

The steady glide phase is the longest, most important and most complex reentry flight phase, and the guidance scheme of the phase rapidly plans a reference trajectory meeting the constraint of a no-fly area and arrival time on line by using a three-dimensional reentry trajectory analytic solution and a flight time analytic solution. The guidance process in the balanced gliding phase is as follows:

s51, designing a reference attack angle section and determining a corresponding reference lift-drag ratio section;

s52 reasonably designing parameterized longitudinal lift-drag ratio in consideration of energy management requirements and arrival time constraints

Profile while allowing determination of a parameterized transverse lift-to-drag ratio

A section;

s53 is for compensating the effect of the earth rotation

And

on the basis of the profile, the equivalent vertical and horizontal lift-drag ratios (

And

) A section;

s54, solving longitudinal profile parameters by using a longitudinal analytical solution in the document I and the flight time analytical solution obtained in the step 2 according to energy management and arrival time requirements;

s55 updates the roll reversal sequence every 60S according to the no-fly zone and the terminal position constraint using the longitudinal analysis solution and the lateral analysis solution in the document one.

S56, adjusting the reference roll angle to track the equivalent longitudinal lift-drag ratio section and reversing according to the planned roll;

s57 introduces ballistic damping feedback into the reference angle of attack and roll in order to suppress ballistic oscillations, resulting in guidance commands.

S58 limiting the roll angle command according to the process constraint;

s59 repeats the above process from step S52 until the steady glide phase ends.

(1) Reference angle of attack and reference lift-to-drag ratio

To maximize aircraft capacity, the following reference angle of attack α is used_bslSection plane

Wherein E is_α＝-5.55×10⁷J/kg is in the meanThe steady glide and the vicinity of the shift point of the altitude adjustment phase. E_TAEMIs the desired terminal energy of the reentry segment α₁In the height adjustment phase, the reference angle of attack is gradually reduced to angle of attack α in order to obtain the desired terminal height₂Here, α₂Is initially set to 6 and is fine tuned before entering the height adjustment phase, in order to maintain glide in high disturbance conditions, an artificial limit α is imposed_bslAngle of attack of not less than 5 degree and reference α degree_bslCorresponding reference lift-to-drag ratio

The cross section is shown in figure 2.

(2) Improved reference profile

It is not easy to control the terminal energy and time using only the longitudinal lift-to-drag ratio profile, due to the coupling and mutual constraints of flight energy management and arrival time control. After a great deal of research, the invention finds the following form of longitudinal lift-drag ratio profile

The effective adjustment of the arrival time can be realized on the premise of ensuring that the energy requirement of the terminal is met.

Where E represents the energy at the current time, x_ERepresenting the energy at any time. Function f_L1/D(1)(x_E)、f_L1/D(2)(x_E)、f_L1/D(3)(x_E) Is defined as follows

And

the longitudinal lift-drag ratio corresponding to the three energy nodes is also three parameters of the profile, wherein the first two parameters are used for managing energy and controlling flight time, and an analytic solution is used for solving hereinafter. To make the roll angle converge to zero finally, let the third parameter

Wherein the content of the first and second substances,

is the ideal lift-drag ratio, coefficient corresponding to the expected terminal state

For roughly compensating the effect of the pneumatic pull-out,

is the current actual lift-drag ratio obtained by using pneumatic identification technology, and

is the ideal lift-drag ratio corresponding to the current state. In the following, when not specifying a function value corresponding to an energy node, the writing of an argument is omitted for convenience, e.g.

Is abbreviated as

After the reference and longitudinal lift-to-drag ratio profiles are determined, a parameterized transverse lift-to-drag ratio profile may be determined from geometric relationships

As follows

Wherein the content of the first and second substances,n_Ris the number of reversals that have occurred, sgn is a sign function, used to determine the initial roll direction of the aircraft,

is that

Is estimated by the following equation

In order to compensate for the effects of earth rotation, inertial forces are combined with aerodynamic forces to form equivalent aerodynamic forces. By reasonably analyzing the change rule of the inertia force

And

on the basis of the profile, the equivalent longitudinal lift-drag ratio profile in the form of the combination of the following inverse proportional functions can be established

Transverse lift-to-drag ratio profile

Here, the first and second liquid crystal display panels are,

wherein h is_mAs shown in the formula (46),

and

is an approximate formula of some complex inertia force combination terms, as follows

Coefficient k of approximation formula_h1(0)、k_h1(1)、k_h2(0)、k_h2(1)、k_h3(0)、k_h3(1)、k_h3(2)、k_h4(0)、k_h4(1)The calculation method is determined by the current state and the expected terminal state as follows.

1) Calculating the coefficient k_h1(0)And k_h1(1)

Wherein

Is the current longitudinal lift-drag ratio obtained in the previous guidance period, is the desired terminal lift-drag ratio, phi_TAEMIs the latitude of the target and the latitude of the target,

is the generalized equatorial azimuth heading angle of the aircraft,

is the generalized equatorial azimuth heading angle of the aircraftThe value at the end point.

2) Calculating the coefficient k_h2(0)And k_h2(1)

Wherein the content of the first and second substances,

3) calculating the coefficient k_h3(0)、k_h3(1)And k_h3(2)

4) Calculating the coefficient k_h4(0)And k_h4(1)

Wherein

(3) Longitudinal profile parameter calculation

Here, the invention is based on energy management and time of arrivalDetermining the parameters in equation (62) using a longitudinal and time-of-flight analytic solution

And

let A be a general form of a longitudinal lift-to-drag ratio polynomial, i.e., an array of coefficients for equation (37), as follows

Then, using the formula of the longitudinal analytic solution in the literature, the longitudinal range of the aircraft from the energy node E0 to the energy node E under the control of the longitudinal lift-drag ratio can be obtained as

Wherein the function f_xD(1)(E，E₀，A)，f_xD(2)(E，E₀A) and f_xD(3)(E，E₀The definition of A) is as follows

Note A₁、A₂、A₃And A₄Sets of polynomial coefficients associated with equations (62) - (63), respectively, are as follows

Then for the longitudinal lift-drag ratio profile shown in equation (62), the terminal longitudinal extent can be predicted to be

Meanwhile, for the vertical lift-drag ratio profile shown in the formula (62), the terminal time can be predicted as follows using the formula (53)

The parameters in the formula (62)

And

can be obtained by solving the following system of linear equations

Wherein s is_goIs the remaining flight distance, t_TAEMThe expected time of arrival.

However, on one hand, the solution of the equation set (90) is limited by the restriction of energy management requirements, on the other hand, the influence of lateral maneuver cannot be fully considered, and the solution of the equation set causes the roll angle to greatly oscillate, so that the forbidden flight zone is prevented from being avoided, and the practical application is not facilitated. The invention therefore proposes to overcome the time offset caused by the transverse manoeuvre by greatly adjusting the longitudinal lift-to-drag ratio profile. Through exploration, the profile parameters of the longitudinal lift-drag ratio are corrected to ensure that

Wherein the content of the first and second substances,

and

is the longitudinal profile parameter corresponding to the nominal trajectory in the ideal case,

and

it is the parameter that needs to be corrected according to the actual situation. Substituting the formula (91) into (88) and (89), and arranging to obtain

Wherein x is_D(bsl)(E_TAEME) and t (E)_TAEME) is the portion corresponding to the nominal trajectory, and Δ x_D(E_TAEME) and Δ t (E)_TAEMAnd E) is the correction of the actual situation, and the formula is as follows

In the nominal case of the standard,

and

the longitudinal extent of the aircraft should be such as to satisfy terminal constraints, i.e.

x_D(bsl)(E_TAEM，E)＝s_go-S_TAEM(97)

Then order

Substituting equation (93) into equation (97) can solve the problem

Note that the expression in the denominator of the above equation takes advantage of the following integral property

On the other hand, in the case of a liquid,

and

for eliminating terminal time error delta t_fBut not cause additional voyage. Using equations (94) and (96) one can obtain

Get it solved

Wherein, the notational symbol D₀Is expressed as follows

To predict the terminal time error Δ t_fIt is necessary to plan the nominal trajectory offline in advance and obtain two reference profiles, namely the range profile s_go(ref)(E) And time profile t_(ref)(E) Wherein the longitudinal lift-drag ratio profile of the nominal trajectory is determined only by energy management requirements, and the terminal time is determined by the time of launchAnd (4) performing inter-control. But the actual profile may be s_go(E) And t (E). Thus predicting the terminal time error Δ t_fIs Δ t_f(1)＝t(E)-t_(ref)(E) The second component is the range error Δ s_go＝s_go(E)-s_go(ref)(E) Induced Δ t_f(2). Now it is deduced that under reference control Δ s_goInduced time error Δ t_f(2). Substituting equation (98) into equation (95) yields

From Δ s by using the differential formula_goInduced time error Δ t_f(2)

Thus, the longitudinal section parameter calculation work is completed.

(4) Roll reversal sequence planning

Since the hypersonic gliding aircraft tracks the longitudinal reference profile mainly by adjusting the roll angle module value, after the longitudinal lift-drag ratio profile is determined, the roll angle module value profile is substantially determined. Thus, the aircraft can only control lateral maneuvers by timely changing the roll direction in order to avoid the no-fly zone and eventually reach the target. The invention adopts an analytic iteration scheme to plan a roll reversal sequence. Because the prior reversal points are preferentially adjusted to correspond to the no-fly zone constraint, the strategy needs fewer reversal times, can fully exert the transverse maneuverability of the aircraft and can correspond to a large amount of no-fly zone constraints. The policy flow is described below.

1) For interference resistance, the radius of all the no-fly zones is artificially enlarged, so that a safe distance is reserved between the aircraft and the boundary of the real no-fly zone. Considering the cumulative effect of ballistic errors, the radius increment is set as a linear function of the distance from the aircraft to the no-fly zone and gradually decreases as the aircraft approaches the target;

2) planning two reversal points by using a traverse analysis solution only according to the terminal position constraint;

3) starting from the nearest no-fly zone, searching a route point nearest to the no-fly zone by using an analytic solution to judge whether the constraint of the no-fly zone is met. If yes, continuing to judge the next no-fly zone, otherwise, executing 4);

4) and judging whether the no-fly zone is positioned before the last inversion point, between the last two inversion points or after the last inversion point. If it is before the penultimate inversion point, then execute 5); if it is between the last two inversion points, then 6) is performed; if it is after the last inversion point, then execute 7);

5) if the no-fly zone is located before the penultimate inversion point, the no-fly zone is far from the terminal point, so that the terminal position constraint can be not considered. The reverse point is adjusted by using the no-fly zone avoiding strategy to deal with the current no-fly zone constraint, and then the subsequent reverse point is adjusted by using the terminal transverse control strategy 1 to ensure that the terminal position requirement is met. Jump execution 8). A description will be given below of a terminal cross-flow control strategy 2 in a no-fly zone avoidance strategy, a terminal cross-flow control strategy 1, and 6);

6) if the no-fly zone is located between the last two reversal points, the planning of the roll reversal sequence needs to take the constraints of the no-fly zone and the requirements of the terminal position into consideration at the moment. Here, the no-fly zone avoidance strategy is also used for dealing with the no-fly zone constraint, but since the distance to the target is closer, the terminal traverse control strategy 2 is used for planning the subsequent inversion point so as to meet the terminal traverse requirement. Jump execution 8);

7) if the no-fly zone is located after the last inversion point, the no-fly zone constraint is ignored and the original inversion point sequence is retained because the no-fly zone is too close to the end point. Jump execution 8);

8) determining the next no-fly zone, returning to 3) and continuously judging whether the no-fly zone constraint is met or not until all the no-fly zone constraints are detected.

The following explains the principles of the no-fly zone avoidance maneuver, the terminal cross-flow control maneuver 1, and the terminal cross-flow control maneuver 2.

1) Flight-forbidden zone avoidance strategy

The strategy can be used for dealing with the constraint of the no-fly zone in two cases: the current no-fly zone is positioned before the last second reversal point; the current no-fly zone constraint lies between the last two inversion points.

In the first case, the terminal location constraint need not be considered in handling the current no-fly zone constraint, since it is far from the target. The method utilizes the analytic solution to preferentially adjust the existing reversal points so that the aircraft bypasses the no-fly zone from one side close to the no-fly zone. If the result is unsuccessful, further attempt is made to bypass the no-fly zone from the other side, or a new inversion point is added to deal with the no-fly zone constraint. After successfully avoiding the current no-fly zone, it is necessary to check whether the previous no-fly zone constraint is still satisfied, since the inversion sequence adjustment would change the previous trajectory. If not, further judgment is made as to which threat level is higher to decide whether to adopt the new inversion point sequence.

In the second case, due to the close distance to the target, the reverse point sequence planning needs to consider both the no-fly zone avoidance and the terminal position constraint. The invention preferentially adjusts the existing reversal points to correspond to the restraint of the no-fly zone, but the difference is that the aircraft bypasses the no-fly zone along the side of the no-fly zone close to the target. The subsequent process is the same as above.

2) Terminal cross-flow control strategy 1

This strategy is primarily directed to the case where the current no-fly zone is located before the penultimate inversion point. After the current no-fly zone constraint is responded by using the no-fly zone avoidance strategy result, the subsequent inversion point sequence does not meet the terminal position constraint due to the track adjustment. Therefore, the terminal traversal control strategy 1 will abandon the subsequent inversion points and re-plan the last two inversions according to the terminal position constraints.

3) Terminal cross-flow control strategy 2

This strategy is mainly directed to the case where the current no-fly zone constraint is located between the last two inversion points. Likewise, subsequent reversal points are abandoned after the no-fly zone avoidance maneuver is successfully invoked. Since the target is close, the last reversal point is adjusted to meet the terminal position constraint, and then whether the reversal point meets the shift condition of the steady glide phase and the altitude adjustment phase is judged. If not, then 1-2 reversals are added appropriately by analyzing the lateral maneuver rules.

(5) Reference roll angle in steady glide phase

The reference roll angle can be determined using the geometric relationship between the longitudinal lift-to-drag ratio and the total lift-to-drag ratio. Note n_REnergy corresponding to the secondary roll reversal is E_BR(nR). When E is_BR(nR)+ΔE≥E≥E_BR(nR+1)At + Δ E, the reference roll angle σ_bslIs composed of

Where the compensation value deltae causes the roll reversal to be properly advanced to compensate for the response lag caused by the roll rate limit.

(6) Commanded angle of attack and roll angle for steady glide phase

To suppress ballistic oscillations, a three-dimensional ballistic damping control technique is introduced, as follows

α_cmd＝α_bsl+cos(σ_bsl)k_γ(γ_SG-γ) (106)

Wherein, α_cmdAnd σ_cmdIs an order angle of attack and roll angle, k_γIs the feedback gain factor, gamma_SGIs the angle of inclination of the smooth gliding trajectory, which can be determined from the current state.

In order to meet the constraints of heat flux density, incoming flow pressure, overload and the like, the process constraint can be converted into the tilt angle constraint, and the size of the tilt angle instruction is further limited.

Step 6: design of guidance method in height adjustment phase

At the last reversal point (note that the energy of the aircraft corresponding to this node is E_BR(nTR)) Thereafter, the aircraft enters the altitude adjustment phase, but before the last reversal pointIs given (as corresponding aircraft energy E)_ST) The aircraft starts the preparation for the altitude adjustment phase. The invention designs a multi-target numerical iterative planning scheme based on online trajectory simulation. When E ═ E_STAnd then, obtaining a longitudinal lift-drag ratio, a residual flight distance and a reference profile of relative energy of flight time by using a multi-target numerical value iterative programming scheme, and accurately fine-tuning a subsequent track so as to meet the requirements of terminal time, speed and height. Since only a few simulations of the subsequent small segment of trajectory are needed, no significant computational burden is incurred. When E is_st＞E＞E_BR(nTR)When the guidance method tracks in E ═ E_STThe longitudinal lift-drag ratio, the remaining flight distance, and the reference profile of the time-of-flight versus-energy obtained from the multi-objective numerical iterative planning scheme, and when E_BR(nTR)And when the angle is more than E, determining the reference roll angle by using a proportional guidance law to eliminate the course error, and tracking the flight time and the residual flight distance reference profile by finely adjusting the attack angle to ensure that the requirements of terminal time and speed are met. The specific guidance scheme of the height adjustment phase is as follows:

(1) multi-target numerical iterative programming scheme

The multi-objective numerical iterative programming scheme takes into account terminal time, speed and altitude constraints, wherein the terminal time and speed are defined by parameters

And E_BR(nTR)Control, terminal height by reference angle of attack parameter α₂And (5) controlling. As shown in fig. 3, the multi-objective numerical iterative planning scheme flow is as follows:

s611, predicting the terminal time and the terminal state by using trajectory simulation;

s612, judging whether the terminal speed or the terminal time meets the requirement; if yes, skipping to step S613, otherwise, skipping to step S614;

s613, adjusting the longitudinal lift-drag ratio profile and the last inversion point by using a quasi-Newton iteration method according to the terminal speed and the terminal time requirement, and jumping to the step S611;

s614, judging whether the height of the terminal meets the requirement; if yes, jumping to step S616, otherwise jumping to step S615;

s615, correcting the reference attack angle parameter according to the height error of the terminal, and skipping to the step S611;

and S616, finishing iterative planning.

In step S613, the coping strategies of the present invention are as follows:

1) if the terminal speed and the terminal time are found to be greater than expected values through ballistic simulation prediction, the terminal speed and the terminal time can be increased appropriately

Reduce

And E_BR(nTR)；

2) If the terminal speed is greater than the desired value and the terminal time is less than the desired value, it may be reduced appropriately

And E_BR(nTR)；

3) If the terminal speed is less than the desired value and the terminal time is greater than the desired value, it may be increased appropriately

And E_BR(nTR)；

4) If both terminal speed and terminal time are less than desired, then the reduction may be appropriate

Increase of

And E_BR(nTR)。

Due to the fact that

And E_BR(nTR)Change in synchronism with the change of

As an auxiliary adjusting means, and designing the following relational expression

Wherein the content of the first and second substances,

and Δ E_BR(nTR)Respectively represent

And E_BR(nTR)Amount of change, k_f＝2×10^-7To adjust the coefficients. The above formula represents if E_BR(nTR)Change 0.1 × 10⁷kJ/kg, then

The corresponding change was 0.2. Thus, one search dimension is reduced, and the complexity of the algorithm is reduced.

Following Newton iteration method search for satisfying terminal time, speed and altitude constraints

And E_BR(nTR). Defining a vector x and a vector function G

Wherein the content of the first and second substances,

and

are respectively a parameter

And E_BR(nTR)And the terminal speed and time obtained by corresponding online trajectory simulation. Solving equation G (x) equal to 0Newton's iterative formula is

Where the superscript "(k)" represents the kth iteration, the matrix J_GIs the Jacobian matrix of the vector function G, as shown in equation (81), and

is J_GThe inverse matrix of (c).

Due to the fact that

And

needs to be obtained by integration, J_GThe solution cannot be resolved. If the conventional differential method is used for solving, 3 ballistic integrations need to be performed at each step, which causes a large computational burden.

In order to reduce the computational burden, the invention provides a J based on a directional derivative_GAnd solving the scheme. Definition vector p^(k)＝x^(k)-x^(k-1)Defining a vector q^(k)Is perpendicular to p^(k)A small vector with a mode length of 0.01. There is the following approximate calculation formula for the directional derivative

Wherein, | | p^(k)| and | q^(k)Respectively representing the vector p^(k)And q is^(k)Die length of (2). The relationship between the directional derivative and the partial derivative can be used

Wherein, theta_x(p^(k)) And theta_y(p^(k)) Is a vector p^(k)Angle of the x-axis and y-axis, respectively, and theta_x(q^(k)) And theta_y(q^(k)) Then is the vector q^(k)Respectively forms an included angle with the x axis and the y axis. Can be solved by the formula (113)

In the current iteration, only G (x) needs to be calculated^(k)+q^(k)) And G (x)^(k)) Here, only 2 ballistic simulations have to be performed, thereby greatly reducing the amount of computation. Since at the time of each iteration of a step,

amount of change of

Is in the order of about 0.1 to 1, and Δ E_BR(nTR)Of the order of 10⁶The above equations are usually pathological, but the problem can be solved by using an element exchange method.

In step S615, the reference angle of attack profile parameter α is fine-tuned using the following equation (85)₂Because of the higher accuracy of the formula, only one correction is needed here α₂The required change Δ α₂Is composed of

Wherein, C_Lf(est)Is a terminal lift coefficient obtained by a pneumatic identification technology,

is the estimated slope of the lift line, S_refIs the aircraft reference area, E_f、q_f、V_f、H_f、φ_f、ψ_fRespectively the terminal energy, dynamic pressure, speed, height, etc. obtained after the on-line trajectory simulation,Latitude, value of course angle, q_TAEM、V_TAEMRespectively desired terminal dynamic pressure and velocity values.

(2) Commanded angle of attack and roll angle at altitude adjustment stage

When E is_st＞E＞E_BR(nTR)The aircraft is in a preparation phase before the altitude adjustment phase. At the moment, the longitudinal lift-drag ratio section, the residual flight distance section and the flight time section which are obtained by the multi-target numerical value iterative programming scheme are tracked by the guidance scheme and are respectively recorded as

s_go(ref)(E)、t_ref(E) In that respect Reference roll angle sigma_bslIs calculated as follows

Wherein the content of the first and second substances,

here longitudinal lift-to-drag ratio

Mainly composed of

To ensure energy management and terminal time requirements in the presence of interference, a residual flight distance profile is designed and tracked that takes into account time error correction, as follows

s_go(ref)(2)(E)＝s_go(ref)(E)-k_sgoV(t-t_ref(E)) (117)

Wherein k is_sgo＝0.2s_go/s_go(ST)。s_go(ST)The residual flight distance corresponding to the moment when the multi-target numerical value iterative planning scheme is executed.

Thereafter, by substituting equation (116) into equation (107) in combination with equation (106), the angle of attack and roll can be commanded.

When E < E_BR(nTR)When the aircraft enters the height adjustment stage, the reference attack angle is still as shown in a formula (61), and the reference roll angle is determined by a proportional guidance law. The rate of change of the azimuth of the line of sight is

Where Δ ψ is the heading error. The proportional steering law generates the required lateral maneuvering acceleration of

Wherein, Delta L₂Is the inertia force edge L₂Component of direction, second term on right of above equation- Δ L₂The term/m is used to compensate for the effects of earth rotation. To prevent the initial roll angle from saturating, the effective guide ratio k is made_PNGradually changing from 2 to 4 with the remaining flying distance. On the other hand, in the case of smooth gliding, the lift acceleration a in the longitudinal plane_L1Counterbalanced with gravitational, centrifugal and inertial forces, i.e.

The reference roll angle is

Further, the design attack angle and roll angle commands are as follows

α_cmd＝α_bsl+k_α[s_go-s_go(ref)(2)(E)](122)

And 7: multi-aircraft arrival time coordination scheme design

Steps 3-6 complete the analytic guidance method design for a single aircraft, considering multiple no-fly zones and terminal time constraints, where arrival time coordination schemes are further studied for multiple aircraft. In general, the time for an aircraft to reach a destination is affected by two important factors: the time control during the launch and the reentry flight, the former for a wide range of time adjustments and the latter for interference avoidance and fine time adjustments. However, since the present invention does not consider the flight process of the boost section and the push-down section, the present invention replaces the transmission time with the start time of the reentry section and replaces the arrival time with the end time of the reentry section. Note that: in practical situations, the three sections of flight processes need to be comprehensively considered. Although uncertainty factors during the boosting process may cause deviations in the starting time and starting state of the re-entry phase, these deviation factors can be effectively overcome during the re-entry flight by using the guidance method described above.

Common temporal coordination problems are: synchronous to the same target, asynchronous to the same target, synchronous to different targets, asynchronous to different targets, and combinations thereof. Since the coping strategies for these problems are similar, the following describes the coordination scheme by taking the problem of synchronous arrival at the same target as an example.

Suppose there is n_HGVAn aircraft, respectively denoted as HGV_i，i＝1，2，…，n_HGV. After the launch point and the target point are given, the initial state of the reentry section of the aircraft can be determined according to the performance of the booster rocket and the selected active section guidance scheme. Then, the longitudinal lift-drag ratio profile parameter is expressed by the same formula (98)

The corresponding longitudinal lift-drag ratio profile parameters are determined using a longitudinal analytical solution based only on energy management requirements, as shown in equation (94). Here, the subscript "hgv (i)" refers to the ith aircraft.

Further, the time of flight t for each aircraft may be predicted using an offline ballistic simulation_EF(HGVi)，i＝1，2，…，n_HGVThe expected time of arrival t can be determined by selecting the longest one of the times and reserving the transmit preparation and boost flight times_TAEMAnd further use t_TAEM-t_EF(HGVi)The re-entry period start time for each aircraft may be determined. Then, according to the selected guidance method for the boost phase, the corresponding launch time can be determined.

Example (b):

example one

In this embodiment, the flight time analytic solution of the present invention is compared with the conventional method and the numerical value trajectory simulation result, and the prediction accuracy of the flight time analytic solution of the present invention is verified.

In the second document, the formula of the flight time under the condition of constant longitudinal lift-drag ratio is

Setting an initial longitude λ of an aircraft₀0deg, initial latitude phi₀Initial energy E of 50deg₀＝-3.8602×10⁴kJ/kg and terminal energy E_f＝-5.5×10⁴kJ/kg. Consider 5 different flight directions: psi₀100deg, 180deg, -100deg, 20deg, -20 deg. Setting the longitudinal lift-drag ratio of the second document to

Transverse lift-to-drag ratio of

The simulation results are shown in table 1. As can be seen from the simulation results, the flight times of the aircraft in different directions are not the same due to the effect of the earth rotation, but the analytic solution of the second document fails to take into account the time difference caused by the earth rotation. By comparing with the trajectory simulation result, it can be seen that: because the influence of the analytic solution on the earth rotation is compensated, the precision is higher, and the prediction error is kept within 3 percent. In terms of computational efficiency, the computation time of the time resolution solution is at least 5 orders of magnitude smaller than that of the trajectory simulation.

TABLE 1

Example two

The embodiment is an ideal and non-interference situation, and solves the problem of cooperative flight of guiding a plurality of hypersonic gliding aircrafts (V1, V2, V3) to be launched from different places but simultaneously reach the same target in a multi-no-fly zone environment. 64 circular no-fly zones with the radius of 200km are arranged in the flight area, the initial conditions of the aircraft are shown in table 1, and the target point position is longitude lambda_T130deg, latitude phi_T-20 deg. T can be selected by using the multi-aircraft arrival time coordination scheme of step 7_TAEM2900s and determines the re-entry start time for each aircraft, 109.29s, 358.69s and 751.62s respectively.

TABLE 2

The horizontal distance from the aircraft to the target in the reentry section is S_TAEMAnd ends when the speed is 50 km. The desired terminal height at this time is H_TAEM25km, terminal speed V_TAEM2000m/s, the terminal heading error satisfies | Δ ψ_TAEMLess than or equal to 5deg and the terminal inclination angle satisfies sigma_TAEM|≤30deg。

Fig. 4 shows that the aircraft successfully avoids all the no-fly zones and reaches the target under the control of the guidance method of the invention. FIG. 5 is a height-velocity curve, where H_minThe curve is the lower bound of height determined by the process constraints. As can be seen, the aircraft is susceptible to flying near the lower altitude boundary due to the harsh thermal environment during the high speed phase. It can also be seen that the terminal speed and height meet the terminal requirements. Fig. 6 shows the angle of attack law of an aircraft, from which it can be seen that three aircraft have different re-entry starting moments. Fig. 7 shows a roll angle profile. FIG. 8 shows the reference remaining flight distances for three casesAn off-plane. Fig. 9 shows the reference time-of-flight profiles for three cases.

EXAMPLE III

In the embodiment, Monte Carlo simulation verifies the robustness of the guidance method under the condition that the aircraft model has the pull bias. Wherein, the pneumatic coefficient adopts the following linear deviation model

Wherein the content of the first and second substances,_CLand_CDthe lift coefficient and drag coefficient deflection percentages, respectively, vary with angle of attack α and mach number Ma._CL0、

Is a relevant zero-mean normal distribution random parameter. The wind power and atmospheric density perturbation model is as follows

_ρ＝k_{ρ ρ(max)}(130)

Wherein the content of the first and second substances,

is the wind speed in the east-west direction,

is the wind speed in the north-south direction,_ρis the percent deviation from atmospheric density. V_wind(max)Is the maximum possible wind speed, varies with altitude, and can reach 170m/s at high altitude._ρ(max)Is the maximum possible deviation percentage of the atmospheric density and can be realized at high altitudeUp to 50%.

And k_ρAre the corresponding zero mean normal distribution random parameters. The statistical properties of the random parameters and initial condition errors are shown in table 3.

TABLE 3

Consider here the scenario of two aircraft (V1, V2) launched from two different locations for coordinated flight, nominal conditions as shown in Table 4, with a target position of longitude λ_T120deg and latitude phi _T0 deg. The no-fly zone distribution is shown in fig. 10. Setting the expected arrival time to 2900s, nominal re-entry segment start times 426.68s and 404.29s may have been determined using a multi-aircraft arrival time coordination scheme. Note that the nominal trajectory, as well as the remaining flight distance and time-of-flight reference profile required for guidance laws, are obtained based on the ideal aerodynamic model and the ideal atmospheric model, and do not take into account, nor do they take into account, any on-line measured pull bias data, since they are calculated off-line in advance.

TABLE 4

Fig. 10 shows the reentry trajectory of the aircraft, and it can be seen that in an environment with uncertain disturbance, the guidance method can successfully guide the aircraft to avoid the no-fly zone and successfully reach the target. FIG. 11 shows a height-velocity curve, where H_minIs a lower bound on height determined by process constraints such as heat flow density. As can be seen from the figure, the guidance method can ensure the flight safety. Fig. 12-13 illustrate an angle of attack curve and a roll angle curve. FIG. 14 shows the distribution of terminal velocity and terminal heading error, both within the required range. FIG. 15 shows the statistics of the arrival time distribution, the percentage of the arrival time error of the aircraft V1 within + -5 s is 79.01%; the ratio of the arrival time error of the aircraft V2 within +/-5 s is74.36％。

Claims (3)

1. A collaborative analysis reentry guidance method considering multi-forbidden flight zone constraint is characterized in that: the method comprises the following steps:

step 1: reentry guidance problem description

The method comprises the steps of adopting a spherical three-dimensional reentry particle dynamic model for rotating the earth, wherein longitude lambda, latitude phi and altitude H are used for describing aircraft position information, and velocity V, ballistic inclination angle gamma and heading angle psi are used for describing velocity information;

in order to ensure that each subsystem of the aircraft works normally, the reentry flight trajectory needs to meet the heat flux density

Incoming flow pressure q and overload n are restrained; considering that the flight control system has limited capability, the range and the change rate of the attack angle and the roll angle need to be limited; in addition to the conventional constraints, the case of multiple circular no-fly zone constraints needs to be considered;

the horizontal distance from the aircraft to the target in the reentry section is S_TAEMThen the process is terminated; the desired terminal height at this time is H_TAEMTerminal velocity is V_TAEMTerminal heading error | Δ ψ_TAEMAnd terminal roll angle σ_TAEMI satisfies a constraint condition; in addition, the desired terminal time is also designated as t_TAEM；

Step 2: based on the time-of-flight analytic solution of the rotating earth model,

to deal with the arrival time constraint, here a time-of-flight resolution is derived; the energy is defined as follows

Wherein V is velocity, H is height, R_eIs the earth's mean radius, μ is the gravitational constant; the derivative of the energy with respect to time t is

Wherein g is gravity acceleration; in the context of the rotating earth, there are complex non-linear equations

Where γ is the ballistic inclination, D is the drag, m is the mass, ω is_eIs the earth rotation angular velocity, phi is the latitude, psi is the heading angle;

right item 3 in formula (4)

Item 4

Is an inertial force component caused by the rotation of the earth; the component of the inertial force in the direction of the velocity is regarded as an additional resistance, as follows

Defining an equivalent resistance of

Substituting the formulas (3) to (5) into the formula (2), and taking the reciprocal thereof to obtain

Defining a longitudinal lift-to-drag ratio of

L therein₁L cos (σ) is the component of lift in the longitudinal plane, plan

Is composed of

Wherein a is₀,a₁,a₂Is a quadratic polynomial coefficient;

now make use of

The profile estimates the equivalent drag in the steady glide state, which is recorded

Time derivative of trajectory inclination as follows

Suppose γ ≈ 0 and

estimate the steady glide longitudinal lift as follows

Wherein, Delta L₁Is the component of the inertia force along the vertical direction and is regarded as the additional longitudinal lift force, and the formula is as follows

At the same time, with the above assumptions, Δ D is also reduced to

Then

Estimated from the following equation

Substituting formula (12) into formula (6), and substituting E for V to obtain

Because H < R_eSo that H is the intermediate value of glide height H^*Further let R^*＝R_e+H^*(ii) a Although the inertial force is a small quantity, when the velocity is close to the first cosmic velocity, the inertial force may cause the denominator of the above equation to be zero, so that singularity occurs; to avoid the above singularity, the inertial force is considered as a small quantity and a first order Taylor expansion is performed as follows

Wherein the simplification symbol h_z1、h_mIs defined as follows

h_m＝2E+μ/R^*(16)

According to Δ L₁And Δ D, expression of_z1Is divided into two parts h_PAF1And h_PAF2As follows

h_z1＝h_PAF1+h_PAF2(17)

Wherein the content of the first and second substances,

h_PAF1＝-2R^*Vω_ecos(φ)sin(ψ) (18)

due to h_PAF2Is much less than h_PAF1Thus, a linear function pair h is used_PAF2Carry out an approximation as follows

h_PAF2(E)≈k_PAF2(1)V+k_PAF2(0)(20)

Wherein the coefficient k_PAF2(1)And k_PAF2(0)Determined by two endpoints, as follows

Wherein h is_PAF2(E) Calculated by substituting the current state into equation (19), and h_PAF2(E_TAEM) The expected terminal state is obtained by substituting the expected terminal state into the formula (19);

actually, the hypersonic gliding aircraft flies to the target around one great circle marked as the generalized equator, the aircraft transversely maneuvers around the generalized equator, and the azimuth angle of the generalized equator

Considered as some weighted average of the aircraft heading angle, and therefore, using the azimuth of the generalized equator with the subsequent heading angle curve not yet known

Substitution of h_PAF2(E) And h_PAF2(E_TAEM) Course angle in the formula; up to this point, in the formula (14), except for the independent variable E, other parameters are constant values, and the formula (14) is integrated; through derivation and arrangement, an expression of a flight time analytic solution is obtained as follows:

wherein the coefficient k_t(1)、k_t(2)、k_t(3)、k_t(4)、k_t(5)、k_t(6)、k_t(7)Is expressed as follows

And step 3: analytic re-entry guidance method considering multiple no-fly zones and arrival time constraints

Aiming at a single aircraft, designing a reentry guidance method considering multiple no-fly zones and arrival time constraints; according to ballistic characteristics, the reentry process is generally divided into three phases: a descending stage, a stable gliding stage and a height adjusting stage;

s31: descending section

In order to avoid overlarge heat flux density caused by falling into a dense atmosphere, the aircraft glides at the maximum available attack angle and zero inclination angle in the descent segment; when the lift force is enough to support the aircraft to glide stably, the attack angle is smoothly transited to the reference attack angle, and then the aircraft enters a stable gliding stage;

s32: steady glide phase

The steady gliding stage is the longest, most important and most complex reentry flight stage, and a guidance scheme in the stage rapidly plans a reference trajectory meeting the constraint of a no-fly area and arrival time on line by utilizing a three-dimensional reentry trajectory analytic solution and a flight time analytic solution;

s33: height adjustment stage

After the last reversal point, the aircraft enters the altitude adjustment phase, but at some node before the last reversal point, the aircraft starts the altitude adjustment phase preparation; wherein, the energy of the aircraft corresponding to the last reversal point is recorded as E_BR(nTR)The energy of the aircraft corresponding to a node before the last reversal point is denoted as E_ST；

Designing a multi-target numerical iteration planning scheme based on-line trajectory simulation; in the scheme, a quasi-Newton method iteration algorithm is improved by using a directional derivative;

when E ═ E_STThen, the subsequent track is finely adjusted accurately by using the online trajectory simulation-based multi-target numerical value iterative programming scheme so as to meet the requirements of terminal time, speed and height; when E is_st＞E＞E_BR(nTR)When the method is used, a collaborative analysis re-entry guidance method considering the constraint of multiple forbidden flight areas tracks the longitudinal lift-drag ratio, the residual flight distance and the reference profile of the flight time relative energy obtained in the iterative planning, and when the E is used_BR(nTR)When the angle is more than E, determining a reference roll angle by using a proportional guidance law to eliminate course errors, and tracking the flight time and the residual flight distance reference profile by finely adjusting an attack angle to ensure that the requirements of terminal time and speed are met;

and 4, step 4: multi-aircraft arrival time coordination scheme

Suppose there is n_HGVAn aircraft, denoted as HGV_i,i＝1,2,…,n_HGV(ii) a After the launching point and the target point are given, the rocket is boosted according to the performance of the rocket and the selected initiativeDetermining the initial state of a reentry section of the aircraft according to a section guidance scheme; then, determining corresponding longitudinal lift-drag ratio profile parameters by using a longitudinal analytical solution only according to energy management requirements; furthermore, the off-line trajectory simulation is used for predicting the flight time t of all the aircrafts_EF(HGVi)，i＝1,2,…,n_HGVSelecting a flight time t with the longest time_EF(max)And reserve a transmission preparation time t_preAnd boost time of flight t_boostDetermining the expected time of arrival t_TAEM＝t_pre+t_boost+t_EF(max)And further use t_TAEM-t_EF(HGVi)Determining the re-entry period starting time of each aircraft; then, according to the selected boosting section guidance method, the corresponding emission time is determined.

2. The cooperative analytic re-entry guidance method considering multiple no-fly zone constraints as set forth in claim 1, wherein: step S32: the guidance process in the steady gliding stage is as follows:

s321, designing a reference attack angle section, and determining a corresponding reference lift-drag ratio section;

s322, reasonably designing parameterized longitudinal lift-drag ratio by considering energy management requirements and arrival time constraints

While determining a parameterized lateral lift-to-drag ratio

Cross section of (a);

s323, in order to compensate the effect of the earth rotation

And

on the basis of the section, the equivalent longitudinal and transverse lift-drag ratios are established

And

cross section of (a);

s324, solving longitudinal profile parameters by utilizing a longitudinal analytical solution and the flight time analytical solution obtained in the step 2 according to energy management and arrival time requirements;

s325, updating a tilting reversal sequence by utilizing a longitudinal analysis solution and a transverse analysis solution at intervals of 60S according to the no-fly zone and the terminal position constraint;

s326, adjusting a reference roll angle, tracking an equivalent longitudinal lift-drag ratio section, and reversing according to the planned roll;

s327, in order to restrain ballistic oscillation, ballistic damping feedback is introduced into a reference attack angle and a roll angle, and therefore a guidance instruction is obtained;

s328, limiting the roll angle instruction according to process constraints;

and S329, repeating the process from the step S322 until the smooth gliding phase is finished.

3. The cooperative analytic re-entry guidance method considering multiple no-fly zone constraints as set forth in claim 1, wherein: the multi-target numerical iterative programming scheme based on the online trajectory simulation described in step S33 includes the following specific steps:

s331, predicting the terminal time and the terminal state by using trajectory simulation;

s332, judging whether the terminal speed or the terminal time meets the requirement; if yes, skipping to step S333, otherwise skipping to step S334;

s333, adjusting a longitudinal lift-drag ratio profile and a last inversion point by using a quasi-Newton iteration method according to the terminal speed and the terminal time requirement, and jumping to the step S331;

s334, judging whether the height of the terminal meets the requirement; if yes, jumping to step S336, otherwise jumping to step S335;

s335, correcting the reference attack angle parameter according to the terminal height error, and skipping to the step S331;

and S336, finishing the iterative planning.

Images