CN106295218B

CN106295218B - A kind of quick determining energetic optimum intercepts the numerical optimization of predicted set-forward position

Info

Publication number: CN106295218B
Application number: CN201610694774.8A
Authority: CN
Inventors: 周浩; 杜文豪; 陈万春; 杨良
Original assignee: Beihang University
Current assignee: Beihang University
Priority date: 2016-08-19
Filing date: 2016-08-19
Publication date: 2019-04-19
Anticipated expiration: 2036-08-19
Also published as: CN106295218A

Abstract

The invention discloses the numerical optimizations that a kind of quick determining energetic optimum intercepts predicted set-forward position, the specific steps are as follows: (1), gives the exploration value of independent variable X, and set coefficient matrix u and v；(2), partial derivative matrix FF ' (X) is calculated by formula (22) and formula (24)；(3), new independent variable value X is calculated by formula (25)_k+1；(4), judge whether to restrain；If Δ XX=XX_k+1‑XX_kLess than the limits of error, calculating terminates；Otherwise step (2) are gone to, step (2)-(4) are repeated.Based on Gauss problem of the present invention by Kepler orbit theory and based on P iterative method, target is had that the interception of orbit maneuver becomes the amendment of rail minimum speed, is converted into optimum programming problem, and obtain and Δ v by optimum programming theory_minCorresponding KKT condition, while being rapid solving KKT condition, Newton iteration is further devised, the gradient information of parsing is obtained.Relative to general searching algorithm, computational accuracy of the present invention is high, and speed is fast, and insensitive to initial value.

Description

A kind of quick determining energetic optimum intercepts the numerical optimization of predicted set-forward position

Technical field

The present invention relates to spacecrafts to intercept intersection control field, and in particular to a kind of quick determining energetic optimum interception prediction The numerical optimization of the point of impact.

Background technique

In general Orbit Interception, a basic problem for needing to solve is in the given time, how from sky Between in a little reach in addition, the two-point boundary value problem in this kind of astrodynamics is referred to as Lambert problem.Mathematician Gauss gives the classical of Lambert problem and solves, as shown in Figure 1, Gauss problem is defined as: given initial position vector r₁, terminal location vector r₂, and from r₁To r₂Flight time t_FAnd the direction of motion, seek initial and terminal velocity vector ν₁With ν₂。

Under the precondition of Kepler track movement, position vector r₁、r₂And velocity vector ν₁、ν₂It is coplanar, so having

Wherein

For intercept problems, it is generally known that initial point needs great speed, and aircraft can be made along a Kepler Track slides to terminal.By formula (1), initial velocity can be expressed as

So the solution of Gauss problem can simplify as coefficient f, g, andWithSolution.By the limit of transcendental equation System, it is necessary to be solved, the general solution of Gauss problem can be summarized as follows with successive approximation method:

1, the initial value of any one unknown number in P, a or Δ E is given；

2, the value of another two unknown number is calculated by formula (2) and (4)；

3, by formula (3) solve time t, and with given flight time t_FCompare, checking computation result；

If 4, calculated result t and given time t_FIt is inconsistent, iteration variable is corrected, and repeat step 2 and step 3.

When changing tracks for target, the numerical optimization problem of energetic optimum predicted set-forward position is determined, It is exactly to intercept to become rail minimum speed amendment problem, as shown in Fig. 2, t₁At the time of indicating kinetic interception weapon E powered phase Burnout, M indicates that Maneuvering Missile, T indicate the position to strike target on ground.P₁It indicates if Maneuvering Missile is in t₁Moment without become rail, Then interceptor meeting and bullet bump against in the point, i.e. P₁Indicate former predicted set-forward position.If bullet is in t₁Moment is by certain anti-plan of dashing forward Slightly carry out orbit maneuver, it is assumed that the track orb1 that target M is indicated from solid line changes the track orb2 indicated to chain-dotted line, hits mesh Mark is also modified therewith to T ', then current time t₁Blocker E can not be with present speed v₀, blocked on new target track orb2 Cut target M.Assuming that new point of impact position is P₂, then the initial velocity of aircraft must be modified to v₁, v₁Solution can be passed through Gauss problem obtains.

Interception different from the past intersects problem, and when target is there are when orbit maneuver, blocker selects different predictions to order Midpoint P₂Different speed will be brought to correct Δ v=v₁-v₀.For the limited blocker of fuel, speed correction amount Δ v=| Δ v | the smallest maneuver strategy is only optimal, predicted set-forward position P at this time₂With flight time t_FIt also is optimal respectively Intercepting position and interception time.

For the numerical optimization problem of above-mentioned determining energetic optimum predicted set-forward position, traditional method for solving is clock synchronization Between t be iterated search.

P iterative method assumes that the exploration value of a P, and thus value calculates other two unknown number a and Δ E, then solves T out, and it and given flight time t_FCompare, examines exploration value whether suitable with this.P iterative method is it can be concluded that t becomes with P Change the analytical expression of the slope of curve, therefore the exploration value of available Newton iteration method amendment P, improves iteration speed.

Define three constants

By Kepler orbit theory, a can be expressed as

Further, by formula (2), (3), (4), time t can be found out

In addition, t is with the analytical expression of P change curve slope

So can be expressed as using the step of P solution by iterative method Gauss problem

1, by formula (7), computational constant k, l and m；

2, the exploration value of independent variable P is given；

3, by formula (8), the value of a is calculated；

4, by formula (2), (3), (4), calculate f, g andValue；

5, by formula (9), calculate time t, and with given flight time t_FCompare；If equal, terminate；Not phase Deng going to step 6；

6, in conjunction with formula (10), and using the new exploration value of Newton iteration method calculating independent variable P, and 3 are gone to step.

The present invention to be inquired into the numerical optimization problem for surely measuring the optimum prediction point of impact really, and motion model can The earth's core is represented to be reduced to O in Fig. 3；r₁It is the initial position vector of interceptor, r₂It is the point of impingement of interceptor and target, that is, blocks Cut the terminal location vector of bullet；v₀It is the velocity vector before interceptor is motor-driven, is represented by dashed line；v₁It is by solving Gauss problem Obtained initial velocity, it is indicated by the solid line；r₂With v₁Subscript represent the different point of impingement and corresponding interceptor at the beginning of Beginning speed.

Under conditions of known target orbit information, right ascension of ascending node Ω, orbit inclination angle i, the perigee width of target track Angle ω, semi-major axis a_MAnd eccentric ratio e_MAnd logical time of perigee passage t_pMThe form for orbital tracking can be arranged

[Ω,i,a_M,e_M,ω,t_pM]

Terminal location vector r₂By given flight time t_FLimitation.Hereinafter, flying under conditions of not giving rise to misunderstanding by given The row time referred to as intercepts time t.Wherein, the eccentric anomaly E of target movement_MWith true anomaly f_MWith the iteration letter for intercepting time t Number is

Wherein, the orbital period of target is T_M。

So terminal location vector r₂It can be expressed as

So far, it is already possible to which solving target by conventional method, there are the interception of orbit maneuver change rail minimum speed amendments to ask Topic.Conventional method needs to design double iteration: internal layer iteration can obtain a certain prediction by P solution by iterative method Gauss problem Point of impact P₂Corresponding speed correction amount Δ v；External iteration is then by calculating different P₂Corresponding Δ v, so that it is determined that minimum Speed correction amount Δ v_min.Calculation flow chart is as shown in figure 4, wherein external iteration accelerates to restrain using split-half method.

Although this kind of searching algorithm thinking is simple, computational accuracy is also able to satisfy engineering reality, takes a long time, even if making With split-half method acceleration search, also it is unfavorable for line solver.

Summary of the invention

In order to solve the above technical problem, the present invention provides the numbers that a kind of quick determining energetic optimum intercepts predicted set-forward position It is worth optimization method, it is theoretical based on optimum programming, mathematical modeling is carried out by correcting problem to interception change rail minimum speed, is passed through The a large amount of derivation of equation obtains Δ v_minCorresponding KKT condition, and by design Newton iterative, obtain the ladder of parsing Information is spent, is realized to Δ v_minRapid solving, while also solve Optimal Intercept position and intercept the time.This method calculates Precision is high, and speed is fast, and insensitive to initial value.

The following technical solution is employed by the present invention:

A kind of quick determining energetic optimum intercepts the numerical optimization of predicted set-forward position, the specific steps are as follows:

(1), the exploration value of independent variable X is given, and sets coefficient matrix u and v；

(2), partial derivative matrix FF ' (X) is calculated by formula (22) and formula (24)；

(3), new independent variable value X is calculated by formula (25)_k+1；

(4), judge whether to restrain；If Δ XX=XX_k+1-XX_kLess than the limits of error, calculating terminates；Otherwise step is gone to (2), repeat step (2)-(4)；

Wherein, in step (1), independent variable X=[the t P f_M]^T, introduce two groups of equality constraints

Keep the speed correction amount of interceptor minimum, is described as mathematic(al) representation, then objective function is

J (X)=(v_1x-v_0x)²+(v_1y-v_0y)²+(v_1z-v_0z)²；

So far, target is had that the interception of orbit maneuver becomes the amendment of rail minimum speed, is converted into one and contains equation The optimum programming problem of constraint, the expression formula of the optimum programming problem arrange as follows:

X=[t P f_M]^T

J (X)=(v_1x-v_0x)²+(v_1y-v_0y)²+(v_1z-v_0z)²

If the augmented objective function with Lagrange multiplier is

L=J+ λ₁f₁+λ₂f₂ (14)

By KKT condition, can obtain

By formula (15) and (16), solve

Formula (18) are substituted into formula (17), obtain equation

So KKT condition is eventually converted into three equality constraints

If F (X)=[f₁ f₂ f₃]^T (21)

Then the solution of equation group (21) is exactly to intercept the optimal solution for becoming rail minimum speed amendment problem；

In step (2), the formula (22) is the Jacobi matrix of equation group F (X) specifically:

Wherein, because of f₂In not aobvious contain P, f₃In it is not aobvious contain t, soThe then iterative formula of independent variable X It can be written as: X_K+1=X_K-[F′(X_K)]^-1·F(X_K) (23)

To improve the stability that numerical value calculates, formula (22) can be converted into formula (24),

The formula (24) are as follows:

Wherein, diagonal matrix u and v is coefficient matrix described in step (1)；

The formula (25) are as follows:

Method as described above, it is preferable that in the formula (22), every analytical expression is as follows:

Because

And, it is clear that

So need to only obtain described J, f₁、f₂To P, f_MSingle order local derviation, the pure local derviation of second order and second order mixing local derviation, It can be obtained the analytical expression of formula (22)；

By vector project in earth right angle coordinate system, and arrange

Then

Firstly, providing J to P, f_MSingle order local derviation, the pure local derviation of second order and second order mixing local derviation,

If

r_1h2=r_1xh_2x+r_1yh_2y+r_1zh_2z；

Then

Because

So if

Then

Because

So

Because

So if

Then

Secondly, providing f₁To P, f_MSingle order local derviation, the pure local derviation of second order and second order mixing local derviation；Because

χ₀=(2m-l²)P²+2klP-k²；

Also, it sets

So

If

χ₅₁=χ₆₁·2P²+χ₆₂·2·P-2·r₂·χ₆·χ₆₃；

χ₆=r₁-r_1h2；

Then

Because

So

Because

So

Finally, providing f₂To f_MSingle order local derviation, second order local derviation；

Because

So

So far, by formula (30) to formula (41), and formula (26) and formula (27) are combined, obtains the parsing of formula (22) Expression formula, algorithm completely execute.

Based on Gauss problem of the present invention by Kepler orbit theory and based on P iterative method, by target, there are tracks Motor-driven interception becomes rail minimum speed and corrects problem, is converted into optimum programming problem, and obtained by optimum programming theory and Δv_minCorresponding KKT condition, while being rapid solving KKT condition, Newton iteration is further devised, parsing is obtained Gradient information.Relative to general searching algorithm, this method computational accuracy is high, and speed is fast, and insensitive to initial value.

By comparative test, using a kind of acceleration search algorithm based on split-half method, although should the acceleration based on split-half method Searching algorithm is largely increased relative to general traversal formula searching algorithm, computational efficiency, but is still of the invention time-consuming 30 times.The present invention is also compared with the MATLAB majorized function fmincon carried, is as a result also indicated that and is either calculated essence Degree or calculating speed, the algorithm of this paper are all more advantageous.

Detailed description of the invention

Fig. 1 is Gauss problem schematic diagram.

Fig. 2 is that target has orbit maneuver, the belligerent schematic diagram of the blocker based on minimum speed correction amount.

Fig. 3 is that simplified target has that the interception of orbit maneuver becomes rail minimum speed amendment schematic diagram.

Fig. 4 is to solve flow chart using the tradition of split-half method.

Fig. 5 is the calculation flow chart of method in the present invention.

Specific embodiment

It should be noted that in the absence of conflict, the features in the embodiments and the embodiments of the present application can phase Mutually combination.It turns next to attached drawing and the present invention will be described in detail in conjunction with the embodiments.

For solving target, there are the numerical optimization sides that the quick determining energetic optimum of orbit maneuver intercepts predicted set-forward position Method, that is, intercept and become rail minimum speed amendment problem, an optimal interception time t is certainly existed, so that target and interception Bullet arrives at Optimal Intercept position r after time t₂, and make interceptor in initial position r₁The speed knots modification at place is most It is small.So unique independent variable that time t is minimum speed amendment problem is substantially intercepted, but the derivation of equation is excessively complicated.? Formula (9) and formula (11), variable P therein and variable f_M, it is to solve for Gauss problem respectively and solves terminal location vector r₂ Important iteration variable in the process.

The independent variable of minimum speed amendment problem is extended into X=[t P f_M]^T；

In combination with formula (9) and formula (11), two groups of equality constraints are introduced

J (X)=(v_1x-v_0x)²+(v_1y-v_0y)²+(v_1z-v_0z)²

So far, target is had that the interception of orbit maneuver becomes the amendment of rail minimum speed, is converted into one and contains equation The optimum programming problem of constraint, the expression formula of the optimum programming problem arrange as follows

X=[t P f_M]^T

J (X)=(v_1x-v_0x)²+(v_1y-v_0y)²+(v_1z-v_0z)²

If the augmented objective function with Lagrange multiplier is

L=J+ λ₁f₁+λ₂f₂ (14)

By KKT condition, can obtain

By formula (15) and (16), can solve

Formula (18) are substituted into formula (17), available equation

So KKT condition is eventually converted into three equality constraints

If F (X)=[f₁ f₂ f₃]^T (21)

Then the solution of equation group (21) is exactly to intercept the optimal solution for becoming rail minimum speed amendment problem.

It is introduced below how to design Newton iteration solve system of equation (21).

The Jacobi matrix of equation group F (X) is

Wherein, because of f₂In not aobvious contain P, f₃In it is not aobvious contain t, so

Then the iterative formula of independent variable X can be written as

X_K+1=X_K-[F′(X_K)]^-1·F(X_K) (23)

To improve the stability that numerical value calculates, it can use coefficient matrix u and v, convert formula (24) for formula (22)

Wherein, diagonal matrix u and v is referred to as the unified normalized coefficient of coefficient matrix and variable X of F (X) magnitude Matrix.

Further, iterative formula (23) can be converted into

So far, the numerical optimization that quick determining energetic optimum of the invention intercepts predicted set-forward position can be illustrated Are as follows:

1, the exploration value of independent variable X is given, and sets coefficient matrix u and v；

2, partial derivative matrix FF ' (X) is calculated by formula (22) and formula (24)；

3, new independent variable value X is calculated by formula (25)_k+1。

4, judge whether to restrain.If Δ XX=XX_k+1-XX_kLess than the limits of error, calculating terminates；Otherwise step 2 is gone to.This The method of invention is realized to Δ v_minRapid solving, while also solve Optimal Intercept position and intercept the time.The algorithm meter It is high to calculate precision, speed is fast, and insensitive to initial value.Calculation flow chart is as shown in Figure 5.

It arranges in formula (22) below, every analytical expression.

Because

And, it is clear that

So need to only obtain J, f₁、f₂To P, f_MSingle order local derviation, the pure local derviation of second order and second order mixing local derviation Obtain the analytical expression of formula (22).

By vector project in earth right angle coordinate system, and arrange

Then

Firstly, providing J to P, f_MSingle order local derviation, the pure local derviation of second order and second order mixing local derviation.

If

r_1h2=r_1xh_2x+r_1yh_2y+r_1zh_2z；

Then

Because

So if

Then

Because

So

Because

So if

Then

Secondly, providing f₁To P, f_MSingle order local derviation, the pure local derviation of second order and second order mixing local derviation.Because

χ₀=(2m-l²)P²+2klP-k²；

Also, it sets

So

Ifχ₅₁=χ₆₁·2P²+χ₆₂·2·P-2·r₂·χ₆·χ₆₃；

χ₆=r₁-r_1h2；

Then,

Because

So

Because

So

Finally, providing f₂To f_MSingle order local derviation, second order local derviation.

Because

So

So far, by formula (30) to formula (41), and formula (26) and formula (27) are combined, formula (22) can be obtained Analytical expression, method can be executed completely.

In conjunction with the annotation in Fig. 2, a kind of possible belligerent situation is described as follows: t₁Moment, the position arrow of blocker E Measuring with velocity vector is respectively

Also, if blocker with target without any motor-driven, the two will be in P₁Point collision.At this point, target M is carried out Orbit maneuver, new orbital tracking are

[Ω,i,a_M,e_M,ω,t_pM]=

[1.25442346395840,1.31152621076928,5319700.61120683,

0.516254059012830,4.56809019516139,1959.60308392449]；

The cycle of operation T of track_M=3861.3607137952.Then blocker is needed speed v₀It corrects to suitable speed v₁, could be in minimum fuel consumption, i.e., under the premise of minimum speed correction amount, the interception target M on new target track.

It is the primary condition of problem above, because the present invention in derivation process, converts band for the engineering problem There is the optimum programming problem of equality constraint, so while being compared with traditional searching method, it also will be soft with MATLAB The included solution function fmincon of part is compared.

Quantity of state X=[t P f_M]^TBound be set to

X_min=[t_min P_min f_{M min}]^T=[100 1e6-π]^T

X_max=[t_max P_max f_{M max}]^T=[800 5e6 π]^T

The selection rule of the present invention and fmincon iterative initial value are as follows: t₀=(t_min+t_max)/2；P₀It continues to use blocker E and becomes rail Preceding semi-focal chord of satellite orbit；f_M0By by t₀It substitutes into formula (11) and solves acquisition.So the initial value of the example is

X₀=[t₀ P₀ f_M0]=[450 2.4473e6-2.8802]

By initial value X₀It substitutes into iterative formula (25), completes to solve.

Whole simulated programs is I5-2410M in CPU, in the personal computer of dominant frequency 2.3GHz and MATLAB2014a It is completed under simulated environment, more efficient translation and compiling environment will improve computational efficiency.

Table 1 is the calculated result comparison of distinct methods.

The computational accuracy of 1 distinct methods of table compares

Although as can be seen that the calculated result of three kinds of methods and the minimum value of objective function be substantially it is consistent, It is that calculated result is substituted into KKT condition (20), it is found that the result of search method cannot meet KKT condition well, so it Computational accuracy is far below method of the invention.Secondly, we compare the computational efficiency of three kinds of algorithms, the acceleration based on split-half method is searched Although the computational efficiency of rope algorithm is much higher than fmincon, calculation method of the invention can be further by time-consuming from 13.47 Millisecond is compressed to 0.490 millisecond, and computational efficiency about improves 27.5 times.

Table 2 is the acceleration search iterative process based on split-half method

The iterative process of 2 search method of table

Number	t(s)	P(m)	f_M(rad)	ΔV(m/s)
					1	450	1.1307e6	-2.8802	1698.3849
2	275	7.8614e6	-2.9912	3698.5361
					3	362.5	2.9945e6	-2.9364	780.7929
4	406.25	1.8565e6	-2.9085	937.5283
					5	384.375	2.3607e6	-2.9225	670.6556
6	373.4375	2.6593e6	-2.9294	665.8208
					7	378.9063	2.5057e6	-2.9260	653.2358
8	376.1719	2.5814e6	-2.9277	655.6299
					9	377.5391	2.5433e6	-2.9268	653.4668
10	378.2227	2.5244e6	-2.9264	653.1118
					11	378.5645	2.5150e6	-2.9262	653.1142
12	378.3936	2.5197e6	-2.9263	653.0981
					13	378.3081	2.5221e6	-2.9263	653.1012
14	378.3508	2.5209e6	-2.9263	653.0987
					15	378.3722	2.5203e6	-2.9263	653.0982
16	378.3829	2.5200e6	-2.9263	653.0981

Table 3 is iterative process of the invention.

The iterative process of the invention of table 3

Number	t(s)	P(m)	f_M(rad)	ΔV(m/s)
					1	450	2.4473e6	-2.8802	333.7781
2	377.3969	2.5651e6	-2.9275	655.3164
					3	378.3854	2.5196e6	-2.9263	653.1484
4	378.3863	2.5199e6	-2.9263	653.0981
					5	378.3866	2.5199e6	-2.9263	653.0981

To sum up, the present invention is based on optimum programming theories, are built by becoming rail minimum speed amendment problem progress mathematics to interception Mould obtains Δ v by a large amount of derivation of equation_minCorresponding KKT condition, and by design Newton iterative, it obtains The gradient information of parsing is realized to Δ v_minRapid solving, while also solve Optimal Intercept position and intercept the time.Through Comparison is crossed, the algorithm computational accuracy is higher, and speed is faster.

Claims

1. the numerical optimization that a kind of quick determining energetic optimum intercepts predicted set-forward position, which is characterized in that specific steps are such as Under:

(3), new independent variable value X is calculated by formula (25)_k+1；

(4), judge whether to restrain；If Δ XX=XX_k+1-XX_kLess than the limits of error, calculating terminates；Otherwise it goes to step (2), weight Step (2)-(4) are carried out again；

J (X)=(v_1x-v_0x)²+(v_1y-v_0y)²+(v_1z-v_0z)²；

So far, target is had that the interception of orbit maneuver becomes the amendment of rail minimum speed, is converted into one and contains equality constraint Optimum programming problem, the expression formula of the optimum programming problem arranges as follows:

If the augmented objective function with Lagrange multiplier is

L=J+ λ₁f₁+λ₂f₂ (14)

By KKT condition, can obtain

By formula (15) and (16), solve

Formula (18) are substituted into formula (17), obtain equation

So KKT condition is eventually converted into three equality constraints

If F (X)=[f₁ f₂ f₃]^T (21)

Wherein, because of f₂In not aobvious contain P, f₃In it is not aobvious contain t, soThen the iterative formula of independent variable X can be with It is written as:

X_K+1=X_K-[F′(X_K)]^-1·F(X_K) (23)

The formula (24) are as follows:

Wherein, diagonal matrix u and v is coefficient matrix described in step (1)；

The formula (25) are as follows:

2. the method as described in claim 1, which is characterized in that in the formula (22), every analytical expression is as follows:

Because

And, it is clear that

So need to only obtain described J, f₁、f₂To P, f_MSingle order local derviation, the pure local derviation of second order and second order mixing local derviation Obtain the analytical expression of formula (22)；

By vector project in earth right angle coordinate system, and arrange

Then

If

r_1h2=r_1xh_2x+r_1yh_2y+r_1zh_2z；

Then

Because

So if

Then

Because

So

Because

So if

Then

χ₀=(2m-l²)P²+2klP-k²；

Also, it sets

So

If

χ₅₁=χ₆₁·2P²+χ₆₂·2·P-2·r₂·χ₆·χ₆₃；

χ₆=r₁-r_1h2；

Then,

Because

So

Because

So

Because

So

So far, by formula (30) to formula (41), and formula (26) and formula (27) are combined, obtains the Analytical Expression of formula (22) Formula, algorithm completely execute.