CN113359819A

CN113359819A - Optimal guidance law with collision angle constraint and acceleration limitation

Info

Publication number: CN113359819A
Application number: CN202110583839.2A
Authority: CN
Inventors: 陈万春; 赵石磊; 杨良
Original assignee: Beihang University
Current assignee: Beihang University
Priority date: 2021-05-27
Filing date: 2021-05-27
Publication date: 2021-09-07
Anticipated expiration: 2041-05-27
Also published as: CN113359819B

Abstract

The invention provides an optimal guidance law with collision angle constraint and acceleration limitation. Firstly, establishing an energy optimal control problem of the residual flight time weight, and converting the energy optimal control problem into a low-order problem through reduced order conversion; then, obtaining an analytical expression of the guidance instruction, wherein the analytical expression can be expressed as a linear combination of zero control miss amount and zero control angle deviation; then, obtaining the maximum acceleration and an analytical expression between the total control energy and the weight coefficient; finally, a method for selecting the optimal weight coefficient is provided, the total control energy is minimized under the condition of meeting the acceleration constraint, the optimal weight coefficient is substituted into the guidance instruction, and the optimal guidance law meeting the impact angle constraint and the acceleration constraint is obtained.

Description

Optimal guidance law with collision angle constraint and acceleration limitation

Technical Field

The invention provides an optimal guidance law with collision angle constraint and acceleration limitation, and belongs to the technical field of aerospace technologies and weapons.

Background

The terminal law has a very important influence on the destruction efficiency of the weapon system. Achieving a proper collision angle is an important goal of a guidance system during an end attack. By adopting the collision angle control, the missile can attack a weak point of a target so as to improve the warhead efficiency to the maximum extent and ensure high killing probability. Furthermore, acceleration constraints are another important constraint for guided missiles because acceleration saturation of the missile near the intercept point can result in a large amount of miss-hits. Therefore, the guidance law considering the impact angle constraint and the acceleration limit has important significance for improving the operational efficiency of the weapon system.

Disclosure of Invention

The invention comprises the following steps:

the method comprises the following steps: and modeling the optimal guidance law of the residual flight time.

In a planar missile-mesh engagement geometry, ignoring gravitational acceleration, the relative motion between the missile and the target can be expressed as

In the formula, r is the distance between the missile and the target, and lambda is the visual line angle of the missile target; v, a and gamma represent velocity, lateral acceleration and heading angle, respectively, and subscripts M and T represent missile and target, respectively;

and

are respectively r, lambda and gamma_MAnd gamma_TDerivative with respect to time. The dynamic behavior of a missile can be generally expressed by an arbitrary order linear equation of state as follows

a_M＝C_Mx_M+D_Mu (3)

Wherein x is_MIs a state variable of the interception bullet control system,

for its derivative with respect to time, u and a_MRespectively, the commanded acceleration and the actual acceleration of the interceptor projectile.

And

is a coefficient matrix related to the state quantity of the interceptor projectile,

and D_ME R is a matrix of coefficients related to the commanded acceleration.

The guidance law derivation herein will be based on a linearized model. In the final guidance stage, the deviation of the missile and the target from the interception triangle is small, and then the nonlinear model can be linearized. In the linear model, the distance xi of the target and the interception bomb in the direction vertical to the reference sight line is taken as a state quantity, and the second derivative of the state quantity to time is taken as

In the formula, k_TIs a component coefficient of the target acceleration in the direction perpendicular to the reference sight line, and has an expression of k_T＝cos(γ_T0+λ₀) Wherein λ is₀Is a reference line of sightAngle, gamma_T0Is a reference target heading angle. The state quantity of the linearization problem is defined as

Wherein,

the first derivative of ξ with respect to time. The kinetic equation in matrix form can be expressed as

Wherein A, B and C are coefficient matrixes, and the expressions are respectively

In the formula, [0] is defined as a zero matrix of an appropriate dimension.

The performance functional of the optimal guidance law is defined as

Wherein, t₀And t_fRespectively an initial time and a terminal time; xi (t)_f)、γ_T(t_f) And gamma_M(t_f) Is xi and gamma respectively_TAnd gamma_MA value at a terminal time;

is a command impingement angle. α and β are weights of miss amount and collision angle deviation; n is a weight coefficient; t is t_go＝t_f-t is the remaining time of flight; and t is the time of flight.

Step two: and solving the optimal guidance law.

Defining a new state quantity Z (t) satisfying

Where τ represents the time variable of the integral, Φ (t)_fT) is the corresponding state transition matrix,

is the terminal constraint value of the state quantity x (t), D is a constant matrix, and the specific value is

Then two components Z of Z (t)₁And Z₂Can be expressed as

Wherein D is_ξAnd D_γRespectively the first and second row elements of the matrix D. The state quantity Z can be obtained from the formula (54)₁And Z₂Respectively is

Then, the performance functional expressed by the equation (8) can be expressed as a new state quantity Z (t)

According to the optimal control theory, the optimal guidance instruction can be obtained as

Wherein,N_ZEMand N_ZEAEIs the guidance coefficient. If the missile has ideal dynamic characteristics, i.e. zero delay and alpha, beta → ∞ time, the coefficient is

And the state quantity Z₁And Z₂Can be expressed as

The analytical expression of the available guidance instruction about the flight time is

Wherein,

is the dimensionless remaining time of flight,

and

is a dimensionless coefficient expressed as

In the formula, κ₀To characterize the dimensionless parameters of the guidance initiation conditions, the expression is k₀＝Z₁₀/(V_MZ₂₀t_f)，Z₁₀And Z₂₀Are each Z₁And Z₂Is started.

Step three: an analytical expression of the maximum value of the acceleration command is derived and a weight coefficient for minimizing the maximum value is obtained.

For ease of analysis, we define the dimensionless guidance instructions as

According to the theorem of extreme value, N is equal to N_extWhen the temperature of the water is higher than the set temperature,

there is one and only one extreme point t₁. According to the initial conditions κ₀Is different from N_extAre respectively as

Wherein N is₁And N₂Two characteristic values of the weight coefficient N are taken, and the value is equal to k₀In connection with, the expression is

When N belongs to N_extThen, can obtain

At t₁Value of (i.e. extremum)

Which is related to the weight coefficient N, expressed as

Wherein,

is also easy to obtain

The value at t ═ 0, i.e. the initial value

Which is also related to the weight coefficient N, expressed as

We define

Has a maximum value of

When N belongs to N_ext，

By

And

is determined by the absolute value of (a), which can be expressed as

For other weight coefficient values, since there is no extreme point, it can be obtained

According to the initial conditions κ₀In the same manner as in the above, the difference,

the analytical expressions of (a) are respectively:

1) if κ₀≤-1，

2) If-1 < kappa₀＜-1/2，

3) If κ₀≥-1/2，

In the formula,

is a function of

The root of (2). According to the expression, can find

Minimum weight factor N_min：

1) If κ₀≤-1，

N_min＝0 (30)

2) If-1 < kappa₀＜-3/4，

3) If-3/4 ≦ κ₀＜-1/2，

4) If κ₀≥-1/2，

In the formula,

is that

The extreme point of (c).

Step four: and solving an analytic expression of the total control energy.

Definition E represents the total control energy of the missile in the whole guidance process, and the expression is

By substituting the formula (19) into the formula (34) and performing analysis and integration, the product can be obtained

Then, the derivative of E to N is

Wherein,

in the formula, the coefficient a₂，a₁And a₀Respectively, are expressions relating to N. By analysis of Θ_aThe solution of 0 can be found if k₀＞κ_max-0.641, E monotonically increasing for any non-negative N; if κ₀＜κ_maxFor N < N_E1Or N > N_E2E monotonically increases, for N_E1＜N＜N_E2E is monotonically decreasing, where N_E1And N_E2Is theta_aTwo non-negative solutions of 0.

Meanwhile, when N is 0, E has a value of

Then E and E can be obtained from formula (35) and formula (38)₀The difference between

Wherein,

coefficient of formula b₂，b₁And b₀Respectively, are expressions relating to N. For any nonnegative N, Θ_bIf > 0 is always true, then E > E can be obtained from the formula (39)₀I.e. E₀Is the minimum value of E.

Step five: and selecting an optimal weight coefficient.

Suppose acceleration is limited to u_limWhen N is equal to 0, the total control energy E is the minimum value, and the guidance coefficient N is the same_ZEMAnd N_ZEAEIs also minimal, which may increase the robustness of the guidance system to errors, so if u is the case_Max(0)≤u_limWe should choose the weighting factor to be N-0.

If u is_Max(N_min)＜u_lim＜u_Max(0) Satisfying the constraint u can be obtained from an analytical expression of the maximum acceleration with respect to N_Max(N)＜u_limThe set of N of (A) is (N)_lim1,N_lim2) Then we should select the value that minimizes E within this set.

From step four, if κ₀≥κ_maxE monotonically increases with increasing N, then we shouldThe selection weight coefficient is N ═ N_lim1. If κ₀＜κ_maxWhen N is less than N_E1Or N > N_E2While, E monotonically increases; when N is present_E1＜N＜N_E2While E monotonically decreases. Thus in the interval (N)_lim1,N_lim2) Above, the minimum value of E may be located at N_lim1，N_lim2Or N_E2Then we should select the one of these three N that minimizes E, i.e. N ═ N_j|min{E(N_j)},j∈{lim1,lim2,E2}。

If the acceleration limit is more severe, so that

We should first choose the weighting factor N-N_minAnd when κ^*When-1/2, the weight coefficient is reset to N-0, which may reduce the maximum value of the guidance command.

The invention has the advantages that:

1. a residual flight time weight optimal guidance law is provided, and a closed-loop solution of the guidance law under any missile dynamics model is deduced.

2. Deducing the maximum acceleration and an analytical expression between the total control energy and the weight coefficient in the performance functional; the influence of the weight coefficients on the guidance instructions can be analyzed based on this expression.

3. A weight coefficient selection method is provided, so that the acceleration saturation can be avoided while the impact angle constraint is met.

Drawings

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

Fig. 2 is a geometric schematic diagram of a plane bullet warfare.

FIG. 3 is θ_aTrace 0.

FIG. 4a shows an acceleration limit of 50m/s²The trajectory of the missile and the target.

FIG. 4b shows the acceleration limit of 60m/s²The trajectory of the missile and the target.

FIG. 4c shows the acceleration limit of 70m/s²The trajectory of the missile and the target.

FIG. 5a shows an acceleration limit of 50m/s²Acceleration curve of the missile.

FIG. 5b shows the acceleration limit of 60m/s²Acceleration curve of the missile.

FIG. 5c shows the acceleration limit of 70m/s²Acceleration curve of the missile.

In the above figures, the symbols and symbols are as follows:

in FIG. 2, X_I-O_I-Y_IIs a cartesian inertial reference frame, the indices M and T representing the missile and target respectively. V, a and gamma represent speed, lateral acceleration and course angle respectively, r is the distance between missile and target, and lambda is the sight angle, namely LOS and X_IThe angle between the axes. X_RThe axis being defined as the LOS along the initial line of sight₀And xi is the relative position of the missile and the target in the vertical X_RThe component of the axis.

Detailed Description

The invention will be further explained in detail with reference to the drawings and the embodiments.

The invention relates to a residual flight time weight optimal guidance law with impact angle constraint and acceleration limitation, which comprises five steps, wherein the specific flow is shown in figure 1, and the six steps are specifically described below.

Considering the plane missile-target engagement geometry as shown in FIG. 2, the missile and target are perpendicular to the initial line-of-sight LOS₀Is defined as a_MNAnd a_TNAnd satisfy

a_MN＝k_Ma_M；k_M＝cos(γ_M0-λ₀) (41)

The relative intercept angle between the missile and the target is defined as gamma_I＝γ_M+γ_TRequiring it to be equal to a given value at the moment of interception

Neglecting gravitational acceleration, missile and eyeThe relative operation model between the targets is

And the kinematic model of the ballistic inclination of the missile and the target can be expressed as

In the terminal guidance stage, the speed changes of the missile and the target are small and can be ignored, once an interception triangle is formed, the speed of the missile and the target are close to a constant value, and the interception moment is a fixed value. In a guidance law implementation, the remaining time of flight can be estimated using the following equation:

the dynamic behavior of a missile can be generally expressed by an arbitrary order linear equation of state as follows

a_M＝C_Mx_M+D_Mu (48)

Wherein x is_MIs the state variable of the interceptor projectile, u and a_MRespectively, the commanded acceleration and the actual acceleration of the interceptor projectile.

And

and D_ME R is a matrix of coefficients related to the commanded acceleration. If the interceptor projectile has ideal dynamic characteristics, then A_M＝B_M＝C _M0 and D_M1 is ═ 1; if the interceptor projectile has first order dynamics, then A_M＝-1/T，B_M＝1/T，C_M1 and D_MWhere T is a time constant.

The guidance law derivation herein will be based on a linearized model. In the final guidance phase, the deviation of the missile and the target from the interception triangle is small, and then the nonlinear models shown in equations (42) and (43) can be linearized. In the linear model, the kinematic differential equation of xi is

Wherein k is_T＝cos(γ_T0+λ₀) The state quantity of the linearization problem is defined as

The kinematic differential equation in the form of a matrix can be expressed as

Wherein

In the formula, [0] is defined as a zero matrix of an appropriate dimension.

The performance functional of the optimal guidance law is defined as

Wherein α and β are weights of miss amount and impingement angle deviation; n is a weight coefficient; t is t_go＝t_f-t is the remaining time of flight.

Step two: and solving the optimal guidance law.

Defining a new state quantity Z (t) satisfying

Wherein, phi (t)_fT) is a state transition matrix corresponding to equation (51),

As is apparent from the formula (54), Z (t) is a two-dimensional state quantity whose first component Z₁Is zero miss control amount (ZEM), the second component Z₂Is the zero angle control deviation (ZEAE), which refers to the minimum distance and the minimum collision angle deviation which can be achieved between the target and the missile under the conditions that the target keeps the current constant value maneuver and the missile does not maneuver respectively.

According to the properties of the state transition matrix

The kinetic equation of the state quantity Z (t) can be expressed as

Two components Z, in particular to Z (t)₁And Z₂The kinetic equation can be expressed as

Then, the performance functional expressed by the equation (53) can be expressed by the new state quantity Z (t)

According to the optimal control theory, the Hamiltonian of the above problem is

Wherein λ is₁And λ₂Is a covariate. By a reduced order transformation, Z₁And Z₂Is state independent, so that the adjoint equation can be simplified to

Then, the covariates can be represented as

Thus, the coupling equation is satisfied

u^*＝arg_u min H (64)

Is optimally controlled to

Wherein λ is₁And λ₂Is represented by the formula (63). Substituting formula (63) for formula (65), and integrating from t to t_fTwo coupled algebraic equations can be obtained as follows:

wherein,

solving for Z₁(t_f) And Z₂(t_f) And substituting the result into formula (65) to obtain the optimal guidance instruction

Wherein N is_ZEMAnd N_ZEAEIs a guidance coefficient expressed as

If the missile has ideal dynamic characteristics, i.e. zero delay, the state quantity Z₁And Z₂Can be simplifiedInto

Will be provided with

And

substituted for formula (67) to give

Then the guidance coefficient shown in equation (69) can be simplified to

Wherein,

and the state quantity Z₁And Z₂Can be expressed as

As can be seen from equation (73), the guidance coefficient is related to α and β. When let α, β → ∞ this means that the tip distance deviation and the angle deviation are forcibly converged to 0, at which time the guidance coefficient becomes a constant value shown in the following formula

By substituting formula (68) for formula (71), the state quantity Z under optimum guidance can be obtained₁And Z₂Is a time differential equation of

Equation (77) is a first order linear time varying ordinary differential equation. In the case of the constant guidance coefficient shown in the formula (76), the system of ordinary differential equations can be obtained as an analytical expression

Wherein,

and Z₁₀And Z₂₀Are each Z₁And Z₂Is started. We introduce a dimensionless parameter

κ＝Z₁/(V_MZ₂t_go) (80)

Using this parameter, equation (78) can be rewritten as

Wherein,

and

is a dimensionless coefficient expressed as

In the formula, κ₀＝Z₁₀/(V_MZ₂₀t_f) Is the initial value of the parameter k. General formula(81) Substituted into formula (68), the analytical expression of the available guidance instruction is

Wherein eta is_N＝(N+1)(N+2)(N+3)，

Is the dimensionless remaining time of flight.

For the purpose of analysis, we define the dimensionless guidance instruction as

According to the extreme theorem, if a function f is continuous over a closed bounded interval, its maximum may exist at an extreme point or boundary point, an extreme point being where its time derivative is zero

For the dimensionless guidance command curve shown in equation (84), 0 and t are easily obtained_fIs its boundary point.

In the interval [0, t_f]Upper continuous microminiature, with derivative of

Then, according to the extreme point condition

The following two extreme points can be obtained:

where the first solution is coincident with the boundary point. Since the dimensionless remaining time of flight is defined within the bounded interval [0,1], the following condition is satisfied:

the second solution can be considered as a boundary point, which corresponds to a time instant of

We define the set of all weight coefficients N satisfying equation (87) as N_ext. Next, we discuss different initial conditions κ₀In case of the set N_extThe value of (c).

If κ₀Formula (87) can be written as-1

0＜N/(N+3)＜1 (89)

Then N is readily obtained_ext(0, infinity). If κ₀Not equal to-1, formula (87) is equivalent to

Or

Wherein,

then we can get

N_ext＝(0,min{N₁,N₂})∪(max{0,N₁,N₂},+∞) (92)

By comparing N₁And N₂The sign and magnitude of (A), while taking into account kappa₀In the case of-1, N is obtained_extIs composed of

When N belongs to N_extSubstitution of the formula (88) into the formula (85) gives

At t₁Has a value of

Wherein,

will be provided with

And

substituted for formula (85), are likewise available

At t-0 and t-t_fHas a value of

Meanwhile, when N → 0,

has a limit value of

Therefore, we can get

Viewed as a

In the "special" case of (1).

We define

Has a maximum value of

Then the guidance instruction u^*The maximum value of (t) can be expressed as

When N belongs to N_ext，

By

And

is determined by the absolute value of (a), which can be expressed as

To obtain

By introducing a function

For arbitrary κ₀And N ∈ N_extIt is obvious that

Substituting the formulae (94) and (96) into the formula (102), xi (N) may be represented as

We analyze the signs and negatives of xi (N) under different initial conditions and obtain them

The analytical expression of (2).

1) If κ₀Less than or equal to-1, obtainable from formula (82)

At the same time, it is clear that 0 < t_go1< 1, then, for

Can obtain the product

Therefore, the temperature of the molten metal is controlled,

has a maximum value of

2) If-1 < kappa₀< -1/2, wherein N is represented by the formula (93)_ext＝(0,max{0,N₁})∪(N₂,+∞)(N₂,+∞)。

When N is within the range of (0, max {0, N)₁}) is easy to obtain

When N is an element of (N)₂And + ∞) can be obtained

At the same time, is easy to obtain

When N is an element of (N)₂,N₃) (109)

And is

When N is an element of (N)₃,+∞) (110)

Wherein,

N₃＝-(3κ₀+1)/(κ₀+1)＞N₂＞0 (111)

thus, it is possible to obtain

In the formula, function f_add(N) and f_sub(N) is defined as

And N is_sub＝(max{0,N₂},N₃]，N_add＝N_ext-N_sub。

For function f_add(N), when N is within the range of (0, max {0, N)₁) } from formula (107)

Then it is obviously available

When N is an element of (N)₂,N₃]For a fixed N, f_subFor parameter k₀Is a derivative of

From N > N₂Can obtain the product

κ₀＞-(N+3/2)/N+3＝κ_min (117)

Then it is determined that,

thus, f_subFor parameter k₀Is satisfied with

At the same time, when κ₀＝κ_min，f_subHas a value of

f_sub|_κmin＝0

Therefore, for

We can get it

f_sub(N)＞f_sub|_κmin＝0 (121)

When N is an element of (N)₃Infinity), we define a function

The derivative to N is

Wherein,

in the formula, N₀＝(2κ₀+1)/(κ₀+1). When k is₀When > -1, N is easily obtained₀< 2, then f (N) the derivative to N satisfies

Meanwhile, when N → ∞ is satisfied, the value of f (N) is

Therefore, we can get f (N) < 0. Then, it is easy to obtain

When N → ∞ is reached,

then g_addHas a value of

If-1 < kappa₀< -1/3 > having N₃> 0, then g_add(N) has a left boundary value of

If κ₀Not less than-1/3, with N₃Less than or equal to 0, then g_add(N) has a left boundary value of

Meanwhile, the formula (108) and the formula (135) can be used to obtain

Then, according to the zero point theorem,

f_add(N) ═ 0 in the interval (max {0, N)₃Infinity) and only one root is present

And are obviously available

Then, the equations (115), (121) and (131) can be used to obtain

At the same time, when N is equal to R⁺-N_ext＝[N₁,N₂]Due to

No extreme point is present, then its maximum value is

In view of the above, it can be seen that,

has a maximum value of

3) If κ₀≥-1/2，N_extIs (0, ∞), so we can get

Furthermore, if-1/2 < kappa₀< -1/3, obviously having N₃> 0, then

When N belongs to (0, N)₃) (136)

And is

When N is an element of (N)₃,+∞) (137)

If κ₀≥-1/3，N₃If < 0, then

When N ∈ (0, ∞) (138)

Then, the formula (103) can be expressed as

Then it can be obtained

Has a maximum value of

In summary, we have derived different initial conditions

According to which the expression can be derived

Minimum weight factor N_min。

1) If κ₀≦ 1, obtainable from formula (106)

The derivative to N is

For the

Apparently satisfy

Thus making

The smallest weight coefficient is

N_min＝0 (142)

For kappa₀In the case of > -1, we first analyze the function

In the interval (max {0, N)₂}, + ∞).

The derivative to N is

We define

If κ₀Not equal to-1, formula (124) can be rewritten as

Wherein N is₀＝(2κ₀+1)/(κ₀+1). Then the first and second derivatives of Γ (N) over N may be expressed as

D is easily obtained from the formula (147)²Γ/dN²Is less than 0. Meanwhile, according to the formula (146), when N → ∞ sets the value d Γ/dN to 0, d Γ/dN > 0, i.e., Γ (N) monotonically increases. Next we analyze the boundary values of Γ (N)The right boundary value of which is

For the left boundary, according to the difference of the left boundary points, they can be respectively expressed as

Therefore, according to the zero-point theorem, Γ (N) is in the interval (max {0, N)₂' + ∞) and only one root, we define it as

Is also readily available when k₀At the time of > -1,

then

It is obvious that

The minimum value point of (c).

2) If-1 < kappa₀< -1/2, adding N₃By substitution of formula (144) with Γ (N) in N₃Has a value of

Wherein x is ═ k₀+1)/(2κ₀) E (0, 1/2). Since Γ (N) monotonically increases with increasing N, then there must be

Meanwhile, it is apparent from the formula (131)

Thus, it is possible to obtain

Combined type (134) is easy to know

When the temperature of the water is higher than the set temperature,

monotonically decreasing;

when the temperature of the water is higher than the set temperature,

monotonically increasing. When in use

The derivative to N is

It is clear that,

satisfy the requirement of

Thus, it is possible to provide

In the interval

Monotonically increasing, then it is clear

Is that

In the interval N e (N)₂And ∞) minimum point.

For N e [0, N₂]In the case of (1), the formula (134) can give

The derivative to N is

It is apparent that equation (154) exists and that there is only one zero point

Then can be obtained as

When the temperature of the water is higher than the set temperature,

monotonically increasing; when in use

When the temperature of the water is higher than the set temperature,

monotonically decreasing. Thus, in the interval [0, N₂]In the interior of said container body,

the minimum value point of (a) may be N-0 or N-N₂。

And

the difference is

If-1 < kappa₀< -3/4, obtainable from formula (156)

Therefore, in the whole non-negative interval,

may be N-0 or

Then it is available

If-3/4 ≦ κ₀< -1/2, obtainable from formula (156)

Then

Therefore, the temperature of the molten metal is controlled,

3) if κ₀Not less than-1/2, when

When the temperature of the water is higher than the set temperature,

the derivative to N is also of the formula (152), obviously

Therefore, if

When in use

When the temperature of the water is higher than the set temperature,

monotonically decreasing;

when the temperature of the water is higher than the set temperature,

monotonically increasing. If it is not

When in use

When the temperature of the water is higher than the set temperature,

monotonically decreasing;

when the temperature of the water is higher than the set temperature,

monotonically increasing. Thus, we can obtain

Step four: and solving an analytic expression of the total control energy.

Substitution of formula (85) into (162) gives

The derivative of which to the weight coefficient N is

Wherein,

and is

As shown in formula (164), the sign of dE/dN is represented by a quadratic polynomial Θ_aAnd (6) determining. As shown in fig. 3, in two dimensions (N, κ)₀) Equation theta_aThe trace of 0 is a half-open curve that will be in space (N, κ)₀) Two regions are divided into dE/dN < 0 and dE/dN > 0. As can be seen from FIG. 3, in equation Θ_aOn the estimated value of 0, k₀There is a maximum value k_max. Calculate Θ _a0 or so pair of two sides kappa₀Is derived from

It is well known that d κ is and only₀when/dN is 0, k₀Has a value of the maximum value k_maxThus, can obtain

At the same time, since κ_maxIs the trajectory Θ_aWhen the value is 0, the result is

Then, by solving the system of two-dimensional quadratic equations shown in equations (168) and (169), κ can be obtained_max-0.6418. Therefore, if κ₀＞κ_maxFor any non-negative N, one can obtain

dE/dN＞0 (170)

If κ₀＜κ_max，Θ_a(N) there are two non-negative solutions N when 0_E1And N_E2We do not assume N_E1＜N_E2Then can obtain

dE/dN > 0 if N < N_E1Or N > N_E2 (171)

And

dE/dN > 0 if N_E1＜N＜N_E2. (172)

When N is 0, E has a value of

Then E and E are obtained from formula (163) and formula (173)₀The difference between

Wherein,

coefficient of formula b₂，b₁，b₀Are respectively as

Second degree polynomial theta_bSatisfy the discriminant of

Δ＝-(N-1)²(4N²+8N+3)＜0 (177)

Thus, for any nonnegative N, Θ_bIf > 0 is always true, then E > E can be obtained from equation (174)₀I.e. E₀Is the minimum value of E.

Step five: and selecting an optimal weight coefficient.

Suppose acceleration is limited to u_lim. And step five, when N is equal to 0, the total control energy E is the minimum value, and the guidance coefficient N is simultaneously used_ZEMAnd N_ZEAEIs also minimal, which may increase the robustness of the guidance system to errors, so if u is the case_Max(0)≤u_limWe should choose the weighting coefficients as

N＝0 (178)

If u is_Max(N_min)＜u_lim＜u_Max(0) The satisfaction of the constraint u can be obtained from an analytical expression of the maximum acceleration with respect to the weight coefficient N_Max(N)＜u_limThe set of N of (A) is (N)_lim1,N_lim2) Then we should choose to select the value that minimizes E within this set. From step five, if κ₀≥κ_maxE increases monotonically with increasing N, then we should choose the weighting factor as

N＝N_lim1 (179)

If κ₀＜κ_maxWhen N is less than N_E1Or N > N_E2While, E monotonically increases; when N is present_E1＜N＜N_E2While E monotonically decreases. Thus in the interval (N)_lim1,N_lim2) Above, the minimum value of E may be located at N_lim1，N_lim2Or N_E2Then we should select the value of these three N that minimizes E, i.e. the value

N＝N_j|min{E(N_j)},j∈{lim1,lim2,E2} (180)

If the acceleration limit is more severe, so that

We should first choose the weighting coefficients as

N＝N_min (181)

In this case, formula (81) is substituted for formula (80), and the value of parameter κ corresponding to the optimal trajectory is

Let kappa^*When available at this time, it is called-1/2

Has a value of

This time can be obtained by substituting formula (183) into formula (81)

Has a value of

If we reselect N to 0 at this time, then the subsequent guidance command is constant as given by equation (83)

We define

To compare

And

by the size relationship of (1), we introduce a variable

The expression is

The derivative to N is

Eta is obtained by formula (187)_u> 0, therefore for the analysis of d η_uSign of dN, we only need to analyze

The derivative to N is

From equation (190), d Λ/dN has only one non-negative zero, and when N is less than this value, d Λ/dN > 0; when N is larger than this value, d Λ/dN is less than 0. Then, we can conclude

Λ(N)≥min{Λ(0),Λ(∞)}＝min{3/2-ln3,0}＝0 (191)

Thus, it is possible to obtain

η_u≥η_u(0)＝1 (192)

Thus, we can be best at κ^*At-1/2, if we reselect the weighting factor to 0, the maximum value of the guidance command will be reduced, which is beneficial to avoid acceleration saturation. Then for

In this case, we should first choose the weighting factor N ═ N_minAnd when κ^*When the value is-1/2, the weight coefficient is reset to N-0.

Example of the implementation

In order to verify the guidance method provided by the invention, numerical simulation is carried out by taking extra-atmospheric interception as a background. The missile had first order dynamics with a time constant of 0.3s, and the simulated initial conditions are shown in table 1. During simulation, the missile-eye relative motion model adopts nonlinear models shown in formulas (42) and (43), and gamma used for realizing guidance law in the flight process_M0，γ_T0And λ₀The updating is continuous, but the weighting factor N in the performance functional is not updated in real time.

TABLE 1 parameter settings in end-guidance simulation

Parameter(s)	Value taking	Parameter(s)	Value taking
				Initial bullet eye relative distance (km)	400	Target initial velocity (m/s)	5000
Initial line of sight angle (deg)	0	Target initial trajectory inclination (deg)	0
				Missile initial velocity (m/s)	5000	Target maneuvering acceleration (m/s)²)	10
Initial trajectory inclination of missile (deg)	0	End impingement Angle restraint (deg)	0

FIGS. 4 a-4 c show acceleration limits of 50m/s, respectively²，60m/s²，70m/s²The trajectories of the missile and the target, as can be seen in the figure, the missile can directly hit the target even if there is an acceleration limit. 5 a-5 c show the missile acceleration curves for three cases, respectively, and the proposed guidance law ensures that the acceleration in the whole flight is within the limit by selecting a proper weight coefficient N as shown in the figure. Therefore, the guidance law proposed by the patent can improve the effectiveness of direct hits compared to other guidance laws.

Claims

1. An optimal guidance law with collision angle constraints and acceleration limits, comprising the steps of:

the method comprises the following steps: modeling the optimal guidance law of the residual flight time;

in the plane missile-mesh engagement geometry, ignoring gravitational acceleration, the relative motion between the missile and the target is represented as

and

are respectively r, lambda and gamma_MAnd gamma_TA derivative with respect to time; the dynamic characteristics of the missile are expressed by the following linear state equation of any order

a_M＝C_Mx_M+D_Mu (3)

Wherein x is_MIs a state variable of the interception bullet control system,

for its derivative with respect to time, u and a_MRespectively an instruction acceleration and an actual acceleration of the interception bomb;

and

and D_ME R is a coefficient matrix related to the command acceleration;

guidance law derivation will be based on a linearized model; in the final guidance stage, if the deviation of the missile and the target from the interception triangle is very small, the nonlinear model is linearized; in the linear model, the distance xi of the target and the interception bomb in the direction vertical to the reference sight line is taken as a state quantity, and the second derivative of the state quantity to time is taken as

In the formula, k_TIs a component coefficient of the target acceleration in the direction perpendicular to the reference sight line, and has an expression of k_T＝cos(γ_T0+λ₀) Wherein λ is₀To reference the angle of sight, gamma_T0Is a reference target course angle; the state quantity of the linearization problem is defined as

Wherein,

is the first derivative of ξ versus time; the kinetic equation in matrix form is expressed as

Wherein [0] is defined as a zero matrix of suitable dimensions;

the performance functional of the optimal guidance law is defined as

is a command impingement angle; α and β are weights of miss amount and collision angle deviation; n is a weight coefficient; t is t_go＝t_f-t is the remaining time of flight; t is the time of flight;

step two: solving the optimal guidance law;

defining a new state quantity Z (t) satisfying

Then two components Z of Z (t)₁And Z₂Is expressed as

Wherein D is_ξAnd D_γFirst and second row elements of matrix D, respectively; quantity of state Z₁And Z₂Respectively is

Then, the performance functional expressed by the equation (8) is expressed by the new state quantity Z (t)

According to the optimal control theory, the optimal guidance instruction is obtained as

Wherein N is_ZEMAnd N_ZEAEIs the guidance coefficient; if the missile has ideal dynamic characteristics, i.e. zero delay and alpha, beta → ∞ time, the coefficient is

And the state quantity Z₁And Z₂Is shown as

The analytic expression of the obtained guidance instruction about the flight time is as follows

Wherein,

is the dimensionless remaining time of flight,

and

is a dimensionless coefficient expressed as

In the formula, κ₀To characterize the dimensionless parameters of the guidance initiation conditions, the expression is k₀＝Z₁₀/(V_MZ₂₀t_f)，Z₁₀And Z₂₀Are each Z₁And Z₂An initial value of (1);

step three: deducing an analytical expression of the maximum value of the acceleration command and solving a weight coefficient which enables the maximum value to be the minimum;

for the purpose of analysis, dimensionless guidance instructions are defined as

Obtaining the current N belongs to N according to the extreme value theorem_extWhen the temperature of the water is higher than the set temperature,

there is one and only one extreme point t₁(ii) a According to the initial conditions κ₀Is different from N_extAre respectively as

Wherein N is₁And N₂Two features for the weight coefficient N are takenValue, its value and kappa₀In connection with, the expression is

When N belongs to N_extWhen it is obtained

At t₁Value of (i.e. extremum)

Which is related to the weight coefficient N, expressed as

Wherein,

is also easy to obtain

The value at t ═ 0, i.e. the initial value

Which is also related to the weight coefficient N, expressed as

Definition of

Has a maximum value of

When N belongs to N_ext，

By

And

is expressed as

For other weight coefficient values, the extreme point does not exist, and then the value is obtained

the analytical expressions of (a) are respectively:

if κ₀≤-1，

If-1 < kappa₀＜-1/2，

If κ₀≥-1/2，

In the formula,

is a function of

The root of (2); according to the expression, find

Minimum weight factor N_min：

If κ₀≤-1，

N_min＝0 (30)

If-1 < kappa₀＜-3/4，

If-3/4 ≦ κ₀＜-1/2，

If κ₀≥-1/2，

In the formula,

is that

The extreme point of (a);

step four: solving an analytical expression of the total control energy;

Substituting the formula (19) into the formula (34), and analyzing and integrating the result

Then, the derivative of E to N is

Wherein,

in the formula, the coefficient a₂，a₁And a₀Respectively, expressions related to N; by analysis of Θ_aA solution of 0; if κ₀＞κ_max-0.6418, E monotonically increasing for any non-negative N; if κ₀＜κ_maxFor N < N_E1Or N > N_E2E monotonically increases, for N_E1＜N＜N_E2E is monotonically decreasing, where N_E1And N_E2Is theta_aTwo non-negative solutions of 0;

meanwhile, when N is 0, E has a value of

Then E and E are obtained from the formulae (35) and (38)₀Difference of difference

Wherein,

coefficient of formula b₂，b₁And b₀Respectively, expressions related to N; for any nonnegative N, Θ_bIf > 0 is always true, then E > E is obtained from formula (39)₀I.e. E₀Is the minimum value of E;

step five: selecting an optimal weight coefficient;

let the acceleration limit be u_limWhen N is equal to 0, the total control energy E is the minimum value, and the guidance coefficient N is the same_ZEMAnd N_ZEAEIs also minimal, increasing the robustness of the guidance system to errors, so if u is_Max(0)≤u_limThe weight coefficient should be chosen to be N ═ 0;

if u is_Max(N_min)＜u_lim＜u_Max(0) The satisfaction of the constraint u is found from an analytical expression of the maximum acceleration with respect to N_Max(N)＜u_limThe set of N of (A) is (N)_lim1,N_lim2) Then the value that minimizes E within this set should be selected;

from step four, if κ₀≥κ_maxE monotonically increases with increasing N, the weighting factor should be chosen to be N ═ N_lim1(ii) a If κ₀＜κ_maxWhen N is less than N_E1Or N > N_E2While, E monotonically increases; when N is present_E1＜N＜N_E2While, E is monotonically decreasing; thus in the interval (N)_lim1,N_lim2) Above, the minimum value of E may be located at N_lim1，N_lim2Or N_E2Then the value of these three N that minimizes E should be selected, i.e., N-N_j|min{E(N_j)},j∈{lim1,lim2,E2}；

If the acceleration limit is more severe, so that

The weighting factor N-N should first be chosen_minAnd when κ^*When the value is-1/2, the weight coefficient is reset to N-0, which can reduce the maximum value of the guidance command.