CN114935277A - Online planning method for ideal trajectory of gliding extended-range guided projectile - Google Patents

Online planning method for ideal trajectory of gliding extended-range guided projectile Download PDF

Info

Publication number
CN114935277A
CN114935277A CN202210212623.XA CN202210212623A CN114935277A CN 114935277 A CN114935277 A CN 114935277A CN 202210212623 A CN202210212623 A CN 202210212623A CN 114935277 A CN114935277 A CN 114935277A
Authority
CN
China
Prior art keywords
axis
coordinate system
guided projectile
point
projectile
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202210212623.XA
Other languages
Chinese (zh)
Other versions
CN114935277B (en
Inventor
易文俊
孙瑞雪
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing University of Science and Technology
Original Assignee
Nanjing University of Science and Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nanjing University of Science and Technology filed Critical Nanjing University of Science and Technology
Priority to CN202210212623.XA priority Critical patent/CN114935277B/en
Publication of CN114935277A publication Critical patent/CN114935277A/en
Application granted granted Critical
Publication of CN114935277B publication Critical patent/CN114935277B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F41—WEAPONS
    • F41G—WEAPON SIGHTS; AIMING
    • F41G3/00—Aiming or laying means
    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F42—AMMUNITION; BLASTING
    • F42B—EXPLOSIVE CHARGES, e.g. FOR BLASTING, FIREWORKS, AMMUNITION
    • F42B35/00—Testing or checking of ammunition
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Abstract

The invention relates to an online planning method for ideal trajectory of gliding extended range guided projectile. The method comprises the following steps: establishing a ground, translation, projectile body, trajectory and speed coordinate system; establishing a 3dof guided projectile motion equation; carrying out a wind tunnel test to obtain pneumatic parameters; establishing a performance index of a trajectory planning problem according to a hitting requirement, and establishing a constraint condition; planning an ideal trajectory at the initial ascending section or the boosting section after launching, acquiring and processing satellite data at the moment, and substituting the data as a starting point into a guided projectile motion equation; and reconstructing the grids by using a self-adaptive hp method, optimizing the number of the grids and the polynomial order of interpolation, converting the problem of the pseudo-spectral method discrete optimal control into a nonlinear programming problem, and solving to obtain an ideal trajectory curve, a speed curve, a trajectory dip angle curve, an overload curve and a rudder deflection angle curve of the guided projectile centroid motion. The invention is insensitive to the interference of the guided projectile before the online planning of an ideal trajectory, and has high striking precision.

Description

Online planning method for ideal trajectory of gliding extended-range guided projectile
Technical Field
The invention belongs to the field of guided trajectories, and particularly relates to an online planning method for an ideal trajectory of a gliding extended-range guided projectile.
Background
The problem of ballistic planning of guided projectiles is widely applied in the technical field of engineering, and an ideal ballistic refers to a trajectory of the guided projectile from a motion starting point to a motion finishing point, which meets the requirements of striking precision and force. Online planning of ideal trajectories is key to reducing the burden on fire control systems and the difficulty of trajectory tracking. During launching, the deviation of parameters such as firing angle, deflection angle, initial speed and the like of the guided projectile and the interference of a low altitude ascending section can lead the relevant motion parameters of the guided projectile to be different at the moment of starting to plan the trajectory on line. The problem of online planning of ideal trajectory of gliding extended-range guided projectile is a hot spot of research in recent years.
The ballistic programming problem is essentially an optimal control problem with state constraints and control constraints, which satisfies a certain index, and as the system and constraint conditions become more and more complex, the problem is usually solved by a numerical method, which includes an indirect method and a direct method.
In the document, "hypersonic jump-glide trajectory scheme design and optimization" (Zhao Ji Song, and the like, solid rocket technology, vol.32, No. 2, page 123-126, 2009), an indirect method is adopted to convert the trajectory optimization problem into a two-point boundary value problem, and a conjugate gradient method is used for resolving. In the literature, "study on the ascending section track optimization method based on indirect method" (wujia beam, navigation positioning and time service, 3 rd, 2 nd, pages 14-19, 2016), for the problem of optimizing the solid rocket ascending section track, the indirect method is adopted to obtain the optimal track in the real atmospheric environment under the condition of meeting the attack angle constraint. The indirect method mainly comprises the steps of solving a first-order necessary condition of an optimal control problem by applying a Pontryagin extreme value principle to obtain a Hamilton multi-point edge value problem, and enabling a system state, control input and parameters to meet the two-point edge value problem through numerical calculation so as to obtain an optimal solution. However, when the Pontryagin extreme value principle deduces the optimal first-order necessary condition, a covariate is introduced, the method is sensitive to an initial value and has large limitation, and two-point edge values are difficult to solve when a complex optimal control problem is processed, so that the method is limited in practical engineering application.
In the document, "track fast generation under multi-constraint conditions of hypersonic unpowered long-distance gliding aircraft" (Yongmi, etc., astronavigation report, No. 1, page 46-52, 2008), a p method is used for updating a grid, the dimension of an interpolation polynomial is changed to realize convergence, and a serial solving strategy capable of generating an initial value is provided based on a Gaussian pseudo-spectrum method to deal with complex multi-constraint conditions. The direct method comprises a target shooting method and a point matching method, a parameterization method is used for converting a wireless dimensional optimization problem into a finite dimensional nonlinear programming problem, and a nonlinear programming problem solver is used for directly obtaining control variables required by the original optimal control problem, so that the optimal ideal ballistic trajectory is solved and calculated. The target shooting method only discretely controls variables and cannot solve a complex optimal control problem, so that the target shooting method is rarely used; the method is characterized in that a control variable and a state variable are dispersed by a point matching method, wherein a pseudo-spectrum method is widely applied due to good convergence, and a pseudo-spectrum method combined with a p method is difficult to solve discontinuous or non-smooth problems due to node fixation, and dimension disasters can be caused when the dimension is excessively increased.
The existing gliding range-extending guided projectile completes accurate shooting on a target through the control of a fire control system, particularly a carrier-based gliding range-extending guided projectile, completes planning of an ideal trajectory through a carrier-based projectile fire control system before launching, then conducts shooting, achieves a combat mission of remote and accurate shooting on the bank, and generally has accurate requirements on a Circular Error (CEP) of a landing point from the target, such as the requirement of an embodiment that a geostationary index CEP is less than 20m, and the shooting speed is not less than 5 rounds per minute. The number of reserve projectiles on a ship usually reaches the number level of hundreds to thousands according to the battle requirement, when a battle mission is urgent, all ships and cannons on the ship are required to shoot at the allowed maximum shooting speed and finish accurate striking, and under the condition, the planning of an ideal trajectory is carried out on each launched guided projectile, so that the burden is brought to a fire control system.
Disclosure of Invention
The invention aims to provide an online planning method for ideal trajectory of gliding extended-range guided projectile.
The technical solution for realizing the purpose of the invention is as follows: an on-line planning method for ideal trajectory of gliding extended range guided projectile comprises the following steps:
step (1): establishing ground Axyz, translation Oxyz and projectile body Ox 1 y 1 z 1 Trajectory Ox 2 y 2 z 2 Velocity Ox 3 y 3 z 3 A coordinate system;
step (2): establishing a motion equation of the 3dof guided projectile;
and (3): carrying out a wind tunnel test to obtain pneumatic parameters;
and (4): establishing a performance index of a trajectory planning problem according to a hitting requirement, and establishing a constraint condition according to the launching facility capacity of a gun platform, the strength limit of a bullet body and the strength requirement of a hitting target;
and (5): planning an ideal trajectory at the initial ascending section or the boosting section after launching, acquiring and processing satellite data at the moment, and substituting the data as a starting point into a guided projectile motion equation;
and (6): reconstructing grids by using a self-adaptive hp method, optimizing the number of the grids and the polynomial order of interpolation, converting the problem of discrete optimal control by using a pseudo-spectrum method into a nonlinear programming problem, and solving by using a nonlinear programming solver to obtain an ideal ballistic curve, a velocity curve, a ballistic dip angle curve, an overload curve and a rudder deflection angle curve of the guided projectile centroid motion.
Further, the step (1) of establishing the ground Axyz, the translation Oxyz and the projectile body Ox 1 y 1 z 1 Trajectory Ox 2 y 2 z 2 Velocity Ox 3 y 3 z 3 The specific method of the coordinate system is as follows:
the ground coordinate system is fixedly connected with the ground, the origin A is on the center of mass of the projectile body during launching, the Ax axis is a horizontal line connecting the launching device and the target and points to the target, the Ay axis is upward along the vertical line, and the Az axis is vertical to the Ax and Ay;
o is the instantaneous center of mass of the projectile body, and Ox, Oy and Oz are respectively parallel to and in the same direction as Ax, Ay and Az; ox (oxygen) 1 The axis is coincident with the longitudinal axis and points to the warhead as positive, Oy 1 With axes lying in a longitudinal plane, with Ox 1 Vertical, with upward positive, Oz 1 Axis perpendicular to Ox 1 y 1 A plane;
Ox 2 the axis coinciding with the instantaneous velocity V, Oy, of the centre of mass of the guided projectile 2 The axis is located in the vertical plane containing V and the axis x 2 Vertical axis, Oz 2 Axis perpendicular to Ox 2 Axis and Oy 2 A shaft; ox (oxygen) 3 Axis coincident with V, Oy 3 The axis lying in the longitudinal symmetry plane with Ox 3 Vertical axis, Oz 3 Axis perpendicular to Ox 3 y 3 And (4) a plane.
Further, the step (2) of establishing a 3dof guided projectile motion equation comprises the following specific steps:
neglecting the roll of the guided projectile, the kinetic equation set of the center of mass of the guided projectile in the ballistic coordinate system is as follows:
Figure RE-GDA0003758055760000021
the kinematic equation of the center of mass of the guided projectile in the ground coordinate system is as follows:
Figure RE-GDA0003758055760000022
the centroid mass change equation and the governing relation equation are as follows:
Figure RE-GDA0003758055760000031
where m is the mass of the guided projectile,
Figure RE-GDA0003758055760000032
for the rate of change of mass, the boost period is not zero only during operation of the engine, the start and end times of the boost period being dependent on the guided projectile itself, m c The specific value is that the drug loading of the guided cannonball is related to the combustion time of the engine, V is the instantaneous speed of the mass center of the guided cannonball,
Figure RE-GDA0003758055760000033
is the speed change rate, g is the gravity acceleration, P is the average thrust of the engine, and exists only in the boosting period, X is the air resistance borne by the guided projectile, Y is the lift force, Z is the lateral force, and the attack angle alpha is the speed in the projectile coordinate system Ox 1 y 1 z 1 On projection and axis Ox 1 Angle of inclination, alpha c (t) is the value of the attack angle alpha changing with time under the action of the control law, and the sideslip angle beta is the speed and projectile coordinate system Ox 1 y 1 z 1 Angle of (b) of c (t) is the value of the slip angle beta under the action of the control law as a function of time, f α 、f β Based on a relation function of an attack angle, a side slip angle and two rudder deflection angles which are obtained by the instantaneous balance hypothesis expansion and are related to pneumatic parameters, the gravity center of the guided projectile, the pressure center and the pressure center of the control rudder, the guided projectile is expanded to the top point of a trajectory and then the control rudder is expanded for control, and the speed inclination angle gamma is v Is the Oy of the velocity coordinate system 3 The angle between the axis and Ox2y2 of the trajectory coordinate system, the trajectory inclination angle theta is the angle between V and the horizontal plane,
Figure RE-GDA0003758055760000034
is the rate of change of ballistic inclination angle, ballistic declination angle
Figure RE-GDA0003758055760000035
Is the angle between the projection of V on the horizontal plane and the Ox axis,
Figure RE-GDA0003758055760000036
is the rate of change of ballistic inclination, (x, y, z) is the instantaneous position of the center of mass of the guided projectile in the ground coordinate system,
Figure RE-GDA0003758055760000037
the position change rate of the center of mass of the guided projectile on the x axis, the y axis and the z axis of a ground coordinate system is obtained.
Further, the wind tunnel test in the step (3) is carried out to obtain the pneumatic parameters, and the specific method is as follows:
the gliding extended-range guided projectile model is placed in a wind tunnel, a plurality of groups of data of zero lift resistance coefficients, lift coefficient derivatives, pitch rudder and yaw rudder zero lift resistance coefficients, pitch rudder and yaw rudder lift coefficient derivatives, missile wing combination pressure center positions and control rudder pressure center positions corresponding to different Mach numbers are obtained according to the action of gas flow and the guided projectile model, and an induced resistance coefficient and a control rudder induced resistance coefficient are selected.
Further, the step (4) is specifically as follows:
when the range is maximized, selecting the x-axis coordinate of the flight end point of the guided projectile as a final value type performance index, wherein the expression is as follows:
J＝-x(t f )
wherein J represents a performance index, t f Represents the time required to fly off the ideal trajectory, x (t) f ) Representing the coordinates of the ballistic endpoint on the axis of the ground coordinate system Ax;
setting a starting point parameter value of the guided projectile as a starting point constraint, setting the maximum value and the minimum value of a state variable and a control variable in the flight process as a process constraint, and setting the falling angle and the falling speed of the guided projectile when the guided projectile falls to the ground to be not less than a terminal point constraint set by a user.
Further, the step (5) is specifically as follows:
beidou satellite represents the position P point (x) of the guided cannonball based on WGS-84 coordinate system ε ,y ε ,z ε ) Conversion to position data in a geodetic coordinate system
Figure RE-GDA0003758055760000038
Figure RE-GDA0003758055760000041
Wherein a represents the major semi-axis of the WGS-84 ellipsoid and has a value of 6378137.0 m; e is the first eccentricity of WGS-84, squared 0.0066943799; n is half of the curvature of the ellipsoidal unitary fourth of the twelve earthly branches; the origin w of the WGS-84 coordinate system is the earth centroid, the wz axis points to the north pole of an ellipsoid, the wx axis points to the intersection point of the equatorial plane and the Greenwich meridian plane, and the wy axis in the equatorial plane and xwz form a right-hand system coordinate system; the short semi-axis of the ellipsoid of the geodetic coordinate system is coincident with the rotation axis wz of the earth, the longitude lambda is the included angle between the meridian plane of the ellipsoid where the point P is located and the meridian plane of Greenwich mean, and the latitude
Figure RE-GDA0003758055760000042
Is the included angle between the normal line of the ellipsoid passing through the point P and the equatorial plane, and the height h is the distance from the normal line of the ellipsoid to the ellipsoid along the point P;
the coordinates of the guided projectile in the geodetic coordinate system are converted into coordinates (x) in the north-celestial-east coordinate system n ,y n ,z n )：
Figure RE-GDA0003758055760000043
In the formula (x) ε0 ,y ε0 ,z ε0 ) The position of the origin of a north heaven coordinate system in a WGS-84 coordinate system is shown, the origin n of the north heaven coordinate system is a transmitter point, xnz is a ground horizontal plane where the origin is located, an nx axis points to the north of the earth, an nz axis points to the east of the earth, and a ny axis is perpendicular to a xnz plane and points to the air;
the north-east coordinate system to the shooting coordinate system (x, y, z) is:
Figure RE-GDA0003758055760000044
in the formula, α is the shot direction in the north-east-heaven coordinate system, and the shooting coordinate system is the ground coordinate system.
Further, the step (6) is specifically as follows:
obtaining the motion state equation constraint of optimal control according to the dynamics and kinematics equation set and the mass change equation in step (3) and step (2)
Figure RE-GDA0003758055760000045
And (5) obtaining the performance index J and constraint conditions according to the step (4), and finally determining the description model of the optimal control problem as follows:
Figure RE-GDA0003758055760000046
Figure RE-GDA0003758055760000047
Figure RE-GDA0003758055760000048
Figure RE-GDA0003758055760000049
Figure RE-GDA00037580557600000410
wherein J represents a performance index, Φ is a final value type performance index, t f Indicating the time required to fly-off the ideal trajectory,
Figure RE-GDA0003758055760000051
respectively representing vectors formed by all state variables and all control variables, 4, the control variables are derivatives of two rudder deflection angles with respect to time, the initial value of the rudder deflection angle is 0 degree, the numerical value of the rudder deflection angle is obtained through integration,
Figure RE-GDA0003758055760000052
to represent
Figure RE-GDA0003758055760000053
Derivative with respect to time, f is a mapping by which
Figure RE-GDA0003758055760000054
And
Figure RE-GDA0003758055760000055
to represent
Figure RE-GDA0003758055760000056
q is
Figure RE-GDA0003758055760000057
And
Figure RE-GDA0003758055760000058
as a function of (a) or (b),
Figure RE-GDA0003758055760000059
a state vector representing the end point of the trajectory,
Figure RE-GDA00037580557600000510
state vector representing the start of trajectory, phi being an edgeBound event constraint function, including the vector of states x (t) for the starting point 0 ) And an end point state vector x (t) f ) The constraint of (2); c. C 1 And c 2 Are equality and inequality constraints.
Further, the step (6) of reconstructing grids by using a self-adaptive hp method, optimizing the grid number and interpolation polynomial order, and converting the discrete optimal control problem into a nonlinear programming problem by using a pseudo-spectral method specifically comprises the following steps:
step (61): carrying out initial grid division; setting a maximum allowable error epsilon and a maximum allowable curvature r max The initial grid number and the order of the lagrange interpolation polynomial;
will time interval t 0 ,t f ]Divided into K grid subintervals [ t ] k-1 ,t k ]K1, …, K, a new time variable τ e [ -1,1 ] is introduced]Converting the original variable t into a new variable tau through the following formula, namely converting the single-interval nonlinear optimal control problem into a multi-interval nonlinear optimal control problem:
Figure RE-GDA00037580557600000511
step (62): calculating the position of an LGR interpolation point, an LGR weight and a Radau pseudo-spectral differential matrix in each current grid subinterval, approximately expressing an optimal control model by combining a Lagrange interpolation polynomial, and converting the optimal control model into a nonlinear programming problem;
let the kth grid subinterval [ t ] k-1 ,t k ]The state variables and control variables in are denoted X k (τ) and U k (τ)，N k +1 discrete points
Figure RE-GDA00037580557600000512
And end point
Figure RE-GDA00037580557600000513
Form an interval tau epsilon-1, 1]Inner N k +2 Lagrangian interpolation points, wherein
Figure RE-GDA00037580557600000514
Is provided with
Figure RE-GDA00037580557600000515
Is N k A +2 order lagrange interpolation basis function,
Figure RE-GDA00037580557600000516
is N k The +1 order Lagrange interpolation basis function has the expression:
Figure RE-GDA00037580557600000517
Figure RE-GDA00037580557600000518
state variable X (k) In N k +2 discrete points of use
Figure RE-GDA00037580557600000519
The approximation is:
Figure RE-GDA00037580557600000520
in the formula (I), the compound is shown in the specification,
Figure RE-GDA0003758055760000061
is a state variable X k At discrete points
Figure RE-GDA0003758055760000062
Taking the value of (A);
control variable U (k) At N k +1 discrete point application
Figure RE-GDA0003758055760000063
The approximation is:
Figure RE-GDA0003758055760000064
in the formula (I), the compound is shown in the specification,
Figure RE-GDA0003758055760000065
is a state variable U k At discrete points
Figure RE-GDA0003758055760000066
Taking the value of (A);
for state variable X (k) The derivative to get the approximate equation of state is:
Figure RE-GDA0003758055760000067
order to
Figure RE-GDA0003758055760000068
The specific expression of a Radau pseudo-spectral differential approximation matrix representing the kth subinterval is as follows:
Figure RE-GDA0003758055760000069
Figure RE-GDA00037580557600000610
in the formula, g (k) (τ) is P n Function of (τ), P n (τ) is a Legendre orthogonal polynomial, and the differential expression is:
Figure RE-GDA00037580557600000611
then N is k The state equation constraints on +1 interpolation points are converted into algebraic equation constraints expressed as:
Figure RE-GDA00037580557600000612
in the formula F s (k) Is a vector of a function of the state equation with respect to the state variables and the control variables;
the LGR integral is used for approximating an integral term in a Boercus-type performance index, and the approximation is as follows:
Figure RE-GDA00037580557600000613
in the formula (I), the compound is shown in the specification,
Figure RE-GDA00037580557600000614
is a function of the state vector and the control vector,
Figure RE-GDA00037580557600000615
as discrete points of the kth sub-interval
Figure RE-GDA00037580557600000616
The corresponding LGR weights are in the specific form:
Figure RE-GDA0003758055760000071
Figure RE-GDA0003758055760000072
finally, the optimal control problem is discretized into a nonlinear programming problem expressed as:
Figure RE-GDA0003758055760000073
Figure RE-GDA0003758055760000074
Figure RE-GDA0003758055760000075
Figure RE-GDA0003758055760000076
Figure RE-GDA0003758055760000077
Figure RE-GDA0003758055760000078
where φ is a boundary event constraint function, C 1 And C 2 For the equality and inequality constraints,
Figure RE-GDA0003758055760000079
is the state variable at the end of the kth trellis,
Figure RE-GDA00037580557600000710
is a state variable at the start of the (k + 1) th trellis;
step (63) of calculating the maximum relative error of the kth grid subinterval
Figure RE-GDA00037580557600000711
If it is
Figure RE-GDA00037580557600000712
Then the next step is carried out, otherwise, the step (65) is skipped;
let N of the kth grid k +1 samples are:
Figure RE-GDA00037580557600000713
the residuals of the state constraint equation and the process constraint equation at the sample points are expressed as:
Figure RE-GDA00037580557600000714
Figure RE-GDA00037580557600000715
wherein
Figure RE-GDA00037580557600000716
And
Figure RE-GDA00037580557600000717
is a residual error that is a function of the error,
Figure RE-GDA00037580557600000718
is the sampling point of the state variable and the control variable
Figure RE-GDA00037580557600000719
The value of (a);
taking the maximum residual error of the current interval as the maximum relative error:
Figure RE-GDA00037580557600000720
step (64): computing the relative curvature r of the current mesh (k) If r is (k) ≤r max Recalculating the order of the interpolation polynomial; if r (k) ＞r max If yes, subdividing the current grid and calculating the position of a new interpolation point; then, the steps (62) to (64) are carried out until
Figure RE-GDA00037580557600000721
Relative curvature r (k) The expression of (c) is:
Figure RE-GDA0003758055760000081
in the formula (I), the compound is shown in the specification,
Figure RE-GDA0003758055760000082
maximum and average curvatures for all LGR points in the kth grid;
curvature of each point
Figure RE-GDA0003758055760000083
The expression is as follows:
Figure RE-GDA0003758055760000084
in the formula (I), the compound is shown in the specification,
Figure RE-GDA0003758055760000085
and
Figure RE-GDA0003758055760000086
respectively representing the maximum relative error in the kth grid subinterval
Figure RE-GDA0003758055760000087
A first derivative function and a second derivative function of the state function at the corresponding LGR point;
if the order of the polynomial is increased, the increased number Δ N can be expressed as:
Figure RE-GDA0003758055760000088
wherein ceil (·) represents rounding up, B is a positive integer with a value greater than 0 and adjustable;
refining the grid, and then the number n of the subintervals into which the interval is divided again k The following equation is obtained:
Figure RE-GDA0003758055760000089
wherein λ is an integer adjustment factor having a value greater than 0;
step (65): let k be k +1, proceed step (64) and step (65) for the next grid, and go straightTo all grid sub-intervals
Figure RE-GDA00037580557600000810
If the requirements are met, outputting the result and exiting the program.
A glide extended-range guided projectile comprising an on-board computer for on-line planning of an ideal trajectory of the glide extended-range guided projectile using the method of any one of claims 1 to 8.
Compared with the prior art, the invention has the remarkable advantages that:
(1) compared with other ideal trajectory planning methods meeting the requirements of hitting targets, the method for planning the trajectory of each guided projectile on line through the missile-borne computer reduces the burden of a fire control system.
(2) Aiming at the problem that the guided projectile is interfered by gust and the like in an initial ascending section and deviates too far from an ideal trajectory before reaching the top point of the actual trajectory, compared with the situation that the control law corresponding to the ideal trajectory planned by a fire control system is not suitable for interference, and the target cannot be hit, the method provided by the invention is used for carrying out online planning in a boosting section, the interference changes the position and motion state information of the guided projectile during the online planning, so that the error between the ideal trajectory and the actual trajectory planned on the basis is small, the difficulty of tracking the ideal trajectory by the guided projectile is reduced, and the target can be hit accurately.
(3) Compared with the problems of sensitivity to initial values and difficulty in solving two-point edge values caused by an indirect method, the technical method provided by the invention depends on a pseudo-spectrum method, obtains discrete points by a Newton iteration method in a continuous time period, selects a proper polynomial order and then performs equation approximation, avoids two difficulties of the indirect method, and ensures that the pseudo-spectrum method has good convergence and is easy to calculate when solving an infinite smooth problem.
(4) Compared with the problem that the pseudo-spectrum method based on the p method cannot solve the discontinuous or non-smooth problem, the technical method provided by the invention updates the grid based on the self-adaptive hp method, so that the error and curvature of sampling points in the grid meet the requirements, the characteristic that the h method can capture the non-smoothness of the solution is kept, the exponential convergence speed of the p method during solution is also kept, and the self-adaptive hp method improves the calculation efficiency and the resolution precision of the non-linear programming problem.
Drawings
FIG. 1 is a diagram of the angular relationship between a projectile coordinate system and a velocity coordinate system established by the present invention.
Fig. 2 is a diagram of the angular relationship between the translational coordinate system and the ballistic coordinate system established by the invention.
Fig. 3 is a diagram of the angular relationship between the ballistic coordinate system and the velocity coordinate system established by the present invention.
Fig. 4 is a graph of the longitudinal flight of an ideal trajectory of an on-line planned guided projectile.
Fig. 5 is a lateral flight diagram of an ideal trajectory of an online-planned guided projectile.
Fig. 6 is a velocity profile of an ideal trajectory of an on-line planned guided projectile.
Fig. 7 is a ballistic dip graph of an ideal trajectory of an on-line planned guided projectile.
Fig. 8 is an overload graph of an ideal trajectory of an online planned guided projectile.
Fig. 9 is a graph of the change of the pitching rudder deflection angle of the ideal trajectory of the guided cannonball planned on line.
Fig. 10 is a diagram of the change of the yaw rudder deflection angle of the ideal trajectory of the guided projectile planned on line.
FIG. 11 is an ideal ballistic diagram of the guided projectile based on off-line planning and on-line planning during gust disturbance.
Detailed Description
The present invention is described in further detail below with reference to the attached drawing figures.
In order to make the objects, technical solutions and advantages of the present application more apparent, the present application is described in further detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the present application and are not intended to limit the present application.
The method is based on the optimal control theory, considers the initial launching deviation and the interference on the guided projectile before the trajectory is planned on line, can effectively plan the ideal trajectory meeting the requirements of a drop point, a drop speed and a drop angle on line by a self-adaptive pseudo-spectrum method, and provides an efficient and accurate calculation method for researching the on-line planning of the ideal trajectory of the gliding extended-range guided projectile. An on-line planning method for ideal trajectory of gliding extended range guided projectile specifically comprises the following steps:
step (1)' establishing ground Axyz, translation Oxyz and projectile body Ox 1 y 1 z 1 Trajectory Ox 2 y 2 z 2 Velocity Ox 3 y 3 z 3 The specific method of the coordinate system is as follows:
the ground coordinate system is fixedly connected with the ground, the origin A is on the center of mass of the projectile body during launching, the Ax axis is a horizontal line connecting the launching device and the target and points to the target, the Ay axis is upward along the vertical line, and the Az axis is vertical to the Ax and Ay;
o is the instantaneous center of mass of the projectile body, and Ox, Oy and Oz are respectively parallel to and in the same direction as Ax, Ay and Az; ox 1 The axis is coincident with the longitudinal axis and points to the warhead as positive, Oy 1 The axis lying in a longitudinal plane, with Ox 1 Vertical, with upward positive, Oz 1 Axis perpendicular to Ox 1 y 1 A plane;
Ox 2 the axis coinciding with the instantaneous velocity V, Oy, of the centre of mass of the guided projectile 2 The axis is located in the vertical plane containing V and the axis x 2 Vertical axis, Oz 2 Axis perpendicular to Ox 2 Axis and Oy 2 A shaft; ox 3 The axis coinciding with V, Oy 3 The axis lying in the longitudinal symmetry plane with Ox 3 Vertical axis, Oz 3 Axis perpendicular to Ox 3 y 3 And (4) a plane.
The angular relationships between the translational, projectile, ballistic and velocity coordinate systems are shown in fig. 1-3.
Step 2, establishing a 3dof guided projectile motion equation, wherein the specific method comprises the following steps:
neglecting the roll of the guided projectile, the kinetic equation set of the center of mass of the guided projectile in the ballistic coordinate system is as follows:
Figure RE-GDA0003758055760000101
the kinematic equation of the center of mass of the guided projectile in the ground coordinate system is as follows:
Figure RE-GDA0003758055760000102
the mass center mass change equation and the control relation equation are as follows:
Figure RE-GDA0003758055760000103
wherein m is the mass of the guided projectile,
Figure RE-GDA0003758055760000104
for the rate of change of mass, the boost period is not zero only during operation of the engine, the start and end times of the boost period being dependent on the guided projectile itself, m c The specific value is that the drug loading of the guided cannonball is related to the combustion time of the engine, V is the instantaneous speed of the mass center of the guided cannonball,
Figure RE-GDA0003758055760000105
is the speed change rate, g is the gravity acceleration, P is the average thrust of the engine, and exists only in the boosting period, X is the air resistance borne by the guided projectile, Y is the lift force, Z is the lateral force, and the attack angle alpha is the speed in the projectile coordinate system Ox 1 y 1 z 1 On projection and axis Ox 1 Angle of inclination, alpha c (t) is the value of the attack angle alpha changing with time under the action of the control law, and the sideslip angle beta is the speed and projectile coordinate system Ox 1 y 1 z 1 Angle of (b) of c (t) is the value of slip angle beta under the action of the control law as a function of time, f α 、f β Based on the relationship function of the pneumatic parameters, the attack angle and the sideslip angle related to the gravity center and the pressure center of the guided projectile and the pressure center of the control rudder and the deflection angles of the two rudders, which are obtained by the instantaneous balance hypothesis development, the guided projectile is developed to control the control rudder when flying to the top point of the trajectory, and the speed inclination angle gamma is controlled v Oy being a speed coordinate system 3 Ox of axis and trajectory coordinate system 2 y 2 The angle between the ballistic inclination angle theta is the angle between V and the horizontal plane,
Figure RE-GDA0003758055760000106
is the rate of change of ballistic inclination angle, ballistic declination angle
Figure RE-GDA0003758055760000107
Is the angle between the projection of V on the horizontal plane and the Ox axis,
Figure RE-GDA0003758055760000108
is the rate of change of ballistic inclination, (x, y, z) is the instantaneous position of the center of mass of the guided projectile in the ground coordinate system,
Figure RE-GDA0003758055760000109
the position change rate of the center of mass of the guided projectile on the x axis, the y axis and the z axis of a ground coordinate system is obtained.
And (3) performing a wind tunnel test to obtain pneumatic parameters, wherein the specific method comprises the following steps:
the gliding extended-range guided projectile model is placed in a wind tunnel, and according to the action of gas flow and the guided projectile model, a plurality of groups of data of zero lift resistance coefficients, lift coefficient derivatives, pitch rudder and yaw rudder zero lift resistance coefficients, pitch rudder and yaw rudder lift coefficient derivatives, missile wing combination pressure center positions and control rudder pressure center positions corresponding to different Mach numbers are obtained, and an induced resistance coefficient and a control rudder induced resistance coefficient are selected.
The step (4) is specifically as follows:
in many cases of engineering application, a range is required to be specified, so that the shortest flight time is taken as a requirement, the flight time t of the guided projectile is selected as a final value type performance index, and the expression is as follows:
J＝t f
wherein J represents a performance index, t f Representing the time required to fly-off the ideal trajectory.
Setting the initial launching speed of the guided projectile to be 790m/s, the ballistic inclination angle to be 55 degrees, setting the position to be (0km,0km,0km), and setting the maximum value and the minimum value of a state variable and a control variable in the flight process to be process constraints, wherein the embodiment requires that the angles of an attack angle, a side slip angle and a rudder deflection angle in the flight process cannot exceed 6 degrees, the angle of the rudder deflection angle changing every second cannot exceed 2 degrees, the overload absolute value cannot exceed 3, setting a terminal parameter value of the guided projectile to be terminal constraints, and setting the drop angle of the guided projectile to be not less than 65 degrees or a specified drop angle when the guided projectile falls to the ground and setting the drop speed to be not less than 200 m/s.
The step (5) is specifically as follows:
beidou satellite represents the position P point (x) of the guided cannonball based on WGS-84 coordinate system ε ,y ε ,z ε ) Conversion to position data in a geodetic coordinate system
Figure RE-GDA0003758055760000111
Figure RE-GDA0003758055760000112
Wherein a represents the major semi-axis of the WGS-84 ellipsoid and has a value of 6378137.0 m; e is the first eccentricity of WGS-84 with a square of 0.0066943799; and N is the curvature radius of the ellipsoidal unitary fourth of twelve earthly branches. The origin w of the WGS-84 coordinate system is the earth centroid, the wz axis points to the north pole of an ellipsoid, the wx axis points to the intersection point of the equatorial plane and the Greenwich meridian plane, and the wy axis in the equatorial plane and xwz form a right-hand system coordinate system; the short semi-axis of the ellipsoid of the geodetic coordinate system is coincident with the rotation axis wz of the earth, the longitude lambda is the included angle between the meridian plane of the ellipsoid where the point P is located and the meridian plane of Greenwich mean, and the latitude
Figure RE-GDA0003758055760000113
Is the angle between the normal of the ellipsoid passing through point P and the equatorial plane, and the height h is the distance from the normal of the ellipsoid to the ellipsoid along point P.
The coordinates of the guided projectile in the geodetic coordinate system are converted into coordinates (x) in the north heaven coordinate system n ,y n ,z n )：
Figure RE-GDA0003758055760000114
In the formula (x) ε0 ,y ε0 ,z ε0 ) The position of the origin of a north celestial east coordinate system in a WGS-84 coordinate system is shown, the origin n of the north celestial east coordinate system is a transmitter point, xnz is a ground horizontal plane where the origin is located, an nx axis points to the north of the earth, an nz axis points to the east of the earth, and a ny axis is perpendicular to a xnz plane and points to the air.
The north-east coordinate system to the shooting coordinate system (x, y, z) is:
Figure RE-GDA0003758055760000115
in the formula, α is the shot direction under the north-heaven-east coordinate system, and the shooting coordinate system is the ground coordinate system in the step 1.
The step (6) is specifically as follows:
obtaining the motion state equation constraint of optimal control according to the dynamics, the kinematics equation set and the mass change equation in the step (2)
Figure RE-GDA0003758055760000121
And (5) obtaining the performance index J and constraint conditions according to the step (4), and finally determining the description model of the optimal control problem as follows:
Figure RE-GDA0003758055760000122
Figure RE-GDA0003758055760000123
Figure RE-GDA0003758055760000124
Figure RE-GDA0003758055760000125
Figure RE-GDA0003758055760000126
j stands for Performance index,. phi.is the terminal Performance index, t f Indicating the time required to fly-off the ideal trajectory,
Figure RE-GDA0003758055760000127
respectively representing vectors formed by all state variables and all control variables, 4, the control variables are derivatives of two rudder deflection angles with respect to time, the initial value of the rudder deflection angle is 0 degree, the numerical value of the rudder deflection angle is obtained through integration,
Figure RE-GDA0003758055760000128
to represent
Figure RE-GDA0003758055760000129
Derivative with respect to time, f is a mapping by which
Figure RE-GDA00037580557600001210
And
Figure RE-GDA00037580557600001211
to represent
Figure RE-GDA00037580557600001212
q is
Figure RE-GDA00037580557600001213
And
Figure RE-GDA00037580557600001214
as a function of (a) or (b),
Figure RE-GDA00037580557600001215
a state vector representing the end of the trajectory,
Figure RE-GDA00037580557600001216
a state vector representing the starting point of the trajectory, phi being a boundary event constraint function, comprising a state vector x (t) for the starting point 0 ) And an end point state vector x (t) f ) The constraint of (2); c. C 1 And c 2 Equality and inequality constraints.
The step (6) of reconstructing grids by using a self-adaptive hp method, optimizing the number of grids and an interpolation polynomial order, and converting the problem into a nonlinear programming problem by using a pseudo-spectral method to disperse an optimal control problem specifically comprises the following steps:
step (61): carrying out initial grid division; the maximum allowable error epsilon is set to 0.0001, and the maximum allowable curvature r max For 1, the number of initialization grids is 10 and the order of the lagrange interpolation polynomial is 4.
Will time interval t 0 ,t f ]Divided into K grid subintervals [ t ] k-1 ,t k ]K1, …, K, a new time variable τ e [ -1,1 ] is introduced]The original variable t is converted into a new variable tau through the following formula, so that the single-interval nonlinear optimal control problem can be converted into the multi-interval nonlinear optimal control problem
Figure RE-GDA00037580557600001217
Step (62): calculating the position of an LGR interpolation point, an LGR weight and a Radau pseudo-spectral differential matrix in each current grid subinterval, approximately expressing an optimal control model by combining a Lagrange interpolation polynomial, and converting the optimal control model into a nonlinear programming problem;
let the kth grid subinterval [ t ] k-1 ,t k ]The state variables and control variables in are denoted X k (τ) and U k (τ)，N k +1 discrete points
Figure RE-GDA00037580557600001218
And end point
Figure RE-GDA00037580557600001219
Form the interval tau epsilon-1, 1]Inner N k +2 Lagrangian interpolation points, wherein
Figure RE-GDA00037580557600001220
Is provided with
Figure RE-GDA00037580557600001221
Is N k A +2 order lagrange interpolation basis function,
Figure RE-GDA00037580557600001222
is N k The +1 order Lagrange interpolation basis function has the expression:
Figure RE-GDA0003758055760000131
Figure RE-GDA0003758055760000132
state variable X (k) In N k +2 discrete points of use
Figure RE-GDA0003758055760000133
The approximation is:
Figure RE-GDA0003758055760000134
in the formula (I), the compound is shown in the specification,
Figure RE-GDA0003758055760000135
is a state variable X k At discrete points
Figure RE-GDA0003758055760000136
Taking the value of (1);
control variable U (k) In N k +1 discrete point application
Figure RE-GDA0003758055760000137
The approximation is:
Figure RE-GDA0003758055760000138
in the formula (I), the compound is shown in the specification,
Figure RE-GDA0003758055760000139
is a state variable U k At discrete points
Figure RE-GDA00037580557600001310
Taking the value of (1);
for state variable X (k) The derivative to get the approximate equation of state is:
Figure RE-GDA00037580557600001311
order to
Figure RE-GDA00037580557600001312
The specific expression of a Radau pseudo-spectral differential approximation matrix representing the kth subinterval is as follows:
Figure RE-GDA00037580557600001313
Figure RE-GDA00037580557600001314
in the formula, g (k) (τ) is P n Function of (τ), P n (τ) is a Legendre orthogonal polynomial, and the differential expression is:
Figure RE-GDA00037580557600001315
then N is k The state equation constraints at +1 interpolation points are converted into algebraic equation constraints expressed as:
Figure RE-GDA0003758055760000141
in the formula F s (k) With respect to state variables and control variablesA vector of an equation of state function;
the LGR integral is used for approximating an integral term in a Boerch's performance index, and the LGR integral is approximated as:
Figure RE-GDA0003758055760000142
in the formula (I), the compound is shown in the specification,
Figure RE-GDA0003758055760000143
is a function of the state vector and the control vector,
Figure RE-GDA0003758055760000144
as discrete points of the kth sub-interval
Figure RE-GDA0003758055760000145
The corresponding LGR weights are in the specific form:
Figure RE-GDA0003758055760000146
Figure RE-GDA0003758055760000147
finally, the optimal control problem is discretized into a nonlinear programming problem expressed as:
Figure RE-GDA0003758055760000148
Figure RE-GDA0003758055760000149
Figure RE-GDA00037580557600001410
Figure RE-GDA00037580557600001411
Figure RE-GDA00037580557600001412
Figure RE-GDA00037580557600001413
where φ is a boundary event constraint function, C 1 And C 2 For the equality and inequality constraints,
Figure RE-GDA00037580557600001414
is the state vector of the end point of the kth trellis,
Figure RE-GDA00037580557600001415
a state vector that is the starting point of the (k + 1) th trellis;
step (63) of calculating the maximum relative error of the kth grid subinterval
Figure RE-GDA00037580557600001416
If it is
Figure RE-GDA00037580557600001417
Proceeding to the next step, otherwise jumping to the step (65);
let N of the kth grid k +1 sample points are:
Figure RE-GDA00037580557600001418
the residuals of the state constraint equation and the process constraint equation at the sample points are expressed as:
Figure RE-GDA00037580557600001419
Figure RE-GDA0003758055760000151
wherein
Figure RE-GDA0003758055760000152
And
Figure RE-GDA0003758055760000153
is a residual error that is a function of the error,
Figure RE-GDA0003758055760000154
is the sampling point of the state variable and the control variable
Figure RE-GDA0003758055760000155
The value of (a);
taking the maximum residual error of the current interval as the maximum relative error:
Figure RE-GDA0003758055760000156
step (64): computing the relative curvature r of the current mesh (k) If r is (k) ≤r max Recalculating the order of the interpolation polynomial; if r (k) ＞r max If so, subdividing the current grid and calculating the position of a new interpolation point; then, the steps (62) to (64) are carried out until
Figure RE-GDA0003758055760000157
Relative curvature r (k) The expression of (a) is:
Figure RE-GDA0003758055760000158
in the formula (I), the compound is shown in the specification,
Figure RE-GDA0003758055760000159
is the maximum curvature sum of all LGR points in the kth gridAn average curvature;
curvature of each point
Figure RE-GDA00037580557600001510
The expression is as follows:
Figure RE-GDA00037580557600001511
in the formula (I), the compound is shown in the specification,
Figure RE-GDA00037580557600001512
and
Figure RE-GDA00037580557600001513
respectively representing the maximum relative error in the kth grid subinterval
Figure RE-GDA00037580557600001514
A first derivative function and a second derivative function of the state function at the corresponding LGR point;
if the order of the polynomial is increased, the amount of increase Δ N can be expressed as:
Figure RE-GDA00037580557600001515
wherein ceil (·) represents rounding-up, B is a positive integer with a value greater than 0 and adjustable;
refining the grid, and then the number n of the subintervals into which the interval is divided again k The following equation is obtained:
Figure RE-GDA00037580557600001516
wherein λ is an integer adjustment factor having a value greater than 0;
step (65): and k is k +1, and the step (64) and the step (65) are carried out on the next grid until all grid subintervals are reached
Figure RE-GDA00037580557600001517
If the requirements are met, outputting the result and exiting the program.
A certain type of guided cannonball is made to strike a target (65km, 0m and 0m) at a temperature of-65 degrees, the guided cannonball is made to be launched at an initial speed of 790m/s and a shooting angle of 55 degrees, the time of rising to a ballistic vertex is 57.8s when the guided cannonball is launched at a time of 0, and a missile-borne computer is used for shooting at t 0 The parameters of the guided projectile are received at time 20s and the trajectory planning begins. When the flying time of the guided cannonball is not interfered for 20s, the parameter information at the moment is as follows: position (8585m,10170m,0m), velocity 622.94m/s, ballistic inclination angle 44.114 °, ballistic declination angle 0 °. Suppose that the projectile goes from fire to t 0 Unknown disturbance is experienced at the moment, deviation occurs in the state at the moment, four groups of ballistic parameters are arbitrarily selected for simulation, and the parameters are shown in table 1, wherein the unit of distance is (m), the unit of speed is (m/s), and the unit of angle is (DEG). The simulation results are shown in fig. 4-10. Fig. 4-10 show the trajectory and motion parameter curves of the planned trajectory on line, respectively, and it can be seen that the trajectory planned on line with the four sets of parameters in table 1 as the starting points all fulfill the hit objective on the premise of satisfying the constraints of the variables.
Table 1 randomly selected trajectory planning point parameters
Figure RE-GDA0003758055760000161
The projectile was fired at v 790m/s and θ 55 °, hitting the target at x 40km with a drop angle of-65 °. And adding gusts with the direction of 0-50 m/s and the opposite bullet speed at 18-28 s, and drawing the trajectory by using the self-adaptive hp fairway pseudospectral law. The first set of tests used an off-line planned control law, and the second set planned the trajectory on-line at t-30 s based on the technical method of the invention. FIG. 11 shows two sets of trajectories, where the first set has a landing point of 35.41km, a too large error in distance from the target, and a hit failure, and the second set has fulfilled the landing angle and landing point requirements.
A glide extended range guided projectile comprising an on-board computer for on-line planning of an ideal trajectory of the glide extended range guided projectile using the method of any one of claims 1 to 8.

Claims (9)

1. An online planning method for an ideal trajectory of a gliding extended range guided projectile is characterized by comprising the following steps:
step (1): establishing ground Axyz, translation Oxyz and projectile body Ox 1 y 1 z 1 Trajectory Ox 2 y 2 z 2 And velocity Ox 3 y 3 z 3 A coordinate system;
step (2): establishing a motion equation of the 3dof guided projectile;
and (3): performing a wind tunnel test to obtain pneumatic parameters;
and (4): establishing a performance index of a trajectory planning problem according to a hitting requirement, and establishing a constraint condition according to the launching facility capacity of a gun platform, the strength limit of a bullet body and the strength requirement of a hitting target;
and (5): planning an ideal trajectory at the initial ascending section or the boosting section after launching, acquiring and processing satellite data at the moment, and substituting the data as a starting point into a guided projectile motion equation;
and (6): reconstructing grids by using a self-adaptive hp method, optimizing the number of the grids and the polynomial order of interpolation, converting the problem of discrete optimal control by using a pseudo-spectrum method into a nonlinear programming problem, and solving by using a nonlinear programming solver to obtain an ideal ballistic curve, a velocity curve, a ballistic dip angle curve, an overload curve and a rudder deflection angle curve of the guided projectile centroid motion.
2. The method of claim 1, wherein step (1) "establish ground Axyz, translation Oxyz, projectile Ox 1 y 1 z 1 Trajectory Ox 2 y 2 z 2 Velocity Ox 3 y 3 z 3 The specific method of the coordinate system is as follows:
the ground coordinate system is fixedly connected with the ground, the origin A is on the center of mass of the projectile body during launching, the Ax axis is a horizontal line connecting the launching device and the target and points to the target, the Ay axis is upward along the vertical line, and the Az axis is vertical to the Ax and Ay;
o is instantaneous center of mass of the projectile, OxOy and Oz are respectively parallel and equidirectional with Ax, Ay and Az; ox 1 The axis is coincident with the longitudinal axis and points to the warhead as positive, Oy 1 The axis lying in a longitudinal plane, with Ox 1 Vertical, with upward positive, Oz 1 Axis perpendicular to Ox 1 y 1 A plane;
Ox 2 the instantaneous velocity V of the axis coincident with the center of mass of the guided projectile, Oy 2 The axis is located in the vertical plane containing V and the axis x 2 Vertical axis, Oz 2 Axis perpendicular to Ox 2 Axis and Oy 2 A shaft; ox (oxygen) 3 The axis coinciding with V, Oy 3 The axis lying in the longitudinal symmetry plane with Ox 3 Vertical axis, Oz 3 Axis perpendicular to Ox 3 y 3 And (4) a plane.
3. The method of claim 2, wherein the step (2) of establishing the 3dof guided projectile motion equation comprises the following steps:
neglecting the roll of the guided projectile, and the dynamic equation of the center of mass of the guided projectile in a ballistic coordinate system is as follows:
Figure RE-FDA0003758055750000011
the kinematic equation of the center of mass of the guided projectile in the ground coordinate system is as follows:
Figure RE-FDA0003758055750000012
the mass center mass change equation and the control relation equation are as follows:
Figure RE-FDA0003758055750000021
wherein m is the mass of the guided projectile,
Figure RE-FDA0003758055750000022
for rate of change of mass, only during engine operationI.e. the boosting period is not zero, the starting time and the ending time of the boosting period depend on the guided cartridge per se, m c The specific numerical value is that the medicine loading of the guided cannonball is related to the combustion time of the engine, V is the instantaneous speed of the center of mass of the guided cannonball,
Figure RE-FDA0003758055750000023
the speed change rate is g, the gravity acceleration is g, the average thrust of the engine is P, the engine exists only in the boosting period, X is the air resistance borne by the guided projectile, Y is the lift force, Z is the lateral force, and the attack angle alpha is the speed in a projectile coordinate system Ox 1 y 1 z 1 On projection and axis Ox 1 Angle of inclination, alpha c (t) is the value of the attack angle alpha changing with time under the action of the control law, and the sideslip angle beta is the speed and projectile coordinate system Ox 1 y 1 z 1 Angle of (b) of c (t) is the value of the slip angle beta under the action of the control law as a function of time, f α 、f β Based on a relation function of an attack angle, a side slip angle and two rudder deflection angles which are obtained by the instantaneous balance hypothesis expansion and are related to pneumatic parameters, the gravity center of the guided projectile, the pressure center and the pressure center of the control rudder, the guided projectile is expanded to the top point of a trajectory and then the control rudder is expanded for control, and the speed inclination angle gamma is v Is the Oy of the velocity coordinate system 3 Ox of axis and trajectory coordinate system 2 y 2 The angle between the ballistic inclination angle theta is the angle between V and the horizontal plane,
Figure RE-FDA0003758055750000024
is the rate of change of ballistic inclination angle, ballistic declination angle
Figure RE-FDA0003758055750000025
Is the included angle between the projection of V on the horizontal plane and the Ox axis,
Figure RE-FDA0003758055750000026
is the rate of change of ballistic inclination, (x, y, z) is the instantaneous position of the center of mass of the guided projectile in the ground coordinate system,
Figure RE-FDA0003758055750000027
the position change rate of the center of mass of the guided projectile on the x axis, the y axis and the z axis of a ground coordinate system.
4. The method according to claim 3, wherein the wind tunnel test in step (3) is performed to obtain the pneumatic parameters by:
the gliding extended-range guided projectile model is placed in a wind tunnel, a plurality of groups of data of zero lift resistance coefficients, lift coefficient derivatives, pitch rudder and yaw rudder zero lift resistance coefficients, pitch rudder and yaw rudder lift coefficient derivatives, missile wing combination pressure center positions and control rudder pressure center positions corresponding to different Mach numbers are obtained according to the action of gas flow and the guided projectile model, and an induced resistance coefficient and a control rudder induced resistance coefficient are selected.
5. The method according to claim 4, wherein step (4) is specifically:
when the range is maximized, selecting the x-axis coordinate of the flight end point of the guided projectile as a final value type performance index, wherein the expression is as follows:
J＝-x(t f )
wherein J represents a performance index, t f Represents the time required to fly off the ideal trajectory, x (t) f ) Representing the coordinates of the ballistic endpoint on the axis of the ground coordinate system Ax;
setting a starting point parameter value of the guided projectile as a starting point constraint, setting the maximum value and the minimum value of a state variable and a control variable in the flight process as a process constraint, and setting the falling angle and the falling speed of the guided projectile when the guided projectile falls to the ground to be not less than a terminal point constraint set by a user.
6. The method according to claim 5, wherein step (5) is embodied as:
beidou satellite represents the position P point (x) of the guided cannonball based on WGS-84 coordinate system ε ,y ε ,z ε ) Conversion to position data in a geodetic coordinate system
Figure RE-FDA0003758055750000028
Figure RE-FDA0003758055750000031
Wherein a represents the major semi-axis of the WGS-84 ellipsoid and has a value of 6378137.0 m; e is the first eccentricity of WGS-84, squared 0.0066943799; n is half of the curvature of the ellipsoidal unitary ring; the origin w of the WGS-84 coordinate system is the earth centroid, the wz axis points to the north pole of an ellipsoid, the wx axis points to the intersection point of the equatorial plane and the Greenwich meridian plane, and the wy axis in the equatorial plane and xwz form a right-hand system coordinate system; the short semi-axis of the ellipsoid of the geodetic coordinate system is coincident with the rotation axis wz of the earth, the longitude lambda is the included angle between the meridian plane of the ellipsoid where the point P is located and the meridian plane of Greenwich mean, and the latitude
Figure RE-FDA0003758055750000032
Is the included angle between the normal line of the ellipsoid passing through the point P and the equatorial plane, and the height h is the distance from the normal line of the ellipsoid to the ellipsoid along the point P;
the coordinates of the guided projectile in the geodetic coordinate system are converted into coordinates (x) in the north-celestial-east coordinate system n ,y n ,z n )：
Figure RE-FDA0003758055750000033
In the formula (x) ε0 ,y ε0 ,z ε0 ) The method is characterized in that the original point of a north heaven coordinate system is in a WGS-84 coordinate system, the original point n of the north heaven coordinate system is a transmitter point, xnz is a ground horizontal plane where the original point is located, an nx axis points to the north of the earth, an nz axis points to the east of the earth, and a ny axis is perpendicular to a xnz plane and points to the air;
the north-east coordinate system to the shooting coordinate system (x, y, z) is:
Figure RE-FDA0003758055750000034
in the formula, α is the shot direction in the north-east-heaven coordinate system, and the shooting coordinate system is the ground coordinate system.
7. The method according to claim 6, characterized in that step (6) is embodied as:
obtaining the motion state equation constraint of optimal control according to the dynamics and kinematics equation set and the mass change equation in step (3) and step (2)
Figure RE-FDA0003758055750000035
And (5) obtaining the performance index J and constraint conditions according to the step (4), and finally determining the description model of the optimal control problem as follows:
Figure RE-FDA0003758055750000041
Figure RE-FDA0003758055750000042
Figure RE-FDA0003758055750000043
Figure RE-FDA0003758055750000044
Figure RE-FDA0003758055750000045
wherein J represents a performance index, Φ is a final value type performance index, t f Indicating the time required to fly-off the ideal trajectory,
Figure RE-FDA0003758055750000046
respectively representing all state variables, instituteThe vector formed by control variables is provided, the control variables are derivatives of two rudder deflection angles with respect to time, the initial value of the rudder deflection angle is 0 degree, the numerical value of the rudder deflection angle is obtained through integration,
Figure RE-FDA0003758055750000047
to represent
Figure RE-FDA0003758055750000048
Derivative with respect to time, f is a mapping by which
Figure RE-FDA0003758055750000049
And
Figure RE-FDA00037580557500000410
represent
Figure RE-FDA00037580557500000411
q is
Figure RE-FDA00037580557500000412
And
Figure RE-FDA00037580557500000413
as a function of (a) or (b),
Figure RE-FDA00037580557500000414
a state vector representing the end of the trajectory,
Figure RE-FDA00037580557500000415
a state vector representing the starting point of the trajectory, #, being a boundary event constraint function, including a state vector x (t) for the starting point 0 ) And an end point state vector x (t) f ) The constraint of (2); c. C 1 And c 2 Are equality and inequality constraints.
8. The method according to claim 7, wherein the step (6) of reconstructing the grid by using the adaptive hp method, optimizing the grid number and the interpolation polynomial order, and converting the discrete optimal control problem into the nonlinear programming problem by using the pseudo-spectrum method specifically comprises the following steps:
step (61): carrying out initial grid division; the maximum allowable error epsilon and the maximum allowable curvature r are set max The number of initial grids and the order of the lagrange interpolation polynomial;
will time interval t 0 ,t f ]Divided into K grid subintervals [ t ] k-1 ,t k ]K1, …, K, a new time variable τ e [ -1,1 ] is introduced]The original variable t is converted into a new variable tau through the following formula, so that the single-interval nonlinear optimal control problem can be converted into a multi-interval nonlinear optimal control problem:
Figure RE-FDA00037580557500000416
step (62): calculating the position of an LGR interpolation point, an LGR weight and a Radau pseudo-spectral differential matrix in each current grid subinterval, approximately expressing an optimal control model by combining a Lagrange interpolation polynomial, and converting the optimal control model into a nonlinear programming problem;
let the kth grid subinterval [ t ] k-1 ,t k ]The state variables and control variables in are denoted X k (τ) and U k (τ)，N k +1 discrete points
Figure RE-FDA00037580557500000417
And end point
Figure RE-FDA00037580557500000418
Form the interval tau epsilon-1, 1]Inner N k +2 Lagrangian interpolation points, wherein
Figure RE-FDA00037580557500000419
Is provided with
Figure RE-FDA00037580557500000420
Is N k +2 Lagrangian orderThe basis functions of the interpolation are then interpolated,
Figure RE-FDA00037580557500000421
is N k The +1 order Lagrange interpolation basis function has the expression:
Figure RE-FDA0003758055750000051
Figure RE-FDA0003758055750000052
state variable X (k) In N k +2 discrete point applications
Figure RE-FDA0003758055750000053
The approximation is:
Figure RE-FDA0003758055750000054
in the formula (I), the compound is shown in the specification,
Figure RE-FDA0003758055750000055
is a state variable X k At discrete points
Figure RE-FDA0003758055750000056
Taking the value of (A);
control variable U (k) In N k +1 discrete point application
Figure RE-FDA0003758055750000057
The approximation is:
Figure RE-FDA0003758055750000058
in the formula (I), the compound is shown in the specification,
Figure RE-FDA0003758055750000059
is a state variable U k At discrete points
Figure RE-FDA00037580557500000510
Taking the value of (A);
for state variable X (k) The derivative to get the approximate equation of state is:
Figure RE-FDA00037580557500000511
order to
Figure RE-FDA00037580557500000512
The specific expression of a Radau pseudo-spectral differential approximation matrix representing the kth subinterval is as follows:
Figure RE-FDA00037580557500000513
Figure RE-FDA00037580557500000514
in the formula, g (k) (τ) is P n Function of (τ), P n (τ) is a Legendre orthogonal polynomial, and the differential expression is:
Figure RE-FDA00037580557500000515
then N is k The state equation constraints on +1 interpolation points are converted into algebraic equation constraints expressed as:
Figure RE-FDA0003758055750000061
in the formula F s (k) Is a vector of a function of the state equation with respect to the state variables and the control variables;
the LGR integral is used for approximating an integral term in a Boercus-type performance index, and the approximation is as follows:
Figure RE-FDA0003758055750000062
in the formula (I), the compound is shown in the specification,
Figure RE-FDA0003758055750000063
is a function of the state vector and the control vector,
Figure RE-FDA0003758055750000064
as discrete points of the kth sub-interval
Figure RE-FDA0003758055750000065
The corresponding LGR weights are in the specific form:
Figure RE-FDA0003758055750000066
Figure RE-FDA0003758055750000067
finally, the optimal control problem is discretized into a nonlinear programming problem, expressed as:
Figure RE-FDA0003758055750000068
Figure RE-FDA0003758055750000069
Figure RE-FDA00037580557500000610
Figure RE-FDA00037580557500000611
Figure RE-FDA00037580557500000612
Figure RE-FDA00037580557500000613
where φ is a boundary event constraint function, C 1 And C 2 For the equality and inequality constraints,
Figure RE-FDA00037580557500000614
is the state variable at the end of the kth trellis,
Figure RE-FDA00037580557500000615
is a state variable at the start of the (k + 1) th trellis;
step (63) of calculating the maximum relative error of the kth grid subinterval
Figure RE-FDA00037580557500000616
If it is
Figure RE-FDA00037580557500000617
Then the next step is carried out, otherwise, the step (65) is skipped;
let N of the kth grid k +1 samples are:
Figure RE-FDA00037580557500000618
the residuals of the state constraint equation and the process constraint equation at the sample points are expressed as:
Figure RE-FDA0003758055750000071
Figure RE-FDA0003758055750000072
wherein
Figure RE-FDA0003758055750000073
And
Figure RE-FDA0003758055750000074
is a residual error that is a function of the error,
Figure RE-FDA0003758055750000075
is the sampling point of the state variable and the control variable
Figure RE-FDA0003758055750000076
The value of (a);
taking the maximum residual error of the current interval as the maximum relative error:
Figure RE-FDA0003758055750000077
step (64): computing the relative curvature r of the current mesh (k) If r is (k) ≤r max Recalculating the order of the interpolation polynomial; if r is (k) ＞r max If so, subdividing the current grid and calculating the position of a new interpolation point; then, the steps (62) to (64) are carried out until
Figure RE-FDA0003758055750000078
Relative to each otherCurvature r (k) The expression of (a) is:
Figure RE-FDA0003758055750000079
in the formula (I), the compound is shown in the specification,
Figure RE-FDA00037580557500000710
maximum and average curvatures for all LGR points in the kth grid;
curvature of each point
Figure RE-FDA00037580557500000711
The expression is as follows:
Figure RE-FDA00037580557500000712
in the formula (I), the compound is shown in the specification,
Figure RE-FDA00037580557500000713
and
Figure RE-FDA00037580557500000714
respectively representing the maximum relative error in the kth grid subinterval
Figure RE-FDA00037580557500000715
A first derivative function and a second derivative function of the state function at the corresponding LGR point;
if the order of the polynomial is increased, the increased number Δ N can be expressed as:
Figure RE-FDA00037580557500000716
wherein ceil (·) represents rounding-up, B is a positive integer with a value greater than 0 and adjustable;
refining the grid, subdividing the interval into sub-divisionsNumber of intervals n k The following equation is obtained:
Figure RE-FDA00037580557500000717
wherein λ is an integer adjustment factor having a value greater than 0;
step (65): let k be k +1, perform step (64) and step (65) for the next grid until all grid subintervals
Figure RE-FDA00037580557500000718
If the requirements are met, outputting the result and exiting the program.
9. A gliding extended-range guided projectile comprising an on-board computer for on-line planning of an ideal trajectory of the gliding extended-range guided projectile using the method of any one of claims 1 to 8.
CN202210212623.XA 2022-03-05 2022-03-05 Online planning method for ideal trajectory of gliding Cheng Zhidao shell Active CN114935277B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210212623.XA CN114935277B (en) 2022-03-05 2022-03-05 Online planning method for ideal trajectory of gliding Cheng Zhidao shell

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210212623.XA CN114935277B (en) 2022-03-05 2022-03-05 Online planning method for ideal trajectory of gliding Cheng Zhidao shell

Publications (2)

Publication Number Publication Date
CN114935277A true CN114935277A (en) 2022-08-23
CN114935277B CN114935277B (en) 2023-08-04

Family

ID=82862015

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210212623.XA Active CN114935277B (en) 2022-03-05 2022-03-05 Online planning method for ideal trajectory of gliding Cheng Zhidao shell

Country Status (1)

Country Link
CN (1) CN114935277B (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116362163A (en) * 2023-06-01 2023-06-30 西安现代控制技术研究所 Nonsingular multi-constraint trajectory rapid optimization method

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB1118663A (en) * 1965-06-03 1968-07-03 North American Aviation Inc Inertial navigation system error correcting methods
WO2019024303A1 (en) * 2017-08-02 2019-02-07 华南理工大学 Stable flight control method for multi-rotor unmanned aerial vehicle based on finite-time neurodynamics
CN111442697A (en) * 2020-02-07 2020-07-24 北京航空航天大学 Over-emphasis guidance method and trajectory shaping guidance method based on pseudo-spectrum correction
CN111859527A (en) * 2020-06-04 2020-10-30 中国人民解放军国防科技大学 Online planning method for whole-course trajectory of boosting gliding missile
WO2021036778A1 (en) * 2019-08-29 2021-03-04 大连理工大学 Method for directly planning reentry trajectory in altitude-velocity profile
CN112947534A (en) * 2021-04-23 2021-06-11 成都凯天通导科技有限公司 Adaptive pseudo-spectral method trajectory optimization method for depression section of hypersonic aircraft

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB1118663A (en) * 1965-06-03 1968-07-03 North American Aviation Inc Inertial navigation system error correcting methods
WO2019024303A1 (en) * 2017-08-02 2019-02-07 华南理工大学 Stable flight control method for multi-rotor unmanned aerial vehicle based on finite-time neurodynamics
WO2021036778A1 (en) * 2019-08-29 2021-03-04 大连理工大学 Method for directly planning reentry trajectory in altitude-velocity profile
CN111442697A (en) * 2020-02-07 2020-07-24 北京航空航天大学 Over-emphasis guidance method and trajectory shaping guidance method based on pseudo-spectrum correction
CN111859527A (en) * 2020-06-04 2020-10-30 中国人民解放军国防科技大学 Online planning method for whole-course trajectory of boosting gliding missile
CN112947534A (en) * 2021-04-23 2021-06-11 成都凯天通导科技有限公司 Adaptive pseudo-spectral method trajectory optimization method for depression section of hypersonic aircraft

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116362163A (en) * 2023-06-01 2023-06-30 西安现代控制技术研究所 Nonsingular multi-constraint trajectory rapid optimization method
CN116362163B (en) * 2023-06-01 2023-09-01 西安现代控制技术研究所 Multi-constraint trajectory rapid optimization method

Also Published As

Publication number Publication date
CN114935277B (en) 2023-08-04

Similar Documents

Publication Publication Date Title
CN111591470B (en) Aircraft precise soft landing closed-loop guidance method adapting to thrust adjustable mode
CN109933847B (en) Improved active segment trajectory estimation algorithm
CN111721291A (en) Engineering algorithm for strapdown inertial navigation under launching system
CN107367941B (en) Method for observing attack angle of hypersonic aircraft
CN111444603B (en) Method for rapidly planning shortest time off-orbit trajectory of recoverable spacecraft
CN113642122B (en) Single-shot-table-based remote interception emission data acquisition method and system
CN109737812B (en) Air-to-ground guided weapon side attack method and device
Baranowski Equations of motion of a spin-stabilized projectile for flight stability testing
CN108073742B (en) Method for estimating flight state of intercepted missile tail section based on improved particle filter algorithm
CN113602532A (en) Solid carrier rocket in-orbit correction method
CN108646554B (en) Aircraft rapid anti-interference longitudinal guidance method based on designated performance
CN114611416B (en) LS-SVM modeling method for nonlinear unsteady aerodynamic characteristics of missile
CN114935277A (en) Online planning method for ideal trajectory of gliding extended-range guided projectile
Yang et al. An aerodynamic shape optimization study to maximize the range of a guided missile
CN111649734B (en) Particle swarm algorithm-based strapdown seeker target positioning method
CN115574666A (en) Height setting method for grazing cruise target
CN114440707A (en) Aircraft guidance method, device and system for top and side cooperative interception
CN114637304A (en) Scouting and shooting weapon system and follow-up tracking control method
Fariz et al. Missile Initial Engagement Determination and Terminal Phase Guidance
CN102506864B (en) Method for approximately outputting quaternion numbers with arbitrary step size in orthogonal series during extreme flight of aircraft
CN116483109B (en) Sliding mode control-based integrated aircraft guidance control method and system
CN116227340A (en) Multiple maneuver burst prevention method in zero range line direction
de Celis et al. A neural network for sensor hybridization in rocket guidance
He et al. A New Intelligent Recognition Method for Maneuver Modes of Re-entry Gliding Vehicle
CN117706921A (en) High-precision continuous thrust approaching track planning method for spacecraft

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant