CN114935277B - Online planning method for ideal trajectory of gliding Cheng Zhidao shell - Google Patents
Online planning method for ideal trajectory of gliding Cheng Zhidao shell Download PDFInfo
- Publication number
- CN114935277B CN114935277B CN202210212623.XA CN202210212623A CN114935277B CN 114935277 B CN114935277 B CN 114935277B CN 202210212623 A CN202210212623 A CN 202210212623A CN 114935277 B CN114935277 B CN 114935277B
- Authority
- CN
- China
- Prior art keywords
- axis
- coordinate system
- projectile
- guided projectile
- trajectory
- 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.)
- Active
Links
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
Landscapes
- Engineering & Computer Science (AREA)
- General Engineering & Computer Science (AREA)
- Control Of Position, Course, Altitude, Or Attitude Of Moving Bodies (AREA)
- Aerodynamic Tests, Hydrodynamic Tests, Wind Tunnels, And Water Tanks (AREA)
Abstract
The invention discloses an online planning method for ideal trajectory of a gliding Cheng Zhidao projectile. The method comprises the following steps: establishing a ground, translation, projectile, trajectory and speed coordinate system; establishing a motion equation of the 3dof guided projectile; performing wind tunnel test to obtain pneumatic parameters; establishing performance indexes of a ballistic planning problem according to the striking requirements, and establishing constraint conditions; planning an ideal trajectory at the initial ascending section or boosting section after the launching, acquiring and processing satellite data at the moment, and substituting the data into a guided projectile motion equation as a starting point; the grid is reconstructed by a self-adaptive hp method, the number of grids and the order of an interpolation polynomial are optimized, the problem of discrete optimal control by a pseudo-spectrum method is converted into a nonlinear programming problem, and an ideal trajectory curve, a speed curve, a trajectory dip angle curve, an overload curve and a rudder deflection angle curve of the centroid motion of the guided projectile are obtained by solving. The method is insensitive to the interference suffered by the guided projectile before the on-line planning of the ideal trajectory, and has high striking precision.
Description
Technical Field
The invention belongs to the field of guidance trajectory, and particularly relates to an online planning method for an ideal trajectory of a gliding increase Cheng Zhidao projectile.
Background
The problem of trajectory planning of guided projectiles has wide application in the technical field of engineering, and an ideal projectile 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 strength. On-line planning of ideal trajectory is a key to reducing burden of fire control system and difficulty of trajectory tracking. The related motion parameters of the guided projectile at the moment of starting on-line trajectory planning are quite different due to the deviation of parameters such as firing angle, deflection angle, initial speed and the like of the guided projectile during launching and the interference of a low-altitude ascending section. The problem of online planning of the ideal trajectory of gliding increase Cheng Zhidao projectiles is a hot spot of research in recent years.
The ballistic planning problem is essentially an optimal control problem with state constraints and control constraints that meet certain criteria, and as the system and constraint conditions become more complex, it is usually solved by numerical methods, with both indirect and direct methods.
Document hypersonic jump-glide trajectory plan design and optimization (Zhao Jisong et al, solid rocket technology, volume 32, phase 2, pages 123-126, 2009) adopts an indirect method to convert the trajectory optimization problem into a two-point boundary value problem, and adopts a conjugate gradient method to perform solution. Literature 'research on ascending track optimization method based on indirect method' (Wu Jialiang, navigation positioning and time service, 3 rd volume, 2 nd period, 14 th to 19 th pages, 2016 years) is used for solving the problem of optimizing the ascending track of the solid rocket, and the indirect method is adopted to obtain the optimal track under the real atmospheric environment under the condition of meeting the attack angle constraint. The indirect method mainly uses the Pontryagin extremum principle to solve the first-order essential condition of the optimal control problem to obtain the Hamiltonian multipoint boundary problem, and the system state, the control input and the parameters meet the two-point boundary problem through numerical calculation to obtain the optimal solution. However, when the optimal first-order necessary condition is deduced by the Pontryagin extremum principle, a cooperative variable is introduced, the method is sensitive to an initial value, the limitation is large, the two-point edge value is difficult to solve when the complex optimal control problem is processed, and the method is limited in practical engineering application.
The literature ' fast track generation under hypersonic unpowered remote gliding aircraft multi-constraint conditions ' (Yong ' mi and the like, astronautics report, 1 st period, 46-52 pages, 2008) uses a p method to update grids, changes interpolation polynomial dimensions to achieve convergence, and proposes a serial solving strategy capable of generating initial values based on Gaussian pseudo-spectrum method to cope with complex multi-constraint conditions. The direct method comprises a target shooting method and a point matching method, the wireless dimension optimization problem is converted into a finite dimension nonlinear programming problem through a parameterization method, and a nonlinear programming problem solver is used for directly obtaining control variables required by the original optimal control problem, so that an optimal ideal trajectory is solved. The targeting method only disperses the control variable, and cannot solve the complex optimal control problem, so the targeting method is rarely used; the point-matching method is used for dispersing control variables and state variables, wherein the pseudo-spectrum method is widely applied due to good convergence, the pseudo-spectrum method combined with the p method is fixed due to nodes, the problem of discontinuity or non-smoothness is difficult to solve, and dimension disasters are caused when the dimension is excessively increased.
The existing gliding increase Cheng Zhidao shell is used for completing accurate shooting of targets through control of a fire control system, particularly the carrier-based gliding increase Cheng Zhidao shell, ideal trajectory planning is completed through the carrier fire control system before shooting, and then shooting is performed, so that a combat task of remote accurate shooting on shore is realized, and accurate requirements are usually set for Circumferential Errors (CEPs) of landing points from the targets, for example, the embodiment requires that the geostationary index CEP is less than 20m, and the shooting speed is not less than 5 shots/min. The number of the warships on the warships usually reaches several hundred to thousands according to the battle requirements, when the battle task is urgent, all warships on the warships are required to shoot at the allowable maximum shooting speed and finish accurate striking, and under the condition, the planning of 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 the ideal trajectory of a gliding increase Cheng Zhidao projectile.
The technical solution for realizing the purpose of the invention is as follows: an online planning method for ideal trajectory of a gliding increase Cheng Zhidao projectile comprises the following steps:
step (1): establishing the ground Axyz, the translation Oxyz and the projectile 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;
step (3): performing wind tunnel test to obtain pneumatic parameters;
step (4): establishing performance indexes of a ballistic planning problem according to the striking requirements, and establishing constraint conditions according to the capability of a launching facility of a fort, the strength limit of a bullet body and the strength requirement of a striking target;
step (5): planning an ideal trajectory at the initial ascending section or boosting section after the launching, acquiring and processing satellite data at the moment, and substituting the data into a guided projectile motion equation as a starting point;
step (6): the grid is reconstructed by a self-adaptive hp method, the number of grids and the order of an interpolation polynomial are optimized, the discrete optimal control problem of a pseudo-spectrum method is converted into a nonlinear programming problem, and an ideal trajectory curve, a speed curve, a trajectory dip angle curve, an overload curve and a rudder deflection angle curve of the mass center motion of the guided projectile are obtained by solving the problem by a nonlinear programming solver.
Further, "establish ground Axyz, translational Oxyz, elastomer Ox" of step (1) 1 y 1 z 1 Trajectory Ox 2 y 2 z 2 And 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 mass center 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 upwards along the vertical line, and the Az axis is vertical to the Ax and Ay;
o is the instantaneous mass center of the elastomer, and Ox, oy and Oz are respectively parallel to and in the same direction as Ax, ay and Az; ox (Ox) 1 The bullet with the axis coincident with the longitudinal axis is positive, oy 1 The axis being in the longitudinal plane, with Ox 1 Vertical, positive upward, oz 1 The axis being perpendicular to Ox 1 y 1 A plane;
Ox 2 the axis coincides with the instantaneous velocity V of the centre of mass of the guided projectile, oy 2 The axis is positioned in the vertical plane containing V and Ox 2 The axis is vertical, oz 2 The axis being perpendicular to Ox 2 Shaft and Oy 2 A shaft; ox (Ox) 3 The axis coincides with V, oy 3 The axis is positioned in the longitudinal symmetry plane and Ox 3 The axis is vertical, oz 3 The axis being perpendicular to Ox 3 y 3 A plane.
Further, the step (2) establishes a motion equation of the 3dof guided projectile, and the specific method comprises the following steps:
neglecting the roll of the guided projectile, the kinetic equation set of the mass center of the guided projectile under the ballistic coordinate system is as follows:
the kinematic equation of the mass center of the guided projectile in the ground coordinate system is as follows:
the mass center mass change equation and the control relation equation are as follows:
where m is the mass of the guided projectile,for mass change rate, the boosting period is not zero only during engine operation, the starting and ending times of the boosting period being dependent on the guided projectile itself, m c The drug loading of the guided projectile is related to the engine combustion time, and is a specific value, V is the instantaneous speed of the mass center of the guided projectile, < >>Is the rate of change of velocity, g is the acceleration of gravity, P is the average thrust of the engine, which is present only during the boosting period, X is the air resistance of the guided projectile, Y is the lift, Z is the lateral force, and the angle of attack α is the velocity in the projectile coordinate system Ox 1 y 1 z 1 Projection onto and Ox of axis 1 Included angle alpha c (t) is the value of the angle of attack alpha changing with time under the action of control law, and the sideslip angle beta is the velocity and projectile coordinate system Ox 1 y 1 z 1 Included angle beta c (t) is the value of the sideslip angle beta changing along with time under the action of control law, f α 、f β The control rudder is deployed to control the velocity tilt angle gamma when the guided projectile flies to the ballistic vertex based on the relation function of the attack angle, the sideslip angle and the two rudder deflection angles which are obtained by deployment based on the transient balance assumption and are related to the aerodynamic parameters, the center of gravity of the guided projectile, the pressure center and the pressure center of the control rudder v Oy, which is the velocity coordinate system 3 The included angle between the axis and Ox2y2 of the ballistic coordinate system, the ballistic tilt angle θ is the included angle between V and horizontal, and +.>Is the rate of change of the inclination angle of the trajectory, the deviation angle of the trajectory +.>Is the angle between the projection of V on the horizontal plane and the Ox axis, +.>Is the rate of change of the trajectory tilt, (x, y, z) is the instantaneous position of the centre of mass of the guided projectile in the ground coordinate system,/->The position change rate of the mass center of the guided projectile in the x axis, the y axis and the z axis of the ground coordinate system is obtained.
Further, the wind tunnel test is carried out in the step (3) to obtain aerodynamic parameters, and the specific method comprises the following steps:
and placing the gliding increase Cheng Zhidao shell model in a wind tunnel, obtaining data of a plurality of groups of zero-lift resistance coefficient, lift coefficient derivative, pitch rudder and yaw rudder zero-lift resistance coefficient, pitch rudder and yaw rudder lift coefficient derivative, missile wing combined pressure center position and control rudder pressure center position corresponding to different Mach numbers according to the action of gas flow and the guided shell model, and selecting an induced resistance coefficient and a control rudder induced resistance coefficient.
Further, the step (4) specifically comprises:
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 Indicating the time required to fly through the ideal trajectory, x (t f ) Representing the coordinates of the ballistic endpoint on the axis of the ground coordinate system Ax;
setting the starting point parameter value of the guided projectile as a starting point constraint, setting the maximum value and the minimum value of the state variable and the control variable in the flight process as process constraint, and setting the landing angle and the landing speed of the guided projectile when landing to be not smaller than the end point constraint set by a user.
Further, the step (5) specifically comprises:
beidou satellite represents position P point (x ε ,y ε ,z ε ) Conversion to position data in a geodetic coordinate system
Wherein a represents a WGS-84 ellipsoid long half shaft, and the numerical value is 6378137.0m; e is the WGS-84 first eccentricity, square 0.0066943799; n is the curvature half of the ellipsoidal mortise unitary circle; the origin w of the WGS-84 coordinate system is the earth centroid, the wz axis points to the ellipsoidal north pole, the wx axis points to the intersection of the equatorial plane and the Greenner meridian plane, and the wy axis forms a right-hand system coordinate system with xwz in the equatorial plane; the ellipsoid short half axis of the geodetic coordinate system coincides with the earth rotation axis wz, the longitude lambda is the included angle between the ellipsoid meridian plane where the P point is located and the Greenwich meridian plane, and the latitudeThe included angle between the normal line of the ellipsoid passing through the P point and the equatorial plane is the height h, which is the distance from the normal line of the ellipsoid to the ellipsoid along the P point;
converting the coordinates of the guided projectile in the geodetic coordinate system into coordinates (x) n ,y n ,z n )：
In (x) ε0 ,y ε0 ,z ε0 ) The origin of the north-east coordinate system is in the WGS-84 coordinate system, the origin n of the north-east coordinate system is a transmitting device point, xnz is the ground horizontal plane of the place where the origin is located, the nx axis points to the north of the earth, the nz axis points to the east of the earth, and the ny axis points to the air perpendicular to the xnz plane;
the north eastern coordinate system to the shooting coordinate system (x, y, z) are:
in the formula, alpha is the direction under the north-east coordinate system, and the shooting coordinate system is the ground coordinate system.
Further, the step (6) specifically comprises:
the motion state equation constraint of optimal control is obtained according to the dynamics, the kinematic equation set and the mass change equation composition of the step (2)Obtaining a performance index J and constraint conditions according to the step (4), and finally determining a description model of the optimal control problem as follows:
wherein J represents a performance index, phi is a final value performance index, t f Indicating the time required to fly through the ideal trajectory,the vectors respectively representing all state variables and all control variables, the 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 DEG, the integral obtains the rudder deflection angle value, +/->Representation->Derivative with respect to time, f is a mapping by mapping +.>And->Representation->q is->And->Function of->State vector representing ballistic endpoint, +.>The state vector representing the ballistic origin, phi is the boundary event constraint function, comprising the sum of the state vector x (t 0 ) And an endpoint state vector x (t) f ) Is a constraint of (2); c 1 And c 2 For both equality and inequality constraints.
Further, the step (6) 'reconstructing the grid by using the adaptive hp method, optimizing the number of grids and the order of the interpolation polynomial, 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): performing initial grid division; setting a maximum allowable error epsilon and a maximum allowable curvature r max The initial number of meshes, and the degree of the Lagrangian interpolation polynomial;
time interval t 0 ,t f ]Dividing into K grid subintervals [ t ] k-1 ,t k ]K=1, …, K, introducing a new time variable τ e [ -1, 1)]The original is transformed by the following formulaThe quantity t is converted into a new variable tau, so that the single-interval nonlinear optimal control problem can be converted into a multi-interval nonlinear optimal control problem:
step (62): calculating the position of an LGR interpolation point, an LGR weight and a Radau pseudo-spectrum differential matrix in each current grid subinterval, and combining a Lagrange interpolation polynomial to approximate and express an optimal control model, so as to convert the optimal control model into a nonlinear programming problem;
let k-th grid subinterval t k-1 ,t k ]The state variables and control variables in the system are denoted as X k (τ) and U k (τ)，N k +1 discrete pointsAnd end point->Forms a section tau e [ -1,1]N in k +2 Lagrangian interpolation points, where +.>Is provided with->Is N k +2 Lagrangian interpolation basis function, < ->Is N k The +1 Lagrangian interpolation basis function has the expression:
state variable X (k) At N k For +2 discrete pointsThe approximation is:
in the method, in the process of the invention,is a state variable X k At the discrete point->A value of the position;
control variable U (k) At N k For +1 discrete pointsThe approximation is:
in the method, in the process of the invention,is a state variable U k At the discrete point->A value of the position;
for state variable X (k) The derivative of the approximation state equation is derived as:
order theThe Radau pseudo-spectrum differential approximation matrix representing the kth subinterval has the following specific expression:
in the formula g (k) (τ) is P n Function of (τ), P n (τ) is a Legendre orthopolynomial, the differential expression is:
then N k The state equation constraints at +1 interpolation points translate into algebraic equation constraints expressed as:
f in the formula s (k) A state equation function vector for the state variable and the control variable;
the LGR integral approximation expresses the integral term in the bolwalk performance index, which approximates:
in the method, in the process of the invention,is a function of the state vector and the control vector, +.>Is the firstDiscrete points of k subintervals->The corresponding LGR weight is specifically formed by:
finally, the optimal control problem is discretized into a nonlinear programming problem expressed as:
wherein φ is a boundary event constraint function, C 1 And C 2 For the constraint of the equality and inequality,is the state variable of the last point of the kth grid,/->Is the state variable of the start point of the k+1th grid;
step (63) calculating the maximum relative error of the kth grid subintervalIf->Then go to the next step, otherwise jump to step (65);
n provided with the kth grid k +1 samples are:
the residuals of the state constraint equation and the process constraint equation at the sampling points are expressed as:
wherein the method comprises the steps ofAnd->Is the residual error,/->Is the state variable and the control variable at the sampling point +.>The value of the key is taken;
taking the maximum residual error of the current interval as the maximum relative error:
step (64): calculating the relative curvature r of the current grid (k) If r (k) ≤r max Re-computing the order of the interpolation polynomial; if r (k) ＞r max Subdividing the current grid, and calculating the position of a new interpolation point; then, the steps (62) to (64) are carried out again until
Relative curvature r (k) The expression of (2) is:
in the method, in the process of the invention,maximum curvature and average curvature for all LGR points in the kth grid;
curvature of each pointThe expression is:
in the method, in the process of the invention,and->Respectively representing the maximum relative error in the kth grid subinterval +.>A first order derivative function and a second order derivative function of the state function at the corresponding LGR point;
if the order of the polynomial is increased, the number of increases Δn can be expressed as:
wherein ceil (·) represents an upward rounding, and B is a positive integer with a value greater than 0 and adjustable;
refinement of the grid, then the interval is repartitioned by a number n of subintervals k The method is obtained by the following formula:
wherein λ is an integer adjustment factor having a value greater than 0;
step (65): let k=k+1, go through step (64) and step (65) for the next grid until all grid subintervalsAnd if the requirements are met, outputting a result and exiting the program.
A gliding increase Cheng Zhidao projectile comprising a missile-borne computer configured to perform online planning of the ideal trajectory of the gliding increase Cheng Zhidao projectile using the method of any one of claims 1-8.
Compared with the prior art, the invention has the remarkable advantages that:
(1) Compared with other ideal trajectory planning methods meeting the requirement of hitting targets, the on-line ideal trajectory planning method provided by the invention reduces the burden of a fire control system by carrying out on-line trajectory planning on each guided projectile through the missile-borne computer.
(2) Aiming at the problem that the guided projectile is interfered by gusts and the like in an initial ascending section and deviates from an ideal trajectory too far before reaching an actual trajectory vertex, compared with the situation that a target cannot be hit due to inapplicability to interference of a control law corresponding to the ideal trajectory planned by a fire control system, the method provided by the invention is used for carrying out online planning on a boosting section, and the position and motion state information of the guided projectile are changed due to interference during online planning, so that the error between the ideal trajectory planned on the basis and the actual trajectory is smaller, the difficulty in tracking the ideal trajectory by the guided projectile is reduced, and the target can be accurately hit.
(3) Compared with the problems of initial value sensitivity and difficult solution of two-point side values brought by an indirect method, the technical method provided by the invention relies on a pseudo-spectrum method, the discrete points are obtained by using a Newton iteration method in a continuous time period, and the equation approximation is performed after selecting a proper polynomial order, so that two difficulties of the indirect method are avoided, and the pseudo-spectrum method has good convergence and is easy to calculate when solving the problem of infinity.
(4) Compared with the problem that the pseudo-spectrum method based on the p method cannot solve the discontinuous or unsmooth problem, the technical method provided by the invention is based on the self-adaptive hp method to update the grid, so that the errors and curvatures on sampling points in the grid meet the requirements, the characteristic that the h method can capture the unsmooth characteristic of the solution is maintained, the exponential convergence speed during the p method solution is also maintained, and the self-adaptive hp method improves the calculation efficiency and the resolution precision of the nonlinear programming problem.
Drawings
FIG. 1 is a graph of the angular relationship between an projectile coordinate system and a velocity coordinate system established in the present invention.
Fig. 2 is a graph of the angular relationship between the translational coordinate system and the ballistic coordinate system established by the present invention.
FIG. 3 is a graph of the angular relationship between the ballistic coordinate system and the velocity coordinate system established in the present invention.
Fig. 4 is a longitudinal flight plot of an ideal trajectory of a guided projectile for on-line planning.
Fig. 5 is a graph of the lateral flight of an on-line planned ideal trajectory of a guided projectile.
Fig. 6 is a velocity profile of an ideal trajectory of a guided projectile for on-line planning.
Fig. 7 is a graph of the ballistic tilt angle of an ideal trajectory of a guided projectile for on-line planning.
Fig. 8 is an overload plot of an ideal trajectory of a guided projectile for on-line planning.
Fig. 9 is a graph of change in pitch rudder deflection angle for an on-line planned ideal trajectory of a guided projectile.
Fig. 10 is a graph of yaw rudder deflection angle change of an ideal trajectory of a guided projectile planned on-line.
Fig. 11 is an ideal ballistic diagram of a guided projectile based on offline planning and online planning during gust disturbance.
Detailed Description
The invention is described in further detail below with reference to the accompanying drawings.
In order to make the objects, technical solutions and advantages of the present application more apparent, the present application will be further described in detail with reference to the accompanying drawings and examples. It should be understood that the specific embodiments described herein are for purposes of illustration only and are not intended to limit the present application.
Based on an optimal control theory, the method considers the initial emission deviation and the interference of the guided projectile before the trajectory planning on line, can effectively plan the ideal trajectory meeting the requirements of drop points, drop speeds and drop angles on line through a self-adaptive pseudo-spectrum method, and provides a high-efficiency and accurate calculation method for the on-line planning of researching the ideal trajectory of the gliding increase Cheng Zhidao projectile. An online planning method for ideal trajectory of a gliding increase Cheng Zhidao shell specifically comprises the following steps:
setting up the ground Axyz, translational Oxyz and projectile Ox in the step (1) 1 y 1 z 1 Trajectory Ox 2 y 2 z 2 And 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 mass center 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 upwards along the vertical line, and the Az axis is vertical to the Ax and Ay;
o is the instantaneous mass center of the elastomer, and Ox, oy and Oz are respectively parallel to and in the same direction as Ax, ay and Az; ox (Ox) 1 The bullet with the axis coincident with the longitudinal axis is positive, oy 1 The axis being in the longitudinal plane, with Ox 1 Vertical, positive upward, oz 1 The axis being perpendicular to Ox 1 y 1 A plane;
Ox 2 the axis coincides with the instantaneous velocity V of the centre of mass of the guided projectile, oy 2 The axis is positioned in the vertical plane containing V and Ox 2 The axis is vertical, oz 2 The axis being perpendicular to Ox 2 Shaft and Oy 2 A shaft; ox (Ox) 3 The axis coincides with V, oy 3 The axis is positioned in the longitudinal symmetry plane and Ox 3 The axis is vertical, oz 3 The axis being perpendicular to Ox 3 y 3 A plane.
The angular relationships between the translational, ballistic, and velocity coordinate systems are shown in fig. 1-3.
Step 2, establishing a motion equation of the 3dof guided projectile, wherein the specific method comprises the following steps:
neglecting the roll of the guided projectile, the kinetic equation set of the mass center of the guided projectile under the ballistic coordinate system is as follows:
the kinematic equation of the mass center of the guided projectile in the ground coordinate system is as follows:
the mass center mass change equation and the control relation equation are as follows:
where m is the mass of the guided projectile,for mass change rate, only the engine is operatedThe working period, i.e. the boosting period, is non-zero, the starting and ending times of the boosting period being dependent on the guided projectile itself, m c The drug loading of the guided projectile is related to the engine combustion time, and is a specific value, V is the instantaneous speed of the mass center of the guided projectile, < >>Is the rate of change of velocity, g is the acceleration of gravity, P is the average thrust of the engine, which is present only during the boosting period, X is the air resistance of the guided projectile, Y is the lift, Z is the lateral force, and the angle of attack α is the velocity in the projectile coordinate system Ox 1 y 1 z 1 Projection onto and Ox of axis 1 Included angle alpha c (t) is the value of the angle of attack alpha changing with time under the action of control law, and the sideslip angle beta is the velocity and projectile coordinate system Ox 1 y 1 z 1 Included angle beta c (t) is the value of the sideslip angle beta changing along with time under the action of control law, f α 、f β Based on the relation function of aerodynamic parameters obtained by instantaneous balance assumption expansion, the attack angle, sideslip angle and two rudder deflection angles related to the center of gravity, the center of pressure and the center of pressure of the control rudder of the guided projectile are controlled by expanding the control rudder when the guided projectile flies to the ballistic vertex, and the velocity tilt angle gamma v Oy, which is the velocity coordinate system 3 Ox of axis and trajectory coordinate system 2 y 2 Included angle between the horizontal plane and the V, ballistic inclination angle theta is the included angle between the horizontal plane and the V +>Is the rate of change of the inclination angle of the trajectory, the deviation angle of the trajectory +.>Is the angle between the projection of V on the horizontal plane and the Ox axis, +.>Is the rate of change of the trajectory tilt, (x, y, z) is the instantaneous position of the centre of mass of the guided projectile in the ground coordinate system,/->The position change rate of the mass center of the guided projectile in the x axis, the y axis and the z axis of the ground coordinate system is obtained.
And (3) performing wind tunnel test to obtain pneumatic parameters, wherein the specific method comprises the following steps:
and placing the gliding increase Cheng Zhidao shell model in a wind tunnel, obtaining data of a plurality of groups of zero-lift resistance coefficient, lift coefficient derivative, pitch rudder and yaw rudder zero-lift resistance coefficient, pitch rudder and yaw rudder lift coefficient derivative, missile wing combined pressure center position and control rudder pressure center position corresponding to different Mach numbers according to the action of gas flow and the guided shell model, and selecting an induced resistance coefficient and a control rudder induced resistance coefficient.
The step (4) is specifically as follows:
many engineering applications require a specified range, so here, the minimum flight time is taken as a requirement, and 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 Indicating the time required to fly through the ideal trajectory.
The initial firing speed of the guided projectile is set to 790m/s, the ballistic inclination angle is set to 55 degrees, the positions are set to 0km,0km and 0km, the maximum value and the minimum value of state variables and control variables in the flying process are set to be process constraint, the angle of attack angle, sideslip angle and rudder deflection angle in the flying process cannot exceed 6 degrees, the angle of change of the rudder deflection angle per second cannot exceed 2 degrees, the overload absolute value cannot exceed 3 degrees, the end point parameter value of the guided projectile is set to be end point constraint, the falling angle of the guided projectile when falling to the ground is not less than 65 degrees or is set to be a designated falling angle, and the falling speed is not less than 200m/s.
The step (5) comprises the following steps:
beidou satellite represents position P point (x ε ,y ε ,z ε ) Conversion to position data in a geodetic coordinate system/>
Wherein a represents a WGS-84 ellipsoid long half shaft, and the numerical value is 6378137.0m; e is the WGS-84 first eccentricity, square 0.0066943799; n is the curvature radius of the ellipsoidal mortise unitary circle. The origin w of the WGS-84 coordinate system is the earth centroid, the wz axis points to the ellipsoidal north pole, the wx axis points to the intersection of the equatorial plane and the Greenner meridian plane, and the wy axis forms a right-hand system coordinate system with xwz in the equatorial plane; the ellipsoid short half axis of the geodetic coordinate system coincides with the earth rotation axis wz, the longitude lambda is the included angle between the ellipsoid meridian plane where the P point is located and the Greenwich meridian plane, and the latitudeIs the included angle between the normal of the ellipsoid passing through the point P and the equatorial plane, and the height h is the distance from the normal of the ellipsoid to the ellipsoid along the point P.
Converting the coordinates of the guided projectile in the geodetic coordinate system into coordinates (x) n ,y n ,z n )：
In (x) ε0 ,y ε0 ,z ε0 ) The origin of the north-east coordinate system is the position of the origin in the WGS-84 coordinate system, the origin n of the north-east coordinate system is a transmitting device point, xnz is the ground level of the place where the origin is located, the nx axis points to the north of the earth, the nz axis points to the east of the earth, and the ny axis points to the air perpendicular to the xnz plane.
The north eastern coordinate system to the shooting coordinate system (x, y, z) are:
in the formula, alpha is the direction of the light under the north-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 kinematic equation set and the mass change equation set in the step (2)Obtaining a performance index J and constraint conditions according to the step (4), and finally determining a description model of the optimal control problem as follows:
j represents a performance index, phi is a final value performance index, t f Indicating the time required to fly through the ideal trajectory,the vectors respectively representing all state variables and all control variables, the 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 DEG, the integral obtains the rudder deflection angle value, +/->Representation->Derivative with respect to time, f is a mapping by mapping +.>And->Representation->q is->And->Function of->State vector representing ballistic endpoint, +.>The state vector representing the ballistic origin, phi is the boundary event constraint function, comprising the sum of the state vector x (t 0 ) And an endpoint state vector x (t) f ) Is a constraint of (2); c 1 And c 2 For both equality and inequality constraints.
The step (6) 'reconstructing the grids by using the adaptive hp method, optimizing the number of grids and the order of the interpolation polynomial, 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): performing initial grid division; setting the maximum allowable error epsilon to 0.0001 and the maximum allowable curvature r max 1, the number of initialized meshes is 10, and the order of the Lagrangian interpolation polynomial is 4.
Time interval t 0 ,t f ]Dividing into K grid subintervals [ t ] k-1 ,t k ]K=1, …, K, introducing new onesTime variable τ e [ -1,1]The original variable t is converted into a new variable tau by the following formula, so that the single-interval nonlinear optimal control problem can be converted into a multi-interval nonlinear optimal control problem
Step (62): calculating the position of an LGR interpolation point, an LGR weight and a Radau pseudo-spectrum differential matrix in each current grid subinterval, and combining a Lagrange interpolation polynomial to approximate and express an optimal control model, so as to convert the optimal control model into a nonlinear programming problem;
let k-th grid subinterval t k-1 ,t k ]The state variables and control variables in the system are denoted as X k (τ) and U k (τ)，N k +1 discrete pointsAnd end point->Forms a section tau e [ -1,1]N in k +2 Lagrangian interpolation points, where +.>Is provided with->Is N k +2 Lagrangian interpolation basis function, < ->Is N k The +1 Lagrangian interpolation basis function has the expression:
state variable X (k) At N k For +2 discrete pointsThe approximation is:
in the method, in the process of the invention,is a state variable X k At the discrete point->A value of the position;
control variable U (k) At N k For +1 discrete pointsThe approximation is:
in the method, in the process of the invention,is a state variable U k At the discrete point->A value of the position;
for state variable X (k) The derivative of the approximation state equation is derived as:
order theThe Radau pseudo-spectrum differential approximation matrix representing the kth subinterval has the following specific expression:
in the formula g (k) (τ) is P n Function of (τ), P n (τ) is a Legendre orthopolynomial, the differential expression is:
then N k The state equation constraints at +1 interpolation points translate into algebraic equation constraints expressed as:
f in the formula s (k) A state equation function vector for the state variable and the control variable;
the LGR integral approximation expresses the integral term in the bolwalk performance index, which approximates:
in the method, in the process of the invention,is a function of the state vector and the control vector, +.>Discrete point +.>The corresponding LGR weight is specifically formed by:
finally, the optimal control problem is discretized into a nonlinear programming problem expressed as:
/>
wherein φ is a boundary event constraint function, C 1 And C 2 For the constraint of the equality and inequality,is the state vector of the last point of the kth grid,/->A state vector that is the start of the k+1th trellis;
step (63) calculating the maximum relative error of the kth grid subintervalIf->Then go to the next step, otherwise jump to step (65);
n provided with the kth grid k +1 samples are:
the residuals of the state constraint equation and the process constraint equation at the sampling points are expressed as:
wherein the method comprises the steps ofAnd->Is the residual error,/->Is a state variableAnd control variable at sample point +.>The value of the key is taken;
taking the maximum residual error of the current interval as the maximum relative error:
step (64): calculating the relative curvature r of the current grid (k) If r (k) ≤r max Re-computing the order of the interpolation polynomial; if r (k) ＞r max Subdividing the current grid, and calculating the position of a new interpolation point; then, the steps (62) to (64) are carried out again until
Relative curvature r (k) The expression of (2) is:
in the method, in the process of the invention,maximum curvature and average curvature for all LGR points in the kth grid;
curvature of each pointThe expression is:
in the method, in the process of the invention,and->Respectively representing the maximum relative error in the kth grid subinterval +.>A first order derivative function and a second order derivative function of the state function at the corresponding LGR point;
if the order of the polynomial is increased, the number of increases Δn can be expressed as:
wherein ceil (·) represents an upward rounding, and B is a positive integer with a value greater than 0 and adjustable;
refinement of the grid, then the interval is repartitioned by a number n of subintervals k The method is obtained by the following formula:
wherein λ is an integer adjustment factor having a value greater than 0;
step (65): let k=k+1, go through step (64) and step (65) for the next grid until all grid subintervalsAnd if the requirements are met, outputting a result and exiting the program.
The guided projectile of a certain model is made to strike the target (65 km,0 m) at the angle of-65 DEG, the guided projectile is set to be launched at the initial speed 790m/s at the angle of 55 DEG, the time of launching is 0 time, the time of rising to the ballistic vertex is 57.8s, and the missile-borne computer is used for launching at t 0 The parameters of the guided projectile are received at time=20s and the trajectory is planned. When the 20s flight before the guided projectile is undisturbed, the parameter information at the moment should be: position (8585 m,10170m,0 m), velocity 622.94m/s, ballistic tilt angle 44.114 °, ballistic deflection angle 0 °. Assuming that the projectile is launched from t 0 Time of day is unknownThe disturbance, the state of which is biased, was simulated by arbitrarily selecting four sets of ballistic parameters, the parameters are shown in table 1, the unit of distance is (m), the unit of speed is (m/s), and the unit of angle is (°). The simulation results are shown in fig. 4-10. Fig. 4 to 10 show the trajectory and the motion parameter curves of the on-line planned trajectory, respectively, and it can be seen that the on-line planned trajectory with the four sets of parameters of table 1 as the starting points achieves the striking purpose on the premise of meeting the constraint of each variable.
Table 1 randomly selected ballistic programming point parameters
The projectile was fired at v=790 m/s, θ=55°, striking the target at x=40 km at a falling angle of-65 °. And adding gusts with the direction opposite to the spring speed of 0-50 m/s at 18-28 s, and planning a trajectory by using a self-adaptive hp channel pseudo-spectrum method. The first set of experiments uses the control law of offline planning, and the second set plans the trajectory online at t=30s based on the technical method of the invention. FIG. 11 shows two sets of trajectories, where the first set of drop points is 35.41km, the distance error from the target is too large, the strike fails, and the second set meets the requirements of drop angle and set point.
A gliding increase Cheng Zhidao projectile comprising a missile-borne computer configured to perform online planning of the ideal trajectory of the gliding increase Cheng Zhidao projectile using the method of any one of claims 1-8.
Claims (8)
1. An online planning method for ideal trajectory of a gliding increase Cheng Zhidao projectile is characterized by comprising the following steps:
step (1): establishing the ground Axyz, the translation Oxyz and the projectile 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;
step (3): performing wind tunnel test to obtain pneumatic parameters;
step (4): establishing performance indexes of a ballistic planning problem according to the striking requirements, and establishing constraint conditions according to the capability of a launching facility of a fort, the strength limit of a bullet body and the strength requirement of a striking target;
step (5): planning an ideal trajectory at the initial ascending section or boosting section after the launching, acquiring and processing satellite data at the moment, and substituting the data into a guided projectile motion equation as a starting point;
step (6): reconstructing grids by using a self-adaptive hp method, optimizing the number of grids and the order of an interpolation polynomial, converting the discrete optimal control problem into a nonlinear programming problem by using a pseudo-spectrum method, and solving by using a nonlinear programming solver 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 mass center motion of the guided projectile;
step (2) establishing a motion equation of the 3dof guided projectile, wherein the specific method comprises the following steps:
neglecting the roll of the guided projectile, the kinetic equation set of the mass center of the guided projectile under the ballistic coordinate system is as follows:
the kinematic equation of the mass center of the guided projectile in the ground coordinate system is as follows:
the mass center mass change equation and the control relation equation are as follows:
where m is the mass of the guided projectile,for mass change rate, only send outDuring operation of the engine, i.e. the boosting period is not zero, the starting and ending times of the boosting period being dependent on the guided projectile itself, m c The drug loading of the guided projectile is related to the engine combustion time, and is a specific value, V is the instantaneous speed of the mass center of the guided projectile, < >>Is the rate of change of velocity, g is the acceleration of gravity, P is the average thrust of the engine, which is present only during the boosting period, X is the air resistance of the guided projectile, Y is the lift, Z is the lateral force, and the angle of attack α is the velocity in the projectile coordinate system Ox 1 y 1 z 1 Projection onto and Ox of axis 1 Included angle alpha c (t) is the value of the angle of attack alpha changing with time under the action of control law, and the sideslip angle beta is the velocity and projectile coordinate system Ox 1 y 1 z 1 Included angle beta c (t) is the value of the sideslip angle beta changing along with time under the action of control law, f α 、f β The control rudder is deployed to control the velocity tilt angle gamma when the guided projectile flies to the ballistic vertex based on the relation function of the attack angle, the sideslip angle and the two rudder deflection angles which are obtained by deployment based on the transient balance assumption and are related to the aerodynamic parameters, the center of gravity of the guided projectile, the pressure center and the pressure center of the control rudder v Oy, which is the velocity coordinate system 3 Ox of axis and trajectory coordinate system 2 y 2 Included angle between the horizontal plane and the V, ballistic inclination angle theta is the included angle between the horizontal plane and the V +>Is the rate of change of the inclination angle of the trajectory, the deviation angle of the trajectory +.>Is the angle between the projection of V on the horizontal plane and the Ox axis, +.>Is the rate of change of the trajectory tilt, (x, y, z) is the instantaneous position of the centre of mass of the guided projectile in the ground coordinate system,/->The position change rate of the mass center of the guided projectile in the x axis, the y axis and the z axis of the ground coordinate system is obtained.
2. The method of claim 1, wherein the steps of (1) establish the ground Axyz, translational Oxyz, projectile Ox 1 y 1 z 1 Trajectory Ox 2 y 2 z 2 And 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 mass center 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 upwards along the vertical line, and the Az axis is vertical to the Ax and Ay;
o is the instantaneous mass center of the elastomer, and Ox, oy and Oz are respectively parallel to and in the same direction as Ax, ay and Az; ox (Ox) 1 The bullet with the axis coincident with the longitudinal axis is positive, oy 1 The axis being in the longitudinal plane, with Ox 1 Vertical, positive upward, oz 1 The axis being perpendicular to Ox 1 y 1 A plane;
Ox 2 the axis coincides with the instantaneous velocity V of the centre of mass of the guided projectile, oy 2 The axis is positioned in the vertical plane containing V and Ox 2 The axis is vertical, oz 2 The axis being perpendicular to Ox 2 Shaft and Oy 2 A shaft; ox (Ox) 3 The axis coincides with V, oy 3 The axis is positioned in the longitudinal symmetry plane and Ox 3 The axis is vertical, oz 3 The axis being perpendicular to Ox 3 y 3 A plane.
3. The method of claim 2, wherein the wind tunnel test is performed in step (3) to obtain aerodynamic parameters, and the specific method is as follows:
and placing the gliding increase Cheng Zhidao shell model in a wind tunnel, obtaining data of a plurality of groups of zero-lift resistance coefficient, lift coefficient derivative, pitch rudder and yaw rudder zero-lift resistance coefficient, pitch rudder and yaw rudder lift coefficient derivative, missile wing combined pressure center position and control rudder pressure center position corresponding to different Mach numbers according to the action of gas flow and the guided shell model, and selecting an induced resistance coefficient and a control rudder induced resistance coefficient.
4. A method according to claim 3, 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 Indicating the time required to fly through the ideal trajectory, x (t f ) Representing the coordinates of the ballistic endpoint on the axis of the ground coordinate system Ax;
setting the starting point parameter value of the guided projectile as a starting point constraint, setting the maximum value and the minimum value of the state variable and the control variable in the flight process as process constraint, and setting the landing angle and the landing speed of the guided projectile when landing to be not smaller than the end point constraint set by a user.
5. The method of claim 4, wherein step (5) is specifically:
beidou satellite represents position P point (x ε ,y ε ,z ε ) Conversion to position data in a geodetic coordinate system
Wherein a represents a WGS-84 ellipsoid long half shaft, and the numerical value is 6378137.0m; epsilon is WGS-84 first eccentricity, square 0.0066943799; n is the curvature half of the ellipsoidal mortise unitary circle; the origin w of the WGS-84 coordinate system is the earth centroid, the wz axis points to the ellipsoidal north pole, the wx axis points to the intersection of the equatorial plane and the Greenner meridian plane, and the wy axis forms a right-hand system coordinate system with xwz in the equatorial plane; the ellipsoidal minor half axis of the geodetic coordinate system coincides with the earth rotation axis wz, and the longitude lambda is the ellipsoidal meridian plane where the P point is located and the Greenwich mean meridian planeIncluded angle, latitudeThe included angle between the normal line of the ellipsoid passing through the P point and the equatorial plane is the height h, which is the distance from the normal line of the ellipsoid to the ellipsoid along the P point;
converting the coordinates of the guided projectile in the geodetic coordinate system into coordinates (x) n ,y n ,z n )：
In (x) ε0 ,y ε0 ,z ε0 ) The origin of the north-east coordinate system is in the WGS-84 coordinate system, the origin n of the north-east coordinate system is a transmitting device point, xnz is the ground horizontal plane of the place where the origin is located, the nx axis points to the north of the earth, the nz axis points to the east of the earth, and the ny axis points to the air perpendicular to the xnz plane;
the north eastern coordinate system to the shooting coordinate system (x, y, z) are:
in the formula, alpha is the direction under the north-east coordinate system, and the shooting coordinate system is the ground coordinate system.
6. The method according to claim 5, wherein step (6) is specifically:
the motion state equation constraint of optimal control is obtained according to the dynamics, the kinematic equation set and the mass change equation composition of the step (2)Obtaining a performance index J and constraint conditions according to the step (4), and finally determining a description model of the optimal control problem as follows:
wherein J represents a performance index, phi is a final value performance index, t f Indicating the time required to fly through the ideal trajectory,the vectors respectively represent all state variables and all control variables, the control variables are derivatives of two rudder deflection angles with respect to time, the initial value of the rudder deflection angle is 0 DEG, the rudder deflection angle value is obtained by integration, and the values of the rudder deflection angles are calculated by integration>Representation->Derivative with respect to time, f is a mapping by mapping +.>And->Representation->q is->And->Function of->State vector representing ballistic endpoint, +.>The state vector representing the ballistic origin, phi is the boundary event constraint function, comprising the sum of the state vector x (t 0 ) And an endpoint state vector x (t) f ) Is a constraint of (2); c 1 And c 2 For both equality and inequality constraints.
7. The method of claim 6, wherein the step (6) of reconstructing the grid by using the adaptive hp method, optimizing the number of grids and the order of the interpolation polynomial, and converting the discrete optimal control problem to the nonlinear programming problem by using the pseudo-spectrum method comprises the following steps:
step (61): performing initial grid division; setting a maximum allowable error epsilon and a maximum allowable curvature r max The initial number of meshes, and the degree of the Lagrangian interpolation polynomial;
time interval t 0 ,t f ]Dividing into K grid subintervals [ t ] k-1 ,t k ]K=1, …, K, introducing a new time variable τ e [ -1, 1)]The original variable t is converted into a new variable tau by the following formula, so that the single-interval nonlinear optimal control problem can be converted into a multi-interval nonlinear optimal control problem:
step (62): calculating the position of an LGR interpolation point, an LGR weight and a Radau pseudo-spectrum differential matrix in each current grid subinterval, and combining a Lagrange interpolation polynomial to approximate and express an optimal control model, so as to convert the optimal control model into a nonlinear programming problem;
let k-th grid subinterval t k-1 ,t k ]The state variables and control variables in the system are denoted as X k (τ) and U k (τ)，N k +1 discrete pointsAnd end point->Forms a section tau e [ -1,1]N in k +2 Lagrangian interpolation points, whereIs provided with->Is N k +2 Lagrangian interpolation basis function, < ->Is N k The +1 Lagrangian interpolation basis function has the expression:
state variable X (k) At N k For +2 discrete pointsThe approximation is:
in the method, in the process of the invention,is a state variable X k At the discrete point->A value of the position;
control variable U (k) At N k For +1 discrete pointsThe approximation is:
in the method, in the process of the invention,is a state variable U k At the discrete point->A value of the position;
for state variable X (k) The derivative of the approximation state equation is derived as:
order theThe Radau pseudo-spectrum differential approximation matrix representing the kth subinterval has the following specific expression:
in the formula g (k) (τ) is P n Function of (τ), P n (τ) is a Legendre orthopolynomial, the differential expression is:
P 0 (x)＝1
then N k The state equation constraints at +1 interpolation points translate into algebraic equation constraints expressed as:
f in the formula s (k) A state equation function vector for the state variable and the control variable;
the LGR integral approximation expresses the integral term in the bolwalk performance index, which approximates:
in the method, in the process of the invention,is a function of the state vector and the control vector, +.>Discrete point +.>The corresponding LGR weight is specifically formed by:
finally, the optimal control problem is discretized into a nonlinear programming problem expressed as:
wherein φ is a boundary event constraint function, C 1 And C 2 For the constraint of the equality and inequality,is the state variable of the last point of the kth grid,/->Is the state variable of the start point of the k+1th grid;
step (63) calculating the maximum relative error of the kth grid subintervalIf->Then go to the next step, otherwise jump to step (65);
n provided with the kth grid k +1 samples are:
the residuals of the state constraint equation and the process constraint equation at the sampling points are expressed as:
k＝1，2，…，K
k＝1,2,…,K
wherein the method comprises the steps ofAnd->Is the residual error, X s (k) 、U s (k) Is the state variable and the control variable at the sampling point +.>The value of the key is taken;
taking the maximum residual error of the current interval as the maximum relative error:
step (64): calculating the relative curvature r of the current grid (k) If r (k) ≤r max Re-computing the order of the interpolation polynomial; if r (k) >r max Subdividing the current grid, and calculating the position of a new interpolation point; then, the steps (62) to (64) are carried out again until
Relative curvature r (k) The expression of (2) is:
in the method, in the process of the invention,maximum curvature and average curvature for all LGR points in the kth grid;
curvature of each pointThe expression is:
in the method, in the process of the invention,and->Respectively representing the maximum relative error in the kth grid subinterval +.>A first order derivative function and a second order derivative function of the state function at the corresponding LGR point;
if the order of the polynomial is increased, the number of increases Δn can be expressed as:
wherein c εil (·) represents an upward rounding and B is a positive integer with a value greater than 0 and which is adjustable;
refinement of the grid, then the interval is repartitioned by a number n of subintervals k The method is obtained by the following formula:
wherein λ is an integer adjustment factor having a value greater than 0;
step (65): let k=k+1, go through step (64) and step (65) for the next grid until all grid subintervalsAnd if the requirements are met, outputting a result and exiting the program.
8. A gliding increase Cheng Zhidao projectile comprising a missile-borne computer for on-line planning of the ideal trajectory of the gliding increase Cheng Zhidao projectile using the method of any one of claims 1 to 7.
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 CN114935277A (en) | 2022-08-23 |
| CN114935277B true 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) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116362163B (en) * | 2023-06-01 | 2023-09-01 | 西安现代控制技术研究所 | Multi-constraint trajectory rapid optimization method |
Citations (6)
| 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 |
-
2022
- 2022-03-05 CN CN202210212623.XA patent/CN114935277B/en active Active
Patent Citations (6)
| 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 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114935277A (en) | 2022-08-23 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107966156B (en) | Guidance law design method suitable for carrier rocket vertical recovery section | |
| CN108180910B (en) | One kind being based on the uncertain aircraft quick high accuracy method of guidance of aerodynamic parameter | |
| WO2021063073A1 (en) | Method for constructing free trajectory at specified launching elevation angle | |
| CN111591470B (en) | Aircraft precise soft landing closed-loop guidance method adapting to thrust adjustable mode | |
| CN113642122B (en) | Single-shot-table-based remote interception emission data acquisition method and system | |
| CN105180728B (en) | Front data based rapid air alignment method of rotary guided projectiles | |
| CN111444603B (en) | Method for rapidly planning shortest time off-orbit trajectory of recoverable spacecraft | |
| CN110044210B (en) | Closed-circuit guidance on-line compensation method considering arbitrary-order earth non-spherical gravitational perturbation | |
| CN107367941B (en) | Method for observing attack angle of hypersonic aircraft | |
| 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 | |
| CN101694364B (en) | Method for quickly converting perturbation guidance and iteration guidance | |
| CN114935277B (en) | Online planning method for ideal trajectory of gliding Cheng Zhidao shell | |
| Yang et al. | An aerodynamic shape optimization study to maximize the range of a guided missile | |
| CN110220414B (en) | Coincidence method in terminal guided projectile firing plan | |
| CN116719333A (en) | Design method for vertical-launching missile speed vector control turning instruction | |
| CN115574666A (en) | Height setting method for grazing cruise target | |
| CN111649734B (en) | Particle swarm algorithm-based strapdown seeker target positioning method | |
| CN115092421A (en) | Rail forecasting and Lambert orbital transfer method | |
| CN115336429B (en) | Rocket-borne relay measurement and control system phased array antenna beam pointing verification method | |
| Fu et al. | Partial integrated guidance and control method for the interception of nearspace hypersonic target | |
| CN116227340A (en) | Multiple maneuver burst prevention method in zero range line direction | |
| Nobahari et al. | Integrated optimization of guidance and control parameters in a dual spin flying vehicle | |
| Goncharenko et al. | A method to maneuver aeroballistic missiles under counteractions | |
| CN113900448B (en) | Aircraft prediction correction composite guidance method based on sliding mode interference observer |
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 |