CN106843272A - A kind of explicit Guidance rule with terminal velocity, trajectory tilt angle and overload constraint - Google Patents

Abstract

The present invention is a kind of explicit Guidance rule with terminal velocity, trajectory tilt angle and overload constraint, is formed by trajectory shaping two synthesis of Guidance Law and terminal velocity control program.Trajectory shaping Guidance Law can control aircraft from predetermined direction hit, and terminal velocity control program crossrange maneuvering acceleration then by controlling aircraft to go for a stroll, terminal velocity size is adjusted so as to the distance for controlling extension to fly, crossrange maneuvering acceleration magnitude is determined by iterated revision algorithm.The advantage of invention is：The parsing Guidance Law can not only meet terminal trajectory tilt angle, can also meet terminal velocity constraint, and aircraft is gradually decayed to 0 close to the Maneuver Acceleration of target；The method for determining Guidance Law coefficient is further provided, i.e., Guidance Law coefficient is determined by the appropriate characteristic root for choosing linear approximation system；Obtain Guidance Law coefficient stabilization domain, as long as Strict Proof Guidance Law coefficient is in stable region, guidance system stabilization and aircraft is with Low Angle Of Attack hit.

Description

A kind of explicit Guidance rule with terminal velocity, trajectory tilt angle and overload constraint

Technical field

The present invention relates to a kind of explicit Guidance rule with terminal velocity, trajectory tilt angle and overload constraint, belong to space flight skill Art, weapon technologies, Guidance and control field.

Background technology

For the terminal flight of hypersonic glide vehicle, ordered from approximately perpendicular direction it is generally desirable to directing aircraft Middle ground target, so that target seeker has a preferable visual field；Simultaneously, it is necessary to control the terminal velocity of aircraft, because one As in the case of excessive terminal velocity be unfavorable for ensureing heat flow density and to flow pressure constraint that too small terminal velocity is unfavorable for dashing forward The need for anti-.If what aircraft was loaded is earth penetrator, it is required that with less angle of attack hit, otherwise, when with larger During angle of attack hit, body is subject to huge asymmetric power to act in Penetration so that there is large curved in penetration path, And Penetration Depth is reduced significantly.

The Typical Representative of Guidance Law is shaped as trajectory, explicit Guidance rule is carried by Cherry the sixties in 20th century earliest Go out.Cherry G.,A general,explicit,optimizing guidance law for rocket-propelled Spaceflight (a kind of explicit optimal guidance law of broad sense suitable for rocket-powered vehicle) [C] .AIAA/ION Astrodynamics Guidance and Control Conference, 1964. by assuming that instruction thrust vectoring is the time Polynomial curve, be deduced explicit Guidance rule.This Guidance Law can control airship, and vector reaches predetermined at a predetermined rate Point.Explicit Guidance rule is applied on Apollo later.

Lin,C.F.,Tsai,L.L.,Analytical solution of optimal trajectory-shaping Guidance (trajectory shaping Guidance Law analytic solutions) [J] .Journal of Guidance, 1987,11 (1):Will be aobvious in 61-66 Formula Guidance Law is extended to and only provides on the aircraft of controling power on velocity attitude, and combines aircraft characteristic, right Explicit Guidance rule parameter is optimized.Explicit Guidance rule after improvement is used for controlling terminal velocity direction.

Zarchan, P., Tactical and Strategic Missile Guidance, 5th ed. (tactics and strategy Missile guidance is restrained) [M], AIAA Progress in Aeronautics and Astronautics, Virginia, US, With energetic optimum as target in 2007.Chap.25, it is deduced explicit Guidance using Schwartz inequality and restrains.

Ohlmeyer E.J., Phillips C A.Generalized Vector Explicit Guidance (broad sense Speed explicit Guidance is restrained) [J] .Journal of Guidance Control＆Dynamics, 2015,29 (2):In 261-268 Then by improving the object function in optimal control problem, the span that explicit Guidance restrains coefficient is extended.

In practical problem, some also need to carry out deceleration control using aerodynamic force as the reentry vehicle of major control power System.For example, Pan Xing II ballistic missiles reenter bullet in terminal guidance stage early stage by pull-up reentry trajectory come consumed energy, from And terminal velocity size is controlled roughly, in the later stage, then shape Guidance Law directing aircraft using trajectory and hit mesh from vertical direction Mark.During Shuttle reentry, energy management mainly is carried out in gliding section, and also entered in latter end energy management section (TAEM) The a range of terminal velocity control of row.At TAEM sections, S types are motor-driven and regulation is made laterally doing by directing aircraft for Guidance Law Dynamic plate tracks resistance-distance Curve, to control speed.

Traditional explicit Guidance rule has narrower coefficient stabilization domain, it is made up of discrete point range, and Guidance Law is pushed away More constraints is led, such as speed is definite value.Meanwhile, traditional explicit Guidance rule can not meet to terminal speed simultaneously The constraint of degree, trajectory tilt angle and overload, larger overload can be significantly increased miss distance during missile target encounter.The present invention is by spectrum point Solution method obtains the Generalized Analytic solution of explicit Guidance rule, obtains bigger Guidance Law index layer stable region, and integrate end In Guidance Law, foring can be while meets the aobvious of the constraint of terminal velocity, trajectory tilt angle and overload for spot speed control program Formula Guidance Law.

The content of the invention

The invention aims to solve the above problems, propose that one kind has terminal velocity, trajectory tilt angle and overload about The explicit Guidance rule of beam.This Guidance Law can control unpowered vehicle the terminal guidance stage at a predetermined rate, from predetermined Ground target is hit in direction, and is met the motor-driven overload of terminal and tended to 0.

Explicit guidance rule is formed by trajectory shaping two synthesis of Guidance Law and terminal velocity control program.Trajectory is into shape Leading rule can control aircraft from predetermined direction hit, and terminal velocity control program is then gone for a stroll by controlling aircraft Crossrange maneuvering acceleration so that the distance for controlling extension to fly adjusts terminal velocity size, crossrange maneuvering acceleration magnitude Determined by iterated revision algorithm.In order to not disturb trajectory to shape Guidance Law, crossrange maneuvering acceleration is with aircraft close to target Gradually decay to 0.

In order that the motor-driven overload of terminal is zero, the present invention have studied the linear approximation system under trajectory shaping Guidance Law effect System, and the analytic solutions for instructing motor-driven overload are obtained, the relation of Guidance Law coefficient and the stability of a system is analyzed, obtain guidance system Number stable region.

The action effect of this explicit guidance rule is verified by the example of CAV-H.

The present invention is a kind of explicit Guidance rule with terminal velocity, trajectory tilt angle and overload constraint, is specifically included as follows Step：

Step 1：Set up kinematics and dynamics equation

The inertial reference system o-xyH connected firmly with ground is set up, it is longitudinal range to regard aircraft as particle M, x, and y is laterally Range, H is height above sea level.V is the speed of aircraft.γ is the trajectory tilt angle of aircraft, and relative level reference line turns counterclockwise Move as just.ψ is course angle, i.e. the horizontal component of aircraft speed and the angle of x-axis, is rotated counterclockwise as just relative to x-axis.This Invention does not consider the influence of earth curvature and rotation, and its kinematics and dynamics equation group is as follows：

Wherein, t is the aircraft flight time, and m is vehicle mass, and σ is roll angle, and L is lift, and D is resistance, and g attaches most importance to Power acceleration.

Step 2：Explicit Guidance rule summary of the present invention

Under explicit Guidance rule of the present invention effect, guidance system constitutes the command acceleration of generation by two：Section 1 It is the command acceleration a of trajectory shaping Guidance Law generation_TSG(TSG:Trajectory Shaping Guidance), this acceleration Can directing aircraft from predetermined direction hit；Section 2 is the command acceleration that terminal velocity control program is produced a_speed, this accelerated energy control aircraft does crossrange maneuvering, so that velocity magnitude when adjusting hit.Guidance system is produced Command acceleration expression formula it is as follows

a_cmd=a_TSG+a_speed (7)

Step 3：Solve the analytic solutions that trajectory shapes Guidance Law

The command acceleration a that Section 1 trajectory shaping Guidance Law is produced in the step 2_TSGIt is made up of three events, is used respectively With directing aircraft target, velocity attitude during control aircraft hit, balancing gravity acceleration is perpendicular to speed Degree durection component.Acceleration a_TSGDirection perpendicular to aircraft current velocity vector, its expression formula is shown in formula (8)

Wherein, R is remaining flying distance；C₁And C₂It is Guidance Law coefficient；It is aircraft to the unit direction of line of sight Vector；It is the unit direction vector of aircraft speed；Predetermined speed direction vector when being aircraft hit, by formula (9) calculate；g_nIt is component of the acceleration of gravity perpendicular to speed.

Wherein, γ_fThe predetermined trajectory tilt angle of aircraft, is typically in the range of between -70 ° and -90 ° during for hit, is normal Value；ψ_fIt is as follows equal to the azimuth of current flight device to line of sight

Wherein,WithVector is represented respectivelyIn the component in reference axis x, y direction.

Ignore the influence of gravity, then in fore-and-aft plane, the expression formula of aircraft guidance rule can be by geometric vector form It is transformed into trigonometric function form：

Wherein, γ_LOSIt is the visual line angle of bullet, relative level reference line is rotated counterclockwise as just.

Sight line angular rate of change is：

Assuming that γ ≈ γ_LOS≈γ_f, then the trigonometric function line in formula (11) (12), and arrangement can be obtained such as lower linear Time-varying system

Wherein

B₁=[C₂ 0]^T (16)

Wherein γ and γ_LOSIt is state variable, γ_fIt is control variables, subscript " T " representing matrix here and in full text Transposition.

The present invention carries out the solution of linear time varying system using a new method based on spectral factorization.

Definition

Wherein t₀It is the initial time of system, and has

Based on hypothesis (γ-γ above_LOS) ≈ 0, can approximately there are dR=-Vd τ.Substituting into above formula can solve

At formula (13) left and right two ends while premultiplication Q (t, t₀)

Above formula is rewritable to be

It is reverse to be obtained using step integration rule

Both sides integrate simultaneously

Wherein, γ₀It is initial trajectory inclination angle, γ_LOS0It is the visual line angle of initial bullet.exp(-A₁f₂(t,t₀))=exp (0_2×2)=I_2×2It is the unit matrix of 2 × 2,0 in formula_2×2It is 2 × 2 null matrix.

Q(t,t₀) inverse matrix it is as follows

Φ(t,t₀)=[Q (t, t₀)]^-1=exp (A₁f₂(t,t₀)) (24)

Wherein, Φ (t, t₀) it is called state-transition matrix.

At formula (23) two ends while premultiplication Φ (t, t₀) can obtain

Matrix A can be obtained by formula (15)₁Proper polynomial be

|λI-A₁|=λ²+(C₁+C₂-1)λ-C₂ (26)

Can be in the hope of matrix A₁Characteristic value be

Wherein,

On the contrary, Guidance Law coefficient C can also be obtained₁And C₂On matrix A₁Two expression formulas of characteristic value

Define f₃(x,t,t₀)=exp (xf₂(t,t₀)), formula (19) is substituted into can be obtained

Due to C₁And C₂It is real number, so matrix A₁Two characteristic values can only be two not wait real number, two equal real numbers With the one kind in three kinds of situations of complex conjugate.Here by matrix A₁Two characteristic values be divided into two kinds of situations and discuss, one kind is λ₁ ≠λ₂, another kind is λ₁=λ₂。

(I) λ, is worked as₁≠λ₂When

Can be obtained by spectral factorization formula

Wherein, G₁And G₂It is matrix A₁Spectrum matrix, and have

Wherein, I is unit matrix.Substitution formula (15) and formula (27) are solved

Formula (31) substitution formula (25) right-hand member Section 2 is obtained

Then by formula (25) (31) (34) above, can obtain working as λ₁≠λ₂When trajectory tilt angle γ (t) analytic solutions so that ByThe analytic solutions of Maneuver Acceleration can be obtained, it is as follows

Now prove to work as t → t_fWhen, there is the λ met under certain condition₁And λ₂So thatNeed to A₁Two Eigenvalue λ₁And λ₂For two, grade real number, a pair of conjugate complex numbers, two kinds of situations enter line justification respectively.

a)λ₁And λ₂For two not wait real number

Due to terminal juncture t_fWhen, missile target encounter has R (t_f)=0.Following two conclusions then can be obtained by formula (35)：

(1) if λ₁≤-1,λ₂≤ -1 and λ₁≠λ₂, then as t → t_fWhen, a_LDo not dissipate；

(2) if λ₁＜ -1, λ₂＜ -1 and λ₁≠λ₂, then as t → t_fWhen,

b)λ₁And λ₂It is a pair of Conjugate complex roots

If λ₁And λ₂For

Wherein, p and q are real numbers,Above formula is substituted into formula (29) can obtain

Formula (36) (37) substitution formula (35) is obtained into acceleration analytic solutions expression formula is

As can be seen from the above equation, accelerating curve has concussion and occurs, and as t → t_fWhen, concussion frequency tends to infinitely great, Because now R (t_f)=0.But if p ＜ -1, then as t → t_fShi You

Demonstrate,prove, as Conjugate complex roots λ₁And λ₂When meeting real part p ＜ -1, have

(II) λ, is worked as₁=λ₂When, i.e. Δ=0

It is easily verified that

Then may determine that matrix A₁Minimum formula m (x) be

M (x)=(x- λ₁)² (42)

Then obtained by broad sense spectrum formula

Wherein

Formula (43) is substituted into the integral term of formula (25) right-hand member, can be obtained

Wherein,

Can then be obtained by formula (25) (43) (45), work as λ₁=λ₂When trajectory tilt angle γ (t) analytic solutions, so as to byThe analytic solutions of Maneuver Acceleration can be obtained, it is as follows

Prove to work as λ now₁=λ₂During ＜ -1,Obvious this problem equivalent is in the following problem of proof：

Work as λ₁During ＜ -1

Prove：

Due to R (t_fλ is worked as in)=0, easily card₁During ＜ -1,

Second limit is oneThe type limit, can be obtained by L ' Hospital methods

Work as λ so as to demonstrate₁=λ₂During ＜ -1,

Step 4：Analysis trajectory shaping Guidance Law coefficient C₁And C₂Stable region

It is given to make the characteristic root λ of the line approximation system of guidance system stabilization₁And λ₂, trajectory is then determined by characteristic root Shaping Guidance Law coefficient C₁And C₂.To characteristic root λ₁And λ₂Discuss in two kinds of situation.

Ⅰ.λ₁And λ₂It is real number

Analyzed by step 3, work as λ₁≤-1,λ₂≤ -1 and λ₁And λ₂When being asynchronously -1, guidance system stabilization, and machine Dynamic acceleration instruction does not dissipate.Thus formula (29) can be rewritten as following form

If C₂Definite value, then C₁On λ₁Derivative be

It is easy to get, whenWhen, C₁Minimum value is taken, i.e.,

Work as λ₁=-1 or λ₁=C₂When, C₁Maximum is taken, i.e.,

C₁≤2-2C₂ (52)

In sum, λ₁And λ₂For the coefficient value scope of trajectory shaping Guidance Law in the case of real number is

Ⅱ.λ₁And λ₂It is a pair of Conjugate complex roots

Defined function Re (x) is the real part of plural x, and being analyzed by step 3 to obtain, as Re (λ₁)=Re (λ₂) ＜ -1 when, system Guiding systems stabilization, and Maneuver Acceleration instruction finally converges to zero.Can then be obtained by formula (27)

Arrangement can be obtained

I.e.

Easily demonstrate,proveTherefore inequality group (56) can be reduced to

MakePerseverance is set up, then must have C₂＜ -1.

In sum, can obtain：

If C₁=2-2C₂And C₂＜ -1, then TSG guidance systems stabilization, and in aircraft close to acceleration during target Instruction does not dissipate.If 3-C₂＜ C₁＜ 2-2C₂And C₂＜ -1, then TSG guidance systems stabilization, and aircraft close to target mistake Acceleration instruction converges on 0 in journey.

C₁And C₂Span extend to plane domain, and Guidance Law by the zone of dispersion of traditional trajectory shaping Guidance Law Coefficient C₁And C₂Span not by aircraft speed change influenceed.

Step 5：Terminal velocity control program

The command acceleration a that Section 2 terminal velocity control program described in step 2 is produced_speed, this acceleration effect exists In local level, and perpendicular to speed, it is shown below

Wherein, k_VfIt is positive parameter undetermined, g₀It is the acceleration of gravity of sea level altitude, R₀It is start time aircraft and mesh Remaining flying distance between mark, ψ₀It is the course angle of initial time, ψ_LOS0It is the sight line azimuth of initial time.Sgn (x) is Sign function, it is as follows

By adjusting k_VfValue, aircraft crossrange maneuvering amplitude can be controlled, so that when adjusting aircraft hit Speed.a_speedWith a_TSGThe latter end of aircraft flight is acted on simultaneously.a_speedSize with aircraft close to target gradually Weaken, can so avoid a in aircraft soon hit_speedTo a_TSGInterference.

Aircraft predicted the speed V of hit by Ballistic Simulation of Underwater before the terminal guidance stage is entered_f(K), then press Secant methods correct k_VfValue, finally obtain k_VfIdeal value, it is as follows

Wherein, V_fcmdIt is desired terminal velocity, K is iterations.

The advantage of the invention is that：

(1) compared to traditional Guidance Law, the parsing Guidance Law can not only meet terminal trajectory tilt angle, can also meet Terminal velocity is constrained, and causes that aircraft gradually decays to 0 close to the Maneuver Acceleration of target；

(2) the Generalized Analytic solution of explicit Guidance rule is obtained, the method for determining Guidance Law coefficient is further provided, that is, is led to The appropriate characteristic root for choosing linear approximation system is crossed to determine Guidance Law coefficient；

(3) relation of Guidance Law coefficient and the stability of a system is analyzed, Guidance Law coefficient stabilization domain, Strict Proof is obtained As long as Guidance Law coefficient is in stable region, guidance system stabilization and aircraft is with Low Angle Of Attack hit.

Brief description of the drawings

Fig. 1 is inertia reference coordinate system o-xyH and related state variable.

Fig. 2 is the geometrical relationship in fore-and-aft plane.

Fig. 3 (a), (b) are the simulation results that characteristic value is the trajectory of the present invention shaping Guidance Law in the case of a pair of real numbers.

Fig. 4 (a), (b) are the emulation knots that characteristic value is the trajectory of the present invention shaping Guidance Law in the case of a pair of conjugate complex numbers Really.

Fig. 5 (a), (b) are the emulation knots that characteristic value is the trajectory of the present invention shaping Guidance Law in the case of a pair equal real numbers Really.

Fig. 6 is stable region and the traditional explicit Guidance Law coefficient stabilization domain of trajectory shaping Guidance Law coefficient of the present invention.

Fig. 7 (a), (b) are the aerodynamic data curves of aircraft.

Fig. 8 is terminal velocity V_fWith k_vfChange curve.

Fig. 9 (a), (b) (c), (d), (e), (f) are that explicit guidance's rule of the present invention is directed to two kinds of different terminal velocity situations Simulation result.

Figure 10 is that the instruction of explicit Guidance rule of the present invention is produced and uses flow chart.

Specific embodiment

Below in conjunction with drawings and Examples, the present invention is described in further detail.

The present invention is a kind of explicit Guidance rule with terminal velocity, trajectory tilt angle and overload constraint, it is adaptable to air to surface The final guidance rule of aircraft, this Guidance Law enables to aircraft to be wanted in the constraint that end of time meets trajectory tilt angle and speed Ask, and cause that the motor-driven overload of terminal is 0.

This Guidance Law is formed by trajectory shaping two synthesis of Guidance Law and terminal velocity control program.Trajectory shapes Guidance Law It is that can control aircraft from predetermined direction hit, when terminal velocity control program is used for controlling aircraft hit Velocity magnitude.There is an adjustable parameter in terminal velocity control program, aircraft can be adjusted and gone for a stroll distance, so as to control terminal Speed, parameter size is produced by missile-borne computer Ballistic Simulation of Underwater iteration convergence.Fly to not disturb trajectory to shape Guidance Law guiding Row device hit, terminal velocity control program acts predominantly on the early stage of aircraft terminal flight, and as aircraft is close Target and weaken.Meanwhile, studied by the analytic solutions to line system, it is 0 to be met terminal and instruct motor-driven overload Guidance Law coefficient value scope.As long as Guidance Law coefficient value within this range, can just reach with the mesh of Low Angle Of Attack hit 's.As shown in Figure 10, whole process includes following steps：

Step 1：Set up kinematics and dynamics equation

Fig. 1 describes the inertial reference system o-xyH for using of the invention and in the undefined relevant state variables of this coordinate system. O-xyH is mutually connected firmly with ground, and it is longitudinal range to regard aircraft as particle M, x, and y is horizontal range, and H is height above sea level.V is winged The speed of row device.γ is the trajectory tilt angle of aircraft, and relative level reference line is rotated counterclockwise as just.ψ is course angle, that is, fly The horizontal component of device speed and the angle of x-axis, rotate counterclockwise as just relative to x-axis.The present invention do not consider earth curvature and from The influence for turning, its kinematics and dynamics equation group is as follows：

Wherein, t is the aircraft flight time, and m is vehicle mass, and σ is roll angle, and L is lift, and D is resistance, and g attaches most importance to Power acceleration, they are calculated by following formula.

L=C_LqS (7)

D=C_DqS (8)

Wherein, C_LAnd C_DIt is respectively lift coefficient and resistance coefficient.Q flows pressure.S is pneumatic area of reference.μ is normal Number, value 3.96272 × 10¹⁴m³/s²。R_eIt is earth mean radius, size is 6356.766km.

Step 2：Explicit Guidance rule summary of the present invention

Under explicit Guidance rule of the present invention effect, guidance system constitutes the command acceleration of generation by two：Section 1 It is the command acceleration a of trajectory shaping Guidance Law generation_TSG(TSG:Trajectory Shaping Guidance), this acceleration Can directing aircraft from predetermined direction hit；Section 2 is the command acceleration that terminal velocity control program is produced a_speed, this accelerated energy control aircraft does crossrange maneuvering, so that velocity magnitude when adjusting hit.Command acceleration table It is as follows up to formula

a_cmd=a_TSG+a_speed (10)

Step 3：Solve the analytic solutions that trajectory shapes Guidance Law

The command acceleration a that Section 1 trajectory shaping Guidance Law is produced in the step 2_TSGIt is made up of three events, is used respectively With directing aircraft target, velocity attitude during control aircraft hit, balancing gravity acceleration is perpendicular to speed Degree durection component.Acceleration a_TSGDirection perpendicular to aircraft current velocity vector, its expression formula is shown in formula (11)

Wherein, R is remaining flying distance；C₁And C₂It is Guidance Law coefficient；It is aircraft to the unit direction of line of sight Vector；It is the unit direction vector of aircraft speed；Predetermined speed direction vector when being aircraft hit, by formula (12) calculate；g_nIt is component of the acceleration of gravity perpendicular to speed.

Wherein, γ_fThe predetermined trajectory tilt angle of aircraft, is typically in the range of between -70 ° and -90 ° during for hit；ψ_fIt is equal to Current flight device is as follows to the azimuth of line of sight

Wherein,WithVector is represented respectivelyIn the component in reference axis x, y direction.

End is set to overload the Guidance Law coefficient C for 0 in order to calculate₁And C₂Stable region, it is necessary first to derive trajectory shape The analytic solutions form of Guidance Law.Because aircraft is in hit, speed is approximately perpendicular to horizontal plane, so it is non-to work as aircraft During very close to target, g_nMould it is very small, can be ignored.What Fig. 2 was represented is when the influence of gravity is not considered, in longitudinal direction Belligerent geometrical relationship in plane, wherein, M represents aircraft, and T represents target.Guidance Law expression formula can be by geometric vector shape Formula is transformed into trigonometric function form：

Wherein, γ_LOSIt is the visual line angle of bullet, relative level reference line is rotated counterclockwise as just.

Sight line angular rate of change is：

Assuming that γ ≈ γ_LOS≈γ_f, then the trigonometric function line in formula (14) (15), and arrangement can be obtained such as lower linear Time-varying system

Wherein

B₁=[C₂ 0]^T (19)

Wherein γ and γ_LOSIt is state variable, γ_fIt is control variables, subscript " T " representing matrix here and in full text Transposition.

It is clear to, the equalization point of said system meets：γ=γ_LOS=γ_f, so as to prove to assume to set up at equalization point.

For traditional Linear Time-Invariant System, the method that can change with Laplace is solved.But system (16) is Linear time varying system, the method is no longer applicable.Therefore, the present invention is solved using a new method based on spectral factorization.

Definition

Wherein t₀It is the initial time of system, and has

Based on hypothesis (γ-γ above_LOS) ≈ 0, can approximately there are dR=-Vd τ.Substituting into above formula can solve

At formula (16) left and right two ends while premultiplication Q (t, t₀)

Above formula is rewritable to be

It is reverse to be obtained using step integration rule

Both sides integrate simultaneously

Wherein, γ₀It is initial trajectory inclination angle, γ_LOS0It is the visual line angle of initial bullet.exp(-A₁f₂(t,t₀))=exp (0_2×2)=I_2×2It is the unit matrix of 2 × 2,0 in formula_2×2It is 2 × 2 null matrix.

Q(t,t₀) inverse matrix it is as follows

Φ(t,t₀)=[Q (t, t₀)]^-1=exp (A₁f₂(t,t₀)) (27) wherein, Φ (t, t₀) It is called state-transition matrix.

At formula (26) two ends while premultiplication Φ (t, t₀) can obtain

Matrix A can be obtained by formula (18)₁Proper polynomial be

|λI-A₁|=λ²+(C₁+C₂-1)λ-C₂ (29)

Can be in the hope of matrix A₁Characteristic value be

Wherein,

On the contrary, Guidance Law coefficient C can also be obtained₁And C₂On matrix A₁Two expression formulas of characteristic value

Define f₃(x,t,t₀)=exp (xf₂(t,t₀)), formula (22) is substituted into can be obtained

Due to C₁And C₂It is real number, so matrix A₁Two characteristic values can only be two not wait real number, two equal real numbers With the one kind in three kinds of situations of complex conjugate.Here by matrix A₁Two characteristic values be divided into two kinds of situations and discuss, one kind is λ₁ ≠λ₂, another kind is λ₁=λ₂。

1st, λ is worked as₁≠λ₂When

Can be obtained by spectral factorization formula

Wherein, G₁And G₂It is matrix A₁Spectrum matrix, and have

Wherein, I is unit matrix.Substitution formula (18) and formula (30) are solved

Formula (34) substitution formula (28) right-hand member Section 2 is obtained

Then by formula (28) (34) (37) above, can obtain working as λ₁≠λ₂When trajectory tilt angle γ (t) analytic solutions so that ByThe analytic solutions of Maneuver Acceleration can be obtained, it is as follows

Now prove to work as t → t_fWhen, there is the λ met under certain condition₁And λ₂So thatNeed to A₁Two Eigenvalue λ₁And λ₂For two, grade real number, a pair of conjugate complex numbers, two kinds of situations enter line justification respectively.

a)λ₁And λ₂For two not wait real number

Due to terminal juncture t_fWhen, missile target encounter has R (t_f)=0.Following two conclusions then can be obtained by formula (38)：

(3) if λ₁≤-1,λ₂≤ -1 and λ₁≠λ₂, then as t → t_fWhen, a_LDo not dissipate；

(4) if λ₁＜ -1, λ₂＜ -1 and λ₁≠λ₂, then as t → t_fWhen,

In this case, trajectory shaping Guidance Law of the present invention is imitated with the numerical value of traditional trajectory shaping Guidance Law and nonlinear model True Comparative result analysis is as shown in embodiment one.

b)λ₁And λ₂It is a pair of Conjugate complex roots

If λ₁And λ₂For

Wherein, p and q are real numbers,Above formula is substituted into formula (32) can obtain

Formula (39) (40) substitution formula (38) is obtained into acceleration analytic solutions expression formula is

As can be seen from the above equation, accelerating curve has concussion and occurs, and as t → t_fWhen, concussion frequency tends to infinitely great, Because now R (t_f)=0.But if p ＜ -1, then as t → t_fShi You

Demonstrate,prove, as Conjugate complex roots λ₁And λ₂When meeting real part p ＜ -1, have

In this case, trajectory shaping Guidance Law of the present invention is imitated with the numerical value of traditional trajectory shaping Guidance Law and nonlinear model True Comparative result analysis is as shown in embodiment two.

2nd, λ is worked as₁=λ₂When, i.e. Δ=0

It is easily verified that

Then may determine that matrix A₁Minimum formula m (x) be

M (x)=(x- λ₁)² (45)

Then obtained by broad sense spectrum formula

Wherein

Formula (46) is substituted into the integral term of formula (28) right-hand member, can be obtained

Wherein,

Can then be obtained by formula (28) (46) (48), work as λ₁=λ₂When trajectory tilt angle γ (t) analytic solutions, so as to byThe analytic solutions of Maneuver Acceleration can be obtained, it is as follows

Prove to work as λ now₁=λ₂During ＜ -1,Obvious this problem equivalent is in the following problem of proof：

Work as λ₁During ＜ -1

Prove：

Due to R (t_fλ is worked as in)=0, easily card₁During ＜ -1,

Second limit is oneThe type limit, can be obtained by L ' Hospital methods

Work as λ so as to demonstrate₁=λ₂During ＜ -1,

In this case, trajectory shaping Guidance Law of the present invention is imitated with the numerical value of traditional trajectory shaping Guidance Law and nonlinear model True Comparative result analysis is as shown in embodiment three.

Step 4：Analysis trajectory shaping Guidance Law coefficient C₁And C₂Stable region

Defined function Re (x) is the real part of plural x, then by working above, can obtain drawing a conclusion：

Conclusion 1 works as Re (λ₁)≤- 1 and Re (λ₂)≤- 1, but Re (λ₁) and Re (λ₂) it is different when for -1 when, guidance system is steady Determine, and Maneuver Acceleration instruction does not dissipate.

Conclusion 2 works as Re (λ₁) ＜ -1 and Re (λ₂) ＜ -1, guidance system stabilization, and the final convergence of Maneuver Acceleration instruction To zero.

It should be noted that above inference is still set up when aircraft does variable motion.Therefore, it can first to alignment The characteristic root λ of approximation system₁And λ₂, trajectory shaping Guidance Law coefficient C is then determined by characteristic root₁And C₂.This mode is very big Promoted trajectory shape Guidance Law coefficient span.

1.λ₁And λ₂It is real number

Can be obtained according to inequality group (32), work as λ₁≤-1,λ₂≤ -1 and λ₁And λ₂When being asynchronously -1, C₂＜ -1.Then can be by Inequality group (32) is rewritten as

Formula (32) is substituted into above formula, can be obtained

If C₂Definite value, then C₁On λ₁Derivative be

WhenWhen,C₁With λ₁Increase and increase.

WhenWhen,C₁With λ₁Increase and reduce.

So working asWhen, C₁Minimum value is taken, i.e.,

Work as λ₁=-1 or λ₁=C₂When, C₁Maximum is taken, i.e.,

C₁≤2-2C₂ (56)

In sum, coefficient value scope of trajectory shaping Guidance Law is in the case of this

Formula (57) is substituted into formula (30) can release to draw a conclusion：

If the C of conclusion 3₁=2-2C₂And C₂＜ -1, then λ₁=-1, λ₂=C₂；

If conclusion 4And C₂＜ -1, then

If conclusion 5And C₂＜ -1, then

2.λ₁And λ₂It is a pair of Conjugate complex roots

Re (λ are known by conclusion 1₁)=Re (λ₂) ＜ -1, then can be obtained by formula (30)

Arrangement can be obtained

I.e.

Easily demonstrate,proveTherefore inequality group (60) can be reduced to

MakePerseverance is set up, then must have C₂＜ -1.Thus, convolution (30) obtains drawing a conclusion

If conclusion 6And C₂＜ -1, then λ₁And λ₂It is a pair of Conjugate complex roots and their reality Several are less than -1.

In sum, conclusion 1-6 can merge into following two conclusions：

If the C of conclusion 7₁=2-2C₂And C₂＜ -1, then TSG guidance systems stabilization, and in aircraft close to during target Acceleration instruction does not dissipate.

If the 3-C of conclusion 8₂＜ C₁＜ 2-2C₂And C₂＜ -1, then TSG guidance systems are stable, and in aircraft close to target During acceleration instruction converge on 0.

Traditional trajectory shaping Guidance Law coefficient that early stage Cherry G. is proposed is fixed value, is C₁=6 and C₂=-2, its is right The A for answering₁Characteristic value be λ₁=-1 and λ₂=-2, are then extended the coefficient range (GENEX of Guidance Law by Ohlmeyer E.J. Explicit Guidance is restrained), its coefficient C₁And C₂Value formula is as follows

The corresponding characteristic value of above formula is λ₁=-(n+1) and λ₂=-(n+2).It is of the invention then further expanded by analysis C₁And C₂Span, the shadow region in Fig. 6 is determined by conclusion 7-8, wherein empty circles " o " represent shadow region not Containing (- 1,4) this point.Contrasted with the result (as shown in Fig. 6 solid dots) of GENEX simultaneously.From fig. 6 it can be seen that C₁ And C₂Span has obtained larger extension, and Guidance Law coefficient C₁And C₂Span do not changed by aircraft speed Influence.

Step 5：Terminal velocity control program

The command acceleration a that Section 2 terminal velocity control program wherein described in step 2 is produced_speed, this acceleration work In local level, and perpendicular to speed, it is shown below

Wherein, k_VfIt is positive parameter undetermined, g₀It is the acceleration of gravity of sea level altitude, R₀It is start time aircraft and mesh Remaining flying distance between mark, ψ₀It is the course angle of initial time, ψ_LOS0It is the sight line azimuth of initial time.Sgn (x) is Sign function, it is as follows

By adjusting k_VfValue, aircraft crossrange maneuvering amplitude can be controlled, so that when adjusting aircraft hit Speed.a_speedWith a_TSGThe latter end of aircraft flight is acted on simultaneously.a_speedSize with aircraft close to target gradually Weaken, can so avoid a in aircraft soon hit_speedTo a_TSGInterference.

Aircraft predicted the speed V of hit by Ballistic Simulation of Underwater before the terminal guidance stage is entered_f(K), then press Secant methods correct k_VfValue, finally obtain k_VfIdeal value, it is as follows

Wherein, V_fcmdIt is desired terminal velocity, K is iterations.

Specific embodiment

Embodiment one

In the present embodiment, the influence of gravity is ignored, aircraft does deceleration downglide motion：V (t)=1000-5t (m/s), And wish aircraft with the trajectory tilt angle hit of 0deg.The primary condition of aircraft is：x₀=0km, H₀=10km, γ₀= 0deg, the primary condition of target is：x_T=50km, H_T=0km.

By assuming γ-γ_LOS≈ 0,Solve expression formula of the remaining flying distance on the time By solving equation R (t_f)=0 can obtain terminal juncture t_fValue.Characteristic root is set Value be respectively λ₁=-2, λ₂=-2.5, then formula (32) Guidance Law coefficient C can be solved₁=10.5, C₂=-5.Now, this can be obtained Invention trajectory shaping Guidance Law analytical form be

Traditional trajectory as contrast shapes Guidance Law, and its speed is for definite value and the short transverse equation of motion is linear, shape Formula is shown below

Wherein, t₀It is initial time, t_f2It is the final moment estimated, H₀It is elemental height,Be initial velocity height Component on direction,It is terminal juncture desired speed component in the height direction.

By contrast (66) and formula (67), it can be seen that H₀=-R (t₀)γ_LOS0, If ignoring the change that trajectory of the present invention shapes Guidance Law medium velocity size, its reduced form is exactly traditional trajectory shaping guidance Rule.Fig. 3 (a) illustrates the ballistic curve under trajectory of the present invention shaping Guidance Law effect, Fig. 3 (b) illustrate trajectory of the present invention into Shape Guidance Law and traditional trajectory shaping Guidance Law, the acceleration time graph of numerical simulation contrast.Result shows, trajectory of the present invention Shaping Guidance Law can converge to 0 in hit season aircraft acceleration, and because the present invention does not ignore velocity variations Influence, accelerating curve has smaller error with Numerical Simulation Results.

Embodiment two

The present embodiment is consistent with the primary condition of embodiment one, and setting characteristic value is λ₁=-2+2i, λ₂=-2-2i.Can obtain Guidance Law coefficient is C₁=13, C₂=-8, then the analytical form of trajectory of the present invention shaping Guidance Law be

The analytical form of traditional trajectory shaping Guidance Law is

Fig. 4 (a) illustrates the ballistic curve under trajectory shaping Guidance Law effect of the present invention, and Fig. 4 (b) illustrates bullet of the present invention Road shapes Guidance Law with traditional trajectory shaping Guidance Law, the acceleration time graph of numerical simulation contrast.Result shows, of the invention Trajectory shapes Guidance Law and nonlinear model the approximation of height, it was demonstrated that acceleration instruction has vibration really, but finally When close to target, acceleration can tend to 0.

Embodiment three

The present embodiment is consistent with the primary condition of embodiment one, and setting characteristic value is λ₁=λ₂=-2, then can obtain Guidance Law system Number is C₁=9, C₂=-4, then the analytical form of trajectory of the present invention shaping Guidance Law be

The analytical form of traditional trajectory shaping Guidance Law is

Fig. 5 (a) illustrates the ballistic curve under trajectory shaping Guidance Law effect of the present invention, and Fig. 5 (b) illustrates system of the present invention Rule is led with traditional trajectory shaping Guidance Law, the acceleration time graph of numerical simulation contrast.Result shows, in this case this Invention trajectory shaping Guidance Law can highly approach the Numerical Simulation Results of nonlinear model, can control aircraft in hit mesh Target moment, trajectory tilt angle, acceleration instruction meet constraints.

Example IV

The present invention is using CAV-H aircraft as simulation model.Vehicle mass is 907.2kg, and pneumatic area of reference is 0.48387m², flying drilling angle is between [- 25,25] degree.Aerodynamic data is calculated using following fitting formula

Wherein C_LIt is lift coefficient, C_DIt is resistance coefficient, α is Aircraft Angle of Attack, unit：Rad, C_L0It is zero-incidence lift system Number,It is lift coefficient slope, C_D0It is zero-lift drag coefficient, K is induced drag coefficient, C_L、C_D0It is Mach number with K Shown in function, such as Fig. 7 (a), (b).

Aircraft using Bank to Turn (BTT) control mode, the angle of attack and angle of heel, and can used as flight controlled quentity controlled variable Calculated by command acceleration overload：

Wherein x₁It is unit vector, direction is perpendicular to the vertical guide for including velocity vector；x₂It is unit vector, direction is hung down It is straight in velocity vector and positioned at urging in vertical guide.

x₁=[sin (ψ) ,-cos (ψ), 0]^T (75)

x₂=[- sin (γ) cos (ψ) ,-sin (γ) sin (ψ), cos (γ)]^T (76)

The present embodiment requirement aircraft in end of time with desired trajectory tilt angle and speed hit, and cause machine Dynamic overload is 0.The simulation parameter of aircraft sets as follows：x₀=y₀=0, H₀=20km, V_o=2500m/s, γ₀=0deg and ψ₀ =0deg, end point requirements are：x_f=50km, y_f=H_f=0 and γ_f=-80deg, respectively consider terminal velocity be 1400m/s and Two kinds of situations of 1200m/s, in addition, also require the terminal angle of attack in the range of ± 2deg, here by making the motor-driven overload of terminal Level off to zero come meet the terminal angle of attack constraint.Fig. 8 is by largely emulating the terminal velocity V for obtaining_fWith parameter k_vfChange it is bent Line, it is clear that k_vfBigger, trajectory gets over bending, and terminal velocity is smaller.As can be seen from Figure 8, V_fAnd k_vfThere is a kind of approximately linear Relation, this is conducive to being quickly found out the k for meeting terminal velocity requirement using secant methods_vfValue, as long as typically entering to emulate for several times Iteration, it is possible to obtain solution high-precision enough.By calculating, k_vfMeet V_f=1400m/s ideal values are 21.4088, are met V_fThe ideal value of=1200m/s is 54.2326.

Fig. 9 is directed to the simulation result that above two difference terminal velocity situation is obtained.Wherein, Fig. 9 (a) ballistic curves Figure, 9 (b) rate curve, Fig. 9 (c) trajectory tilt angle curves, the motor-driven overload curves of Fig. 9 (d), Fig. 9 (e) rolling angular curves, Fig. 9 (f) To flow line of buckling.It can be seen that under explicit guidance's rule of the present invention effect, the speed and trajectory tilt angle of drop point are accorded with Close and require, the motor-driven overload of terminal is also close to 0.For the less situation of terminal velocity, aircraft needs curved to disappear around bigger Energy consumption, meanwhile, maximum is pressed with very big reduction to flow.From rolling angular curve, aircraft is finally with inverted flight state Hit.

Claims (1)

1. a kind of explicit Guidance rule with terminal velocity, trajectory tilt angle and overload constraint, is characterised by：Specifically include following step Suddenly：

Step 1：Set up kinematics and dynamics equation

The inertial reference system o-xyH connected firmly with ground is set up, it is longitudinal range to regard aircraft as particle M, x, and y is laterally to penetrate Journey, H is height above sea level；V is the speed of aircraft；γ is the trajectory tilt angle of aircraft, and relative level reference line is rotated counterclockwise For just；ψ is course angle, i.e. the horizontal component of aircraft speed and the angle of x-axis, is rotated counterclockwise as just relative to x-axis；Do not examine Consider the influence of earth curvature and rotation, its kinematics and dynamics equation group is as follows：

$\frac{dx}{dt}=Vcos\left(\mathrm{\γ}\right)cos\left(\mathrm{\ψ}\right)---\left(1\right)$

$\frac{dy}{dt}=Vcos\left(\mathrm{\γ}\right)sin\left(\mathrm{\ψ}\right)---\left(2\right)$

$\frac{dH}{dt}=Vsin\left(\mathrm{\γ}\right)---\left(3\right)$

$\frac{dV}{dt}=-\frac{D}{m}-gsin\left(\mathrm{\γ}\right)---\left(4\right)$

$\frac{d\mathrm{\γ}}{dt}=\frac{Lcos\left(\mathrm{\σ}\right)}{mV}-\frac{gcos\left(\mathrm{\γ}\right)}{V}---\left(5\right)$

$\frac{d\mathrm{\ψ}}{dt}=-\frac{Lsin\left(\mathrm{\σ}\right)}{mVcos\left(\mathrm{\γ}\right)}---\left(6\right)$

Wherein, t is the aircraft flight time, and m is vehicle mass, and σ is roll angle, and L is lift, and D is resistance, and g adds for gravity Speed；

Step 2：Explicit Guidance rule summary

Under the effect of explicit Guidance rule, guidance system constitutes the command acceleration of generation by two：Section 1 is trajectory shaping The command acceleration a that Guidance Law is produced_TSG, this accelerated energy directing aircraft is from predetermined direction hit；Section 2 is eventually The command acceleration a that spot speed control program is produced_speed, this accelerated energy control aircraft does crossrange maneuvering, so as to adjust life Velocity magnitude during middle target；The command acceleration expression formula that guidance system is produced is as follows

a_cmd=a_TSG+a_speed (7)

Step 3：Solve the analytic solutions that trajectory shapes Guidance Law

The command acceleration a that Section 1 trajectory shaping Guidance Law is produced in the step 2_TSGIt is made up of three events, is used to lead respectively Draw aircraft target, control velocity attitude during aircraft hit, balancing gravity acceleration is perpendicular to speed side To component；Acceleration a_TSGDirection perpendicular to aircraft current velocity vector, its expression formula is shown in formula (8)

$a_{TSG}=\frac{V^2}{R}\{C_1\[\hat{r}-(\hat{r}\·\hat{v})\hat{v}\]+C_2\[{\hat{v}}_f-({\hat{v}}_f\·\hat{v})\hat{v}\]\}-g_n---\left(8\right)$

Wherein, R is remaining flying distance；C₁And C₂It is Guidance Law coefficient；It is unit direction vector of the aircraft to line of sight；It is the unit direction vector of aircraft speed；Predetermined speed direction vector when being aircraft hit, is counted by formula (9) Calculate；g_nIt is component of the acceleration of gravity perpendicular to speed；

${\hat{v}}_f=\left[\begin{array}{c}\cos \left({\mathrm{\γ}}_f\right)\cos \left({\mathrm{\ψ}}_f\right)\\ \cos \left({\mathrm{\γ}}_f\right)\sin \left({\mathrm{\ψ}}_f\right)\\ \sin \left({\mathrm{\γ}}_f\right)\end{array}\right]---\left(9\right)$

Wherein, γ_fThe predetermined trajectory tilt angle of aircraft, is typically in the range of between -70 ° and -90 ° during for hit, is constant value；ψ_f It is as follows equal to the azimuth of current flight device to line of sight

${\mathrm{\ψ}}_f=arctan\left(\frac{\hat{r}{|}_y}{\hat{r}{|}_x}\right)---\left(10\right)$

Wherein,WithVector is represented respectivelyIn the component in reference axis x, y direction；

Ignore the influence of gravity, then in fore-and-aft plane, the expression formula of aircraft guidance rule can be by geometric vector formal argument Into trigonometric function form：

$V\frac{d\mathrm{\γ}}{dt}=\frac{V^2}{R}(-C_{1}sin(\mathrm{\γ}-{\mathrm{\γ}}_{LOS})-C_{2}sin(\mathrm{\γ}-{\mathrm{\γ}}_f\left)\right)---\left(11\right)$

Wherein, γ_LOSIt is the visual line angle of bullet, relative level reference line is rotated counterclockwise as just；

Sight line angular rate of change is：

$\frac{{\mathrm{d\γ}}_{LOS}}{dt}=-\frac{Vsin(\mathrm{\γ}-{\mathrm{\γ}}_{LOS})}{R}---\left(12\right)$

Assuming that γ ≈ γ_LOS≈γ_f, then can be to the trigonometric function line in formula (11) (12), and arrangement obtains following linear time-varying System

$\left[\begin{array}{c}\stackrel{\·}{\mathrm{\γ}}\\ {\stackrel{\·}{\mathrm{\γ}}}_{LOS}\end{array}\right]=f_1\left(t\right)A_1\left[\begin{array}{c}\mathrm{\γ}\\ {\mathrm{\γ}}_{LOS}\end{array}\right]+f_1\left(t\right)B_1{\mathrm{\γ}}_f---\left(13\right)$

Wherein

$f_1\left(t\right)=\frac{V\left(t\right)}{R\left(t\right)}---\left(14\right)$

$A_1=\left[\begin{array}{cc}-(C_1+C_2)& C_1\\ -1& 1\end{array}\right]---\left(15\right)$

B₁=[C₂ 0]^T (16)

Wherein γ and γ_LOSIt is state variable, γ_fIt is control variables, the transposition of subscript " T " representing matrix here and in full text；

The solution of linear time varying system is carried out using a new method based on spectral factorization；

Definition

$Q(t,t_0)=\exp (-\underset{t_0}{\overset{t}{\∫}}{A}_1{f}_1(\mathrm{\τ}\left)d\mathrm{\τ}\right)=\exp (-A_1{f}_2(t,t_0\left)\right)---\left(17\right)$

Wherein t₀It is the initial time of system, and has

$f_2(t,t_0)=\underset{t_0}{\overset{t}{\∫}}{f}_1\left(\mathrm{\τ}\right)d\mathrm{\τ}=\underset{t_0}{\overset{t}{\∫}}\frac{V\left(\mathrm{\τ}\right)}{R\left(\mathrm{\τ}\right)}d\mathrm{\τ}---\left(18\right)$

Based on the assumption that (γ-γ_LOS) ≈ 0, can approximately there are dR=-Vd τ；Substituting into above formula can solve

$f_2(t,t_0)=-\underset{R\left(t_0\right)}{\overset{R\left(t\right)}{\∫}}\frac{1}{R}dR=\ln \left(\frac{R\left(t_0\right)}{R\left(t\right)}\right)---\left(19\right)$

At formula (13) left and right two ends while premultiplication Q (t, t₀)

$\exp (-A_1{f}_2(t,t_0\left)\right)\left[\begin{array}{c}\stackrel{\·}{\mathrm{\γ}}\\ {\stackrel{\·}{\mathrm{\γ}}}_{LOS}\end{array}\right]-\exp (-A_1{f}_2(t,t_0\left)\right)f_1\left(t\right)A_1\left[\begin{array}{c}\mathrm{\γ}\\ {\mathrm{\γ}}_{LOS}\end{array}\right]=\exp (-A_1{f}_2(t,t_0\left)\right)f_1\left(t\right)B_1{\mathrm{\γ}}_f---\left(20\right)$

Above formula is rewritable to be

$\exp (-A_1{f}_2(t,t_0\left)\right)\frac{d}{dt}\left[\begin{array}{c}\mathrm{\γ}\\ {\mathrm{\γ}}_{LOS}\end{array}\right]+\frac{d}{dt}\[\exp (-A_1{f}_2(t,t_0\left)\right)\]\left[\begin{array}{c}\mathrm{\γ}\\ {\mathrm{\γ}}_{LOS}\end{array}\right]=\exp (-A_1{f}_2(t,t_0\left)\right)f_1\left(t\right)B_1{\mathrm{\γ}}_f---\left(21\right)$

It is reverse to be obtained using step integration rule

$\frac{d}{dt}\left\{\exp (-A_1{f}_2(t,t_0\left)\right)\left[\begin{array}{c}\mathrm{\γ}\\ {\mathrm{\γ}}_{LOS}\end{array}\right]\right\}=\exp (-A_1{f}_2(t,t_0\left)\right)f_1\left(t\right)B_1{\mathrm{\γ}}_f---\left(22\right)$

Both sides integrate simultaneously

$\exp (-A_1{f}_2(t,t_0\left)\right)\left[\begin{array}{c}\mathrm{\γ}\\ {\mathrm{\γ}}_{LOS}\end{array}\right]-\exp (-A_1{f}_2(t,t_0\left)\right)\left[\begin{array}{c}{\mathrm{\γ}}_0\\ {\mathrm{\γ}}_{LOS0}\end{array}\right]=\underset{t_0}{\overset{t}{\∫}}\exp (-A_1{f}_2(t,t_0\left)\right)f_1\left(\mathrm{\τ}\right)B_1{\mathrm{\γ}}_{f}d\mathrm{\τ}---\left(23\right)$

Wherein, γ₀It is initial trajectory inclination angle, γ_LOS0It is the visual line angle of initial bullet；exp(-A₁f₂(t,t₀))=exp (0_2×2)= I_2×2It is the unit matrix of 2 × 2,0 in formula_2×2It is 2 × 2 null matrix；

Q(t,t₀) inverse matrix it is as follows

Φ(t,t₀)=[Q (t, t₀)]-¹=exp (A₁f₂(t,t₀)) (24)

Wherein, Φ (t, t₀) it is called state-transition matrix；

At formula (23) two ends while premultiplication Φ (t, t₀) can obtain

$\left[\begin{array}{c}\mathrm{\γ}\left(t\right)\\ {\mathrm{\γ}}_{LOS}\left(t\right)\end{array}\right]=\mathrm{\Φ}(t,t_0)\left[\begin{array}{c}{\mathrm{\γ}}_0\\ {\mathrm{\γ}}_{LOS0}\end{array}\right]+\underset{t_0}{\overset{t}{\∫}}{\mathrm{\γ}}_f{f}_1\left(\mathrm{\τ}\right)\mathrm{\Φ}(t,\mathrm{\τ})B_{1}d\mathrm{\τ}---\left(25\right)$

Matrix A can be obtained by formula (15)₁Proper polynomial be

|λI-A₁|=λ²+(C₁+C₂-1)λ-C₂ (26)

Can be in the hope of matrix A₁Characteristic value be

$\left\{\begin{array}{c}{\mathrm{\λ}}_1=\frac{-(C_1+C_2-1)+\sqrt{\mathrm{\Δ}}}{2}\\ {\mathrm{\λ}}_2=\frac{-(C_1+C_2-1)-\sqrt{\mathrm{\Δ}}}{2}\end{array}\right.---\left(27\right)$

Wherein,

$\mathrm{\Δ}=C_1^2+C_2^2+1+2C_1{C}_2-2C_1+2C_2---\left(28\right)$

On the contrary, Guidance Law coefficient C can also be obtained₁And C₂On matrix A₁Two expression formulas of characteristic value

$\left\{\begin{array}{c}{C}_1=1-{\mathrm{\λ}}_1-{\mathrm{\λ}}_2+{\mathrm{\λ}}_1{\mathrm{\λ}}_2\\ C_2=-{\mathrm{\λ}}_1{\mathrm{\λ}}_2\end{array}\right.---\left(29\right)$

Define f₃(x,t,t₀)=exp (xf₂(t,t₀)), formula (19) is substituted into can be obtained

$f_3(x,t,t_0)=\exp \left({\mathrm{xf}}_2\right(t,t_0\left)\right)={\left(\frac{R\left(t_0\right)}{R\left(t\right)}\right)}^x---\left(30\right)$

Due to C₁And C₂It is real number, so matrix A₁Two characteristic values can only be two and do not wait real number, two equal real numbers and multiple One kind in three kinds of situations of conjugation；Here by matrix A₁Two characteristic values be divided into two kinds of situations, one kind is λ₁≠λ₂, it is another It is λ₁=λ₂；

(I) λ, is worked as₁≠λ₂When

Can be obtained by spectral factorization formula

$\begin{array}{c}\mathrm{\Φ}(t,t_0)=f_3({\mathrm{\λ}}_1,t,t_0)G_1+f_3({\mathrm{\λ}}_2,t,t_0)G_2\\ ={\left(\frac{R\left(t_0\right)}{R\left(t\right)}\right)}^{{\mathrm{\λ}}_1}{G}_1+{\left(\frac{R\left(t_0\right)}{R\left(t\right)}\right)}^{{\mathrm{\λ}}_2}{G}_2\end{array}---\left(31\right)$

Wherein, G₁And G₂It is matrix A₁Spectrum matrix, and have

$\left\{\begin{array}{c}{G}_1=\frac{A_1-{\mathrm{\λ}}_{2}I}{{\mathrm{\λ}}_1-{\mathrm{\λ}}_2}\\ G_2=\frac{A_1-{\mathrm{\λ}}_{1}I}{{\mathrm{\λ}}_2-{\mathrm{\λ}}_1}\end{array}\right.---\left(32\right)$

Wherein, I is unit matrix；Substitution formula (15) and formula (27) are solved

$\left\{\begin{array}{c}{G}_1=\frac{1}{\sqrt{\mathrm{\Δ}}}\left[\begin{array}{cc}-\frac{C_1+C_2+1-\sqrt{\mathrm{\Δ}}}{2}& C_1\\ -1& \frac{C_1+C_2+1+\sqrt{\mathrm{\Δ}}}{2}\end{array}\right]\\ G_2=\frac{1}{\sqrt{\mathrm{\Δ}}}\left[\begin{array}{cc}\frac{C_1+C_2+1+\sqrt{\mathrm{\Δ}}}{2}& -C_1\\ 1& -\frac{C_1+C_2+1-\sqrt{\mathrm{\Δ}}}{2}\end{array}\right]\end{array}\right.---\left(33\right)$

Formula (31) substitution formula (25) right-hand member Section 2 is obtained

$\begin{array}{c}\underset{t_0}{\overset{t}{\∫}}{\mathrm{\γ}}_f{f}_1\left(\mathrm{\τ}\right)\mathrm{\Φ}(t,\mathrm{\τ})B_{1}d\mathrm{\τ}\\ =\underset{t_0}{\overset{t}{\∫}}{\mathrm{\γ}}_f\[{\left(\frac{R\left(\mathrm{\τ}\right)}{R\left(t\right)}\right)}^{{\mathrm{\λ}}_1}{G}_1+{\left(\frac{R\left(\mathrm{\τ}\right)}{R\left(t\right)}\right)}^{{\mathrm{\λ}}_2}{G}_2\]\frac{1}{R\left(\mathrm{\τ}\right)}\left[\begin{array}{c}{C}_2\\ 0\end{array}\right]\[V\left(\mathrm{\τ}\right)d\mathrm{\τ}\]\\ =\underset{R\left(t_0\right)}{\overset{R\left(t\right)}{\∫}}{\mathrm{\γ}}_f\[\frac{R^{{\mathrm{\λ}}_1-1}\left(\mathrm{\τ}\right)}{R^{{\mathrm{\λ}}_1}\left(t\right)}{G}_1+\frac{R^{{\mathrm{\λ}}_2-1}\left(\mathrm{\τ}\right)}{R^{{\mathrm{\λ}}_2}\left(t\right)}{G}_2\]\left[\begin{array}{c}{C}_2\\ 0\end{array}\right]\[-dR\left(\mathrm{\τ}\right)\]\\ =\[\frac{{\mathrm{\γ}}_f}{{\mathrm{\λ}}_1}(\frac{R^{{\mathrm{\λ}}_1}\left(t_0\right)}{R^{{\mathrm{\λ}}_1}\left(t\right)}-1)G_1+\frac{{\mathrm{\γ}}_f}{{\mathrm{\λ}}_2}(\frac{R^{{\mathrm{\λ}}_2}\left(t_0\right)}{R^{{\mathrm{\λ}}_2}\left(t\right)}-1)G_2\]\left[\begin{array}{c}{C}_2\\ 0\end{array}\right]\end{array}---\left(34\right)$

Then by formula (25) (31) (34) above, can obtain working as λ₁≠λ₂When trajectory tilt angle γ (t) analytic solutions, so as to byThe analytic solutions of Maneuver Acceleration can be obtained, it is as follows

$\begin{array}{c}{a}_L\left(t\right)=\frac{{\mathrm{\λ}}_1{V}^2\left(t\right)}{\sqrt{\mathrm{\Δ}}}\frac{R^{-{\mathrm{\λ}}_1-1}\left(t\right)}{R^{-{\mathrm{\λ}}_1}\left(t_0\right)}\[C_1{\mathrm{\γ}}_{LOS0}-\frac{C_1+C_2+1-\sqrt{\mathrm{\Δ}}}{2}({\mathrm{\γ}}_0+\frac{C_2{\mathrm{\γ}}_f}{{\mathrm{\λ}}_1})\]\\ +\frac{{\mathrm{\λ}}_1{V}^2\left(t\right)}{\sqrt{\mathrm{\Δ}}}\frac{R^{-{\mathrm{\λ}}_2-1}\left(t\right)}{R^{-{\mathrm{\λ}}_2}\left(t_0\right)}\[-C_1{\mathrm{\γ}}_{LOS0}+\frac{C_1+C_2+1+\sqrt{\mathrm{\Δ}}}{2}({\mathrm{\γ}}_0+\frac{C_2{\mathrm{\γ}}_f}{{\mathrm{\λ}}_2})\]\end{array}---\left(35\right)$

Now prove to work as t → t_fWhen, there is the λ met under certain condition₁And λ₂So thatNeed to A₁Two features Value λ₁And λ₂For two, grade real number, a pair of conjugate complex numbers, two kinds of situations enter line justification respectively；

a)λ₁And λ₂For two not wait real number

Due to terminal juncture t_fWhen, missile target encounter has R (t_f)=0；Following two conclusions then can be obtained by formula (35)：

(1) if λ₁≤-1,λ₂≤ -1 and λ₁≠λ₂, then as t → t_fWhen, a_LDo not dissipate；

(2) if λ₁＜ -1, λ₂＜ -1 and λ₁≠λ₂, then as t → t_fWhen,

b)λ₁And λ₂It is a pair of Conjugate complex roots

If λ₁And λ₂For

$\left\{\begin{array}{c}{\mathrm{\λ}}_1=p+iq\\ {\mathrm{\λ}}_2=p-iq\end{array}\right.---\left(36\right)$

Wherein, p and q are real numbers,Above formula is substituted into formula (29) can obtain

$\left\{\begin{array}{c}{C}_1=1-2p+p^2+q^2\\ C_2=-p^2-q^2\end{array}\right.---\left(37\right)$

Formula (36) (37) substitution formula (35) is obtained into acceleration analytic solutions expression formula is

$\begin{array}{c}{a}_L\left(t\right)=\frac{V^2\left(t\right)R^{-p-1}\left(t\right)}{R^{-p}\left(t_0\right)}\sin \[q\ln \left(\frac{R\left(t\right)}{R\left(t_0\right)}\right)\]\frac{\[-{\mathrm{pC}}_1{\mathrm{\γ}}_{LOS0}+(p-p^2+q^2){\mathrm{\γ}}_0+(1-p)C_2{\mathrm{\γ}}_f}{q})\]\\ +\frac{V^2\left(t\right)R^{-p-1}\left(t\right)}{R^{-p}\left(t_0\right)}\cos \[q\ln \left(\frac{R\left(t\right)}{R\left(t_0\right)}\right)\]\[C_1{\mathrm{\γ}}_{LOS0}-(1-2p){\mathrm{\γ}}_0+C_2{\mathrm{\γ}}_f\]\end{array}---\left(38\right)$

As can be seen from the above equation, accelerating curve has concussion and occurs, and as t → t_fWhen, concussion frequency tends to infinitely great, because Now R (t_f)=0；But if p ＜ -1, then as t → t_fShi You

$\underset{t\→t_f}{\lim }|R^{-p-1}\left(t\right)sin\[qln\left(\frac{R\left(t\right)}{R\left(t_0\right)}\right)\]|\≤\underset{t\→t_f}{\lim }\left|R^{-p-1}\left(t\right)\right|=0---\left(39\right)$

$\underset{t\→t_f}{\lim }|R^{-p-1}\left(t\right)\cos \[qln\left(\frac{R\left(t\right)}{R\left(t_0\right)}\right)\]|\≤\underset{t\→t_f}{\lim }\left|R^{-p-1}\left(t\right)\right|=0---\left(40\right)$

Demonstrate,prove, as Conjugate complex roots λ₁And λ₂When meeting real part p ＜ -1, have

(II) λ, is worked as₁=λ₂When, i.e. Δ=0

It is easily verified that

$\left\{\begin{array}{c}{A}_1-{\mathrm{\λ}}_{1}I\≠0\\ {(A_1-{\mathrm{\λ}}_{1}I)}^2=0\end{array}\right.---\left(41\right)$

Then may determine that matrix A₁Minimum formula m (x) be

M (x)=(x- λ₁)² (42)

Then obtained by broad sense spectrum formula

$\begin{array}{c}\mathrm{\Φ}(t,t_0)=f_3({\mathrm{\λ}}_1,t,t_0)I+\frac{\∂f_3(x,t,t_0)}{\∂x}{|}_{x={\mathrm{\λ}}_1}(A_1-{\mathrm{\λ}}_{1}I)\\ ={\left(\frac{R\left(t_0\right)}{R\left(t\right)}\right)}^{{\mathrm{\λ}}_1}I+{\left(\frac{R\left(t_0\right)}{R\left(t\right)}\right)}^{{\mathrm{\λ}}_1}\ln \left(\frac{R\left(t_0\right)}{R\left(t\right)}\right)A_2\end{array}---\left(43\right)$

Wherein

$A_2=(A_1-{\mathrm{\λ}}_{1}I)=\left[\begin{array}{cc}-\frac{C_1+C_2+1}{2}& C_1\\ -1& \frac{C_1+C_2+1}{2}\end{array}\right]---\left(44\right)$

Formula (43) is substituted into the integral term of formula (25) right-hand member, can be obtained

$\begin{array}{c}\underset{t_0}{\overset{t}{\∫}}{\mathrm{\γ}}_f{f}_1\left(\mathrm{\τ}\right)\mathrm{\Φ}(t,\mathrm{\τ})B_{1}d\mathrm{\τ}\\ =\underset{t_0}{\overset{t}{\∫}}{\mathrm{\γ}}_f\[{\left(\frac{R\left(\mathrm{\τ}\right)}{R\left(t\right)}\right)}^{{\mathrm{\λ}}_1}I+{\left(\frac{R\left(\mathrm{\τ}\right)}{R\left(t\right)}\right)}^{{\mathrm{\λ}}_1}\ln \left(\frac{R\left(\mathrm{\τ}\right)}{R\left(t\right)}\right)A_2\]\frac{1}{R\left(\mathrm{\τ}\right)}\left[\begin{array}{c}{C}_2\\ 0\end{array}\right]\[V\left(\mathrm{\τ}\right)d\mathrm{\τ}\]\\ =\underset{R\left(t_0\right)}{\overset{R\left(t\right)}{\∫}}{\mathrm{\γ}}_f\[\frac{R^{{\mathrm{\λ}}_1-1}\left(\mathrm{\τ}\right)}{R^{{\mathrm{\λ}}_1}\left(t\right)}I+\frac{R^{{\mathrm{\λ}}_1-1}\left(\mathrm{\τ}\right)}{R^{{\mathrm{\λ}}_1}\left(t\right)}\ln \left(\frac{R\left(\mathrm{\τ}\right)}{R\left(t\right)}\right)A_2\]\left[\begin{array}{c}{C}_2\\ 0\end{array}\right]\[-dR\left(\mathrm{\τ}\right)\]\\ =\frac{{\mathrm{\γ}}_f}{{\mathrm{\λ}}_1^2}\[1-{\left(\frac{R\left(t_0\right)}{R\left(t\right)}\right)}^{{\mathrm{\λ}}_1}\]A_3\left[\begin{array}{c}{C}_2\\ 0\end{array}\right]+\frac{{\mathrm{\γ}}_f}{{\mathrm{\λ}}_1}{\left(\frac{R\left(t_0\right)}{R\left(t\right)}\right)}^{{\mathrm{\λ}}_1}\ln \left(\frac{R\left(t_0\right)}{R\left(t\right)}\right)A_2\left[\begin{array}{c}{C}_2\\ 0\end{array}\right]\end{array}---\left(45\right)$

Wherein,

$A_3=A_2-{\mathrm{\λ}}_{1}I=\left[\begin{array}{cc}-1& C_1\\ -1& C_1+C_2\end{array}\right]---\left(46\right)$

Can then be obtained by formula (25) (43) (45), work as λ₁=λ₂When trajectory tilt angle γ (t) analytic solutions, so as to byThe analytic solutions of Maneuver Acceleration can be obtained, it is as follows

$\begin{array}{c}{a}_L\left(t\right)=V^2\left(t\right)\frac{R^{-{\mathrm{\λ}}_1-1}\left(t\right)}{R^{-{\mathrm{\λ}}_1}\left(t_0\right)}\[-(C_1+C_2){\mathrm{\γ}}_0+C_2{\mathrm{\γ}}_f+C_1{\mathrm{\γ}}_{LOS0}\]\\ -V^2\left(t\right)\frac{R^{-{\mathrm{\λ}}_1-1}\left(t\right)}{R^{-{\mathrm{\λ}}_1}\left(t_0\right)}\frac{C_1+C_2-1}{2}\ln \left(\frac{R\left(t_0\right)}{R\left(t\right)}\right)(C_1{\mathrm{\γ}}_{LOS0}\\ -\frac{C_1+C_2+1}{2}{\mathrm{\γ}}_0\frac{C_1+C_2+1}{C_1+C_2-1}{C}_2{\mathrm{\γ}}_f)\end{array}---\left(47\right)$

Prove to work as λ now₁=λ₂During ＜ -1,Obvious this problem equivalent is in the following problem of proof：

Work as λ₁During ＜ -1

$\left\{\begin{array}{c}\underset{t\→t_f}{\lim }{R}^{-{\mathrm{\λ}}_1-1}\left(t\right)=0\\ \underset{t\→t_f}{\lim }\left(R^{-{\mathrm{\λ}}_1-1}\left(t\right)ln\left(\frac{R\left(t_0\right)}{R\left(t\right)}\right)\right)=0\end{array}\right.$

Prove：

Due to R (t_fλ is worked as in)=0, easily card₁During ＜ -1,

Second limit is oneThe type limit, can be obtained by L ' Hospital methods

$\underset{t\→t_f}{\lim }\left(R^{-{\mathrm{\λ}}_1-1}\right(t\left)ln\right(\frac{R\left(t_0\right)}{R\left(t\right)}\left)\right)=\underset{t\→t_f}{\lim }\frac{ln\left(\frac{R\left(t_0\right)}{R\left(t\right)}\right)}{\frac{1}{R^{-{\mathrm{\λ}}_1-1}\left(t\right)}}=\underset{t\→t_f}{\lim }\frac{{\[ln\left(\frac{R\left(t_0\right)}{R\left(t\right)}\right)\]}^{\′}}{{\[\frac{1}{R^{-{\mathrm{\λ}}_1-1}\left(t\right)}\]}^{\′}}=-\underset{t\→t_f}{\lim }\frac{R^{-{\mathrm{\λ}}_1-1}\left(t\right)}{({\mathrm{\λ}}_1+1)}=0---\left(48\right)$

Work as λ so as to demonstrate₁=λ₂During ＜ -1,

Step 4：Analysis trajectory shaping Guidance Law coefficient C₁And C₂Stable region

It is given to make the characteristic root λ of the line approximation system of guidance system stabilization₁And λ₂, then determine that trajectory shapes by characteristic root Guidance Law coefficient C₁And C₂；To characteristic root λ₁And λ₂Discuss in two kinds of situation；

Ⅰ.λ₁And λ₂It is real number

Analyzed by step 3, work as λ₁≤-1,λ₂≤ -1 and λ₁And λ₂When being asynchronously -1, guidance system stabilization, and motor-driven acceleration Degree instruction does not dissipate；Thus formula (29) can be rewritten as following form

$\left\{\begin{array}{c}{C}_1=1-{\mathrm{\&lambda;}}_1+C_2/{\mathrm{\&lambda;}}_1-C_2\\ C_2\&le;{\mathrm{\&lambda;}}_1\&le;-1\\ C_2<-1\end{array}\right.---\left(49\right)$

If C₂Definite value, then C₁On λ₁Derivative be

$\frac{{\mathrm{dC}}_1}{{\mathrm{d\&lambda;}}_1}=-(1+\frac{C_2}{{\mathrm{\&lambda;}}_1^2})---\left(50\right)$

It is easy to get, whenWhen, C₁Minimum value is taken, i.e.,

$C_1\&GreaterEqual;1+2\sqrt{-C_2}-C_2---\left(51\right)$

Work as λ₁=-1 or λ₁=C₂When, C₁Maximum is taken, i.e.,

C₁≤2-2C₂ (52)

To sum up, λ₁And λ₂For the coefficient value scope of trajectory shaping Guidance Law in the case of real number is

$\left\{\begin{array}{c}1+2\sqrt{-C_2}-C_2\&le;C_1\&le;2-2C_2\\ C_2<-1\end{array}\right.---\left(53\right)$

Ⅱ.λ₁And λ₂It is a pair of Conjugate complex roots

Defined function Re (x) is the real part of plural x, and being analyzed by step 3 to obtain, as Re (λ₁)=Re (λ₂) ＜ -1 when, guidance system System stabilization, and Maneuver Acceleration instruction finally converges to zero；Can then be obtained by formula (27)

$\left\{\begin{array}{c}{C}_1+C_2-1>2\\ \mathrm{\&Delta;}=C_1^2+C_2^2+1+2C_1{C}_2-2C_1+2C_2<0\end{array}\right.---\left(54\right)$

Arrangement can be obtained

$\left\{\begin{array}{c}{C}_1>3-C_2\\ {(C_1+C_2-1)}^2<-4C_2\end{array}\right.---\left(55\right)$

I.e.

$\left\{\begin{array}{c}{C}_1>3-C_2\\ 1-2\sqrt{-C_2}-C_2<C_1<1+2\sqrt{-C_2}-C_2\end{array}\right.---\left(56\right)$

Easily demonstrate,proveTherefore inequality group (56) can be reduced to

$3-C_2<C_1<1+2\sqrt{-C_2}-C_2---\left(57\right)$

MakePerseverance is set up, then must have C₂＜ -1；

Can to sum up obtain：

If C₁=2-2C₂And C₂＜ -1, then TSG guidance systems stabilization, and in aircraft close to acceleration instruction during target Do not dissipate；If 3-C₂＜ C₁＜ 2-2C₂And C₂＜ -1, then TSG guidance systems stabilization, and in aircraft close to during target Acceleration instruction converges on 0；

C₁And C₂Span extend to plane domain, and Guidance Law coefficient by the zone of dispersion of traditional trajectory shaping Guidance Law C₁And C₂Span not by aircraft speed change influenceed；

Step 5：Terminal velocity control program

The command acceleration a that Section 2 terminal velocity control program described in step 2 is produced_speed, this acceleration effect is in locality In horizontal plane, and perpendicular to speed, it is shown below

$a_{speed}=\mathrm{sgn}({\mathrm{\&psi;}}_0-{\mathrm{\&psi;}}_{LOS0})k_{Vf}{g}_0{\left(\frac{R}{R_0}\right)}^2\left[\begin{array}{c}-sin\left(\mathrm{\&psi;}\right)\\ \cos \left(\mathrm{\&psi;}\right)\\ 0\end{array}\right]---\left(58\right)$

Wherein, k_VfIt is positive parameter undetermined, g₀It is the acceleration of gravity of sea level altitude, R₀Start time aircraft with target it Between remaining flying distance, ψ₀It is the course angle of initial time, ψ_LOS0It is the sight line azimuth of initial time；Sgn (x) is symbol Function, it is as follows

$\mathrm{sgn}\left(x\right)=\left\{\begin{array}{cc}1& x\&GreaterEqual;0\\ -1& x<0\end{array}\right.---\left(59\right)$

By adjusting k_VfValue, aircraft crossrange maneuvering amplitude can be controlled, so that speed when adjusting aircraft hit； a_speedWith a_TSGThe latter end of aircraft flight is acted on simultaneously；a_speedSize gradually weaken close to target with aircraft, The a in aircraft soon hit can so be avoided_speedTo a_TSGInterference；

Aircraft predicted the speed V of hit by Ballistic Simulation of Underwater before the terminal guidance stage is entered_f(K), then press Secant methods correct k_VfValue, finally obtain k_VfIdeal value, it is as follows

$k_{Vf}(K+1)=k_{Vf}\left(K\right)-\frac{\left(V_f\right(K)-V_{fcmd})\left(k_{Vf}\right(K)-k_{Vf}(K-1\left)\right)}{\left(V_f\right(K)-V_f(K-1\left)\right)}---\left(60\right)$

Wherein, V_{f cmd}It is desired terminal velocity, K is iterations.