CN113742643A - Low-energy-consumption trajectory optimization method based on system state constraint - Google Patents

Low-energy-consumption trajectory optimization method based on system state constraint Download PDF

Info

Publication number
CN113742643A
CN113742643A CN202110863258.4A CN202110863258A CN113742643A CN 113742643 A CN113742643 A CN 113742643A CN 202110863258 A CN202110863258 A CN 202110863258A CN 113742643 A CN113742643 A CN 113742643A
Authority
CN
China
Prior art keywords
constraint
fitness
trajectory
terminal
value
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
CN202110863258.4A
Other languages
Chinese (zh)
Other versions
CN113742643B (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.)
Beihang University
Original Assignee
Beihang University
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 Beihang University filed Critical Beihang University
Priority to CN202110863258.4A priority Critical patent/CN113742643B/en
Publication of CN113742643A publication Critical patent/CN113742643A/en
Application granted granted Critical
Publication of CN113742643B publication Critical patent/CN113742643B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/13—Differential equations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/12—Computing arrangements based on biological models using genetic models
    • G06N3/126—Evolutionary algorithms, e.g. genetic algorithms or genetic programming

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Biophysics (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Evolutionary Biology (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Operations Research (AREA)
  • Physiology (AREA)
  • Genetics & Genomics (AREA)
  • Artificial Intelligence (AREA)
  • Biomedical Technology (AREA)
  • Computational Linguistics (AREA)
  • Evolutionary Computation (AREA)
  • General Health & Medical Sciences (AREA)
  • Molecular Biology (AREA)
  • Computing Systems (AREA)
  • Feedback Control In General (AREA)

Abstract

The invention discloses a low-energy-consumption trajectory optimization method based on system state constraint, which comprises the following steps of: establishing a projectile kinematics equation set; solving ballistic parameters by taking the angle of attack as a control quantity; setting a constraint condition and a target fitness function; and optimizing the trajectory by using an adaptive genetic algorithm. According to the method, system state constraint conditions are set for a target position area, the optimized trajectory is a low-energy-consumption trajectory meeting process constraint and terminal constraint, instability is avoided in the flight process, the terminal guidance stage is started in a better system state, and a foundation is laid for accurate hitting of terminal guidance and terminal attack angle constraint.

Description

Low-energy-consumption trajectory optimization method based on system state constraint
Technical Field
The invention relates to a low-energy-consumption trajectory optimization method based on system state constraint, and belongs to the technical field of trajectory optimization of guided rocket projectiles.
Background
The guided rocket projectile is additionally provided with the detection guidance device on the basis of the common uncontrolled rocket projectile, has higher hit precision and lower cost, and can effectively fill the gap between the uncontrolled rocket projectile and the guided missile. As shown in figure 1, the flight of the guided rocket projectile is divided into an uncontrolled section, a middle guidance section and a final guidance section, the projectile body flies in the middle guidance section according to a scheme trajectory, and the final guidance flies according to a guidance law to realize the pursuit of a target. Therefore, the final guidance is the key for realizing accurate striking, and the intermediate guidance and the final guidance are in a connection state in order to ensure smooth running.
Aiming at a target position area, setting terminal altitude and speed constraint, selecting a function with low energy consumption as a target fitness function, ensuring that sufficient energy is available in a terminal guidance stage to chase a target, and realizing accurate guidance; the terminal trajectory inclination angle constraint is set, so that the change of the angle in the terminal guidance stage and the realization of the terminal attack angle constraint during striking are facilitated; the set process is restricted to ensure that the projectile body cannot be unstable, and smoothly enters a terminal guidance stage. And optimizing the trajectory through an optimization algorithm, wherein the optimized trajectory meets the terminal constraint, the process constraint and the optimization target, and lays a foundation for smoothly entering terminal guidance and realizing accurate striking.
Disclosure of Invention
In order to ensure that the guided rocket projectile smoothly enters a terminal guidance stage and has a good terminal guidance initial state, the invention provides a low-energy-consumption trajectory optimization method based on system state constraint, which can reduce energy consumption, meet process constraint and terminal constraint, provide a good initial state for terminal guidance, improve guidance precision and ensure the smooth proceeding of terminal attack angle constraint. The method comprises the steps of establishing a kinematic equation set of a projectile body, solving a differential equation set by adopting a Runge-Kutta method, setting constraint conditions and a target fitness function, and optimizing a trajectory by using a self-adaptive genetic algorithm, wherein the method specifically comprises the following steps:
step one, establishing a ballistic kinematics equation set
Preferably, the scheme flight generally selects a fixed vertical plane flight, so that a ballistic kinematics equation set is established on a vertical plane, and the following equation set is obtained through derivation and simplification:
Figure BDA0003186474110000021
in the formula (1), m is projectile mass, V is projectile velocity, alpha is attack angle, theta is trajectory inclination angle, g is gravitational acceleration, and t is flight time; x, Y are the horizontal position and altitude of the projectile respectively; cx0Zero lift resistance, Cx1To induce resistance;
Figure BDA0003186474110000022
is the derivative of lift coefficient, S is characteristic area, Q is dynamic pressure, and is rho V2And ρ is the atmospheric density.
Step two, solving a differential equation set
In the invention, an attack angle alpha is selected as a control variable, discretization is carried out on the attack angle, a trajectory optimization problem is converted into a nonlinear programming problem, and r equal division is carried out on time to obtain U ═ alpha1,α2,...αr,αr+1]As the control variable to be optimized. In order to ensure the continuity of control, a fifth-order polynomial interpolation method is adopted to smooth the control quantity. And solving a differential equation set by adopting a fourth-order Runge-Kutta method according to alpha (t) and an initial state to obtain parameter information such as the position, the speed, the trajectory inclination angle and the like of the trajectory.
Step three, setting constraint conditions and a target fitness function
The selection constraints are as follows:
(1) angle of attack control allowable constraints:
in order to prevent the instability of the projectile body and prevent the projectile body from flying at an overlarge attack angle, the control allowable range is selected as follows:
|α|≤αmax (2)
in the formula (2) < alpha >, (B)maxTo allow control of the maximum value of the absolute value of the angle of attack.
(2) And (4) terminal speed constraint:
in order to ensure that the guided rocket projectile successfully chases the target in the final guidance stage and the projectile body speed cannot be too low, the constraint is set as follows:
V(tf)≥Vmin (3)
v (t) in formula (3)f) Is the velocity value of the terminal moment, VminIs the minimum value allowed by the terminal speed.
(3) And (4) terminal position constraint:
in order to ensure that the guided rocket projectile successfully chases the target in the final guidance stage, the position of the terminal is close to the target area, and in order to chase the target, the altitude of the terminal cannot be too low, the selection constraint is as follows:
Xmin≤X(tf)≤Xmax,Ymin≤Y(tf)≤Ymax (4)
x (t) in the formula (4)f)、Xmin、XmaxRespectively the minimum and maximum values of the terminal horizontal flight position and its constraints, Y (t)f)、Ymin、YmaxRespectively the minimum and maximum of the terminal altitude and its constraints.
(4) And (3) terminal trajectory inclination angle constraint:
in order to ensure that the terminal attack angle constraint of the guided rocket projectile is realized in the final guidance stage, a smaller trajectory inclination angle constraint is set, and the realization of the final guidance terminal constraint angle is facilitated, the constraint is selected as follows:
θmin≤θ(tf)≤θmax (5)
in the formula (5), θ (t)f)、θmin、θmaxRespectively the minimum and maximum values of the terminal ballistic dip and its constraints.
Selecting a target fitness function as follows:
in order to ensure that the energy consumption of ballistic flight is the lowest, the invention selects an optimization target that the reserved energy of a projectile is the largest, namely the consumed energy is the smallest, and selects a target fitness function as follows:
Figure BDA0003186474110000031
equation (6) is the inverse of energy consumption, J is the fitness function,
Figure BDA0003186474110000032
when alpha is a control amount, J is made to be maximum Y (t)0)、V(t0) Respectively, the initial values for altitude and speed.
Step four, optimizing the trajectory by using an adaptive genetic algorithm
The algorithm flow consists of two parts of ballistic parameter and fitness calculation and a genetic algorithm, wherein the ballistic parameter and fitness calculation comprises the following specific steps:
(1) the quintic polynomial interpolation method is adopted to continuously process the initial generation population, so that discrete control variables can be approximated to continuous control functions;
(2) substituting the control function into a system state differential equation system, namely an equation (1), and calculating the motion parameter of the projectile body at each moment by using a Runge-Kutta method;
(3) calculating a target fitness function according to the solved motion parameters, performing penalty function processing on the fitness of the ballistic individual which does not meet the constraint condition set in the third step, wherein the fitness is positive, and reducing the fitness in a certain proportion can effectively eliminate the individual which does not meet the constraint condition, so as to prevent the occurrence of wrong optimal solution, and the penalty function expression is as follows:
Jc＝cJ (7)
formula (III) J, JcRespectively representing the fitness before and after the penalty function is applied, and c is a proportionality coefficient.
The specific execution steps of the genetic algorithm are as follows:
(1) generating a new generation of population: determining the initial population scale N, and determining initial conditions such as parent number M, evolution algebra D and the like;
(2) selecting a parent: selecting M/2 pairs of matrixes from N individuals of the initial generation population by using a selection operator according to the discrete controlled quantity individuals obtained by calculation and the fitness function values corresponding to the discrete controlled quantity individuals;
(3) and (3) cross operation: for the selected M/2 pairs of parents, according to the cross probability PcPerforming cross operation by using a decimal cross operator to generate M intermediate individuals;
(4) and (3) mutation operation: for M intermediate individuals, according to the cross probability PmPerforming mutation operation by using a decimal mutation operator to generate M candidate individuals;
(5) calculating and selecting the fitness of the M candidate individuals to generate a new generation parent, and executing the evolution process again until a termination condition is met;
(6) and calculating and outputting the optimal individual and the corresponding motion parameter and trajectory thereof.
The invention has the advantages and beneficial effects that: the optimized trajectory can consume less flight energy under the condition of meeting process constraint and terminal constraint conditions, sufficient energy can chase the target in the terminal guidance stage, the instability of the projectile body can be prevented by meeting the process constraint, the terminal constraint limits the system state of the projectile body at the terminal, the initial condition of the terminal guidance can enter the terminal guidance stage, and the smooth proceeding of the terminal guidance and the realization of the terminal attack angle constraint are ensured. Therefore, the invention can obtain low-energy-consumption trajectory beneficial to terminal guidance through an optimization algorithm.
Drawings
FIG. 1 is a schematic view of a guided rocket projectile in flight.
Fig. 2 is a schematic flow chart of a ballistic optimization method.
Figure 3 is a diagram of the optimized ballistic trajectory.
Fig. 4 is a graph of angle of attack change.
Figure 5 is a graph of ballistic inclination change.
Fig. 6 is a graph of speed change.
The specific implementation mode is as follows:
step one, establishing a ballistic kinematics equation set
Preferably, the scheme flight generally selects a fixed vertical plane flight, so that a ballistic kinematics equation set is established on a vertical plane, and the following equation set is obtained through derivation and simplification:
Figure BDA0003186474110000041
in the formula (1), m is projectile mass, V is projectile velocity, alpha is attack angle, theta is trajectory inclination angle, g is gravitational acceleration, and t is flight time; x, Y are the horizontal position and altitude of the projectile respectively; cx0Zero lift resistance, Cx1To induce resistance;
Figure BDA0003186474110000051
is the derivative of lift coefficient, S is characteristic area, Q is dynamic pressure, and is rho V2Q＝ρV2And ρ is the atmospheric density.
Step two, solving a differential equation set
In the invention, an attack angle alpha is selected as a control variable, the discretization of the attack angle is convenient for an optimization algorithm, a trajectory optimization problem is converted into a nonlinear programming problem, and r is equally divided for time to obtain U ═ alpha1,α2,...αr,αr+1]As the control variable to be optimized. In order to ensure the continuity in actual control, the control quantity needs to be subjected to smoothing processing before trajectory parameters are solved, the method adopts quintic polynomial interpolation, and is characterized in that the speed and the acceleration of the interpolated quantity are continuous, the control process after interpolation is ensured to be continuous, and the system instability caused by mutation can not occur. The specific calculation method is as follows:
Figure BDA0003186474110000052
in the formula, alpha and alphav、αaThe angle of attack, the rate of change of the angle of attack, and the magnitude of acceleration of the change of the angle of attack are respectively; t is t0Is the initial time, t is the flight time; b0、b1、b2、b3、b4、b5The intermediate function is calculated specifically as follows:
Figure BDA0003186474110000053
t in formula (3)0、t1、α0、α1、
Figure BDA0003186474110000054
The flight time, the attack angle, the change rate of the attack angle and the state of the acceleration of the change of the attack angle at the starting point and the ending point of each section are respectively represented, T is the time interval between two adjacent control points, and d is the difference value of the starting point time and the ending point time of the attack angle.
For U ═ alpha1,α2,...αr,αr+1]And calculating all adjacent two points by using the method to obtain the final smooth continuous control function alpha (t).
And solving a differential equation set by adopting a fourth-order Runge-Kutta method according to alpha (t) and an initial state to obtain parameter information such as the position, the speed, the trajectory inclination angle and the like of the trajectory. The specific solving method is as follows:
one differential equation and its initial values are known as follows:
y'＝f(t,y),y(t0)＝y0 (4)
where y is a function with respect to t, f (t, y) is a function containing t and y, and y (t0)＝y0Is shown at an initial time t0The initial value of (c).
Let tnDenotes the time at which n time intervals h have elapsed after the initial time, ynIs shown at tnThe value of time y.
Y after the next time interval hn+1The values are solved as follows:
Figure BDA0003186474110000061
wherein y isn+1From the value y at the present momentnAdd h times an estimated slope, which is a weighted average of the following four slopes:
Figure BDA0003186474110000062
in the formula
k1Is the slope at the beginning of time；
k2Is the slope of the midpoint of the time segment, and the slope k is adopted by the Euler method1To calculate y at point tnA value of + h/2;
k3is made of k2To calculate the midpoint slope of the y value;
k4is the slope of the end of the time period, using k3Calculating the y value
And selecting a proper step length h in the computer, calculating the y value at any moment after the initial time, and directly selecting the value at a certain moment if the value at the certain moment needs to be used.
Step three, setting constraint conditions and a target fitness function
The process constraint, the terminal constraint, the optimization target and the initial condition are arranged to establish a constraint state function as follows:
Figure BDA0003186474110000063
y (t) in formula (7)0)、V(t0) The altitude and the speed are respectively initial values, the first term is that J tends to the maximum value, namely the energy consumption is minimum, and the trajectory can be ensured to be minimum under the condition of meeting the constraint condition;
in the formula (7) < alpha >maxTo allow control of the maximum value of the absolute value of angle of attack, VminMinimum value allowed for terminal speed, X (t)f)、Xmin、XmaxRespectively the minimum value and the maximum value of the terminal horizontal flight position and the constraint thereof, Y (t)f)、Ymin、YmaxThe minimum and maximum values of terminal altitude and its constraints, θ (t), respectivelyf)、θmin、θmaxRespectively the minimum and maximum values of the terminal ballistic dip and its constraints.
In the formula X0、Y0、θ0Respectively represents the initial value of the control starting time as Xmin、Xmax、、Ymin、Ymax、θmin、θmaxRespectively, the desired constraint values for the system state.
Step four, optimizing the trajectory by using an adaptive genetic algorithm
The process of optimizing the trajectory by the adaptive genetic algorithm is shown in the attached drawing, and specifically comprises the following steps:
the algorithm flow consists of two parts of ballistic parameter and fitness calculation and genetic algorithm, wherein the specific steps of ballistic parameter and fitness calculation are
(1) The quintic polynomial interpolation method is adopted to continuously process the initial generation population, so that discrete control variables can be approximated to continuous control functions; (ii) a
(2) Substituting the control function into a system state differential equation system, namely an equation (1), and calculating the motion parameter of the projectile body at each moment by using a Runge-Kutta method;
(3) calculating a target fitness function according to the solved motion parameters, performing penalty function processing on the fitness of the ballistic individual which does not meet the constraint condition set in the third step, wherein the fitness is positive, and reducing the fitness in a certain proportion can effectively eliminate the individual which does not meet the constraint condition, so as to prevent the occurrence of wrong optimal solution, and the penalty function expression is as follows:
Jc＝cJ (8)
formula (III) J, JcRespectively representing the fitness before and after the penalty function is applied, and c is a proportionality coefficient.
The specific execution steps of the genetic algorithm are as follows:
(1) generating a new generation of population: determining the initial population scale N, and determining initial conditions such as parent number M, evolution algebra D and the like;
(2) selecting a parent: selecting M/2 pairs of parents from N individuals of the initial generation population by using a selection operator according to the calculated control quantity and the fitness thereof; the invention adopts a fitness roulette method, and can select a better individual for multiple times according to the selected probability of the individual. Wherein each individual is selected at a time with a probability of:
Figure BDA0003186474110000071
in the formula (9), N is the total number of individuals, i isIndividual numbering, JiThe fitness function value of the ith individual.
The roulette method is similar to a turntable, the size of the turntable is distributed according to fitness probability, N times of rotation are selected, and the selected individual is left, and the specific selection method comprises the following implementation flows:
1) generating a group of random numbers with the length of N and between 0 and 1 and sequencing the random numbers from small to large;
2) the probabilities of each individual are summed and sorted, which is equivalent to sequential accumulation;
3) and sequentially summing and sequencing the random number and the individual probabilities, comparing the results, if the random number is smaller than the summation probability, selecting the individual corresponding to the summation probability, namely, each individual corresponds to a probability interval, and if the random number is in the interval, selecting the individual.
(3) And (3) cross operation: for the selected M/2 pairs of parents, according to the cross probability PcAnd performing cross operation by using a decimal cross operator to generate M intermediate individuals. In order to improve the convergence speed of the algorithm, the cross probability is subjected to self-adaptive processing, so that invalid cross is reduced, wherein the cross probability is as follows:
Figure BDA0003186474110000081
p in formula (10)c1Is a lower bound of the cross probability, Pc2Is the upper bound of the crossover probability; j. the design is a squaremFor larger fitness values in two matrices crossing each other, JavgAverage of all maternal fitness values, JmaxIs the maximum of all maternal fitness values.
(4) And (3) mutation operation: for M intermediate individuals, according to the cross probability PmAnd performing mutation operation by using a decimal mutation operator to generate M candidate individuals. In order to improve the convergence speed of the algorithm, the mutation probability is subjected to self-adaptive processing to reduce invalid mutation, wherein the mutation probability is as follows:
Figure BDA0003186474110000082
p in formula (11)m1Is a lower bound of the cross probability, Pm2Is the upper bound of the crossover probability; j. the design is a squaremIs the greater fitness of two parents crossing each other, JavgIs the average of all maternal fitness, JmaxIs the maximum value of all maternal fitness.
(5) Carrying out fitness calculation and selection on the M candidate individuals to generate a new generation of parent, and carrying out the evolution process of cross and variation selection again until a termination condition is met;
(6) and calculating and outputting the optimal individual and the corresponding motion parameter and trajectory thereof.
In practical application of the invention, scheme trajectory planning can be performed for different target positions and initial states of the projectile body, and the optimized trajectory meets the selected constraint. The simulation conditions were selected as follows:
(1) initial state: ballistic inclination angle theta0＝0°，X0＝20000m，Y015000m, angle of attack α 00 DEG, speed V0＝350m/s；
(2) Control action tolerance range: the alpha is less than or equal to 10 degrees;
(3) terminal constraint parameters: terminal velocity VfNot less than 450m/s, terminal trajectory inclination angle-15 degree not more than thetafLess than or equal to 10 degrees, and the terminal altitude is less than or equal to 3000m and less than or equal to Yf≤4000m；
(4) The size of the initial generation population N is 1000, the number of the parent M is 100, the evolution generation number D is 50, Pc1＝0.9，Pc2＝0.4，Pm1＝0.0001，Pm2＝0.3。
The simulation results of optimizing the ballistic scheme by using the ballistic optimization method are shown in fig. 3-6.
As can be seen from the simulation results shown in fig. 3 to 6, the optimized trajectory adopts a flight mode of maintaining high-altitude flight, rapidly descending to approach the target, and finally reducing the trajectory inclination angle to make the flight more gentle.
The final terminal horizontal position is Xf52.01km, altitude Yf4.05km, ballistic inclination θfAt-11.5 °, terminal velocity Vf＝558.2m/s，All satisfy the terminal constraint requirements.
Meanwhile, the attack angle does not have frequent and large sudden change, does not exceed the control allowable range of the attack angle, meets the condition that alpha is less than or equal to 10 degrees, and can not cause the instability of the projectile body. The altitude of the projectile body at the early stage is higher, and the air density of the high-altitude projectile body is lower compared with that of the low altitude projectile body, because Q ═ rho V2The dynamic pressure Q is small, from QS (C) in equation (1) of the projectile motionx0+Cx1α2) And
Figure BDA0003186474110000091
it can be seen that a larger angle of attack is required to achieve an equivalent effect. Because the dynamic pressure of the high-altitude flight is low and the energy consumption is low, the early-stage attack angle of the trajectory is large, the change rate of the inclination angle of the trajectory is reduced, the high-altitude flight is maintained, and the range is extended; in the middle stage, the attack angle is reduced, the change rate of the trajectory inclination angle is increased, the energy consumption is less when the aircraft flies at a small attack angle, and the altitude rapidly decreases to approach the target; and finally, increasing attack angle control, increasing the change rate of the trajectory inclination angle, and rapidly reducing the absolute value of the trajectory inclination angle to approach to a target value. And a small trajectory inclination angle enters a final guidance stage, so that a large change space can be generated in the trajectory inclination angle during guidance, and the constraint condition of fixed line-of-sight angle striking in the guidance law can be realized.
Under the condition of meeting the process constraint and the terminal constraint, the target fitness function value of the trajectory is the maximum, namely the energy consumption is the minimum, and the trajectory optimized by the method is a low-energy-consumption trajectory meeting the system state constraint.
The above description is only an embodiment of the present invention, but the scope of the present invention is not limited thereto, and any person skilled in the art can understand that the modifications or substitutions within the technical scope of the present invention are included in the scope of the present invention, and therefore, the scope of the present invention should be subject to the protection scope of the claims.

Claims (8)

1. A low energy consumption trajectory optimization method based on system state constraint is characterized in that: the method comprises the following specific steps:
step one, establishing a ballistic kinematics equation set
Selecting a fixed vertical plane for flying, establishing a ballistic kinematic equation set on the vertical plane, and obtaining the following equation set through derivation and simplification:
Figure FDA0003186474100000011
wherein m is the projectile mass, V is the projectile velocity, alpha is the angle of attack, theta is the ballistic inclination, g is the gravitational acceleration, and t is the flight time; x, Y are the horizontal position and altitude of the projectile respectively; cx0Zero lift resistance, Cx1To induce resistance;
Figure FDA0003186474100000012
is the derivative of lift coefficient, S is characteristic area, Q is dynamic pressure, and is rho V2ρ is the atmospheric density;
step two, solving a differential equation set
Selecting an attack angle alpha as a control variable, discretizing the attack angle, converting the trajectory optimization problem into a nonlinear programming problem, equally dividing time by r to obtain U ═ alpha1,α2,...αr,αr+1]As a control variable to be optimized; in order to ensure the continuity of control, a fifth-order polynomial interpolation method is adopted to carry out smoothing treatment on the control quantity; solving a differential equation set by adopting a fourth-order Runge-Kutta method according to alpha (t) and an initial state to obtain parameter information of the position, the speed and the trajectory inclination angle of the trajectory;
step three, setting constraint conditions and a target fitness function
The selection constraints are as follows:
3.1 Angle of attack control tolerance constraints:
in order to prevent the instability of the projectile body and prevent the projectile body from flying at an overlarge attack angle, the control allowable range is selected as follows:
|α|≤αmax (2)
wherein alpha ismaxTo allow control of the maximum value of the absolute value of the angle of attack;
3.2 terminal speed constraint:
in order to ensure that the guided rocket projectile successfully chases the target in the final guidance stage and the projectile body speed cannot be too low, the constraint is set as follows:
V(tf)≥Vmin (3)
wherein, V (t)f) Is the velocity value of the terminal moment, VminThe minimum value allowed by the terminal speed;
3.3 terminal position constraint:
in order to ensure that the guided rocket projectile successfully chases the target in the final guidance stage, the position of the terminal is close to the target area, and in order to chase the target, the altitude of the terminal cannot be too low, the selection constraint is as follows:
Xmin≤X(tf)≤Xmax,Ymin≤Y(tf)≤Ymax (4)
wherein, X (t)f)、Xmin、XmaxRespectively the minimum and maximum values of the terminal horizontal flight position and its constraints, Y (t)f)、Ymin、YmaxRespectively the minimum value and the maximum value of the terminal altitude and the constraint thereof;
3.4 terminal ballistic dip constraint:
in order to ensure that the terminal attack angle constraint of the guided rocket projectile is realized in the final guidance stage, a smaller trajectory inclination angle constraint is set, and the realization of the final guidance terminal constraint angle is facilitated, the constraint is selected as follows:
θmin≤θ(tf)≤θmax (5)
wherein, θ (t)f)、θmin、θmaxRespectively the terminal trajectory inclination angle and the minimum value and the maximum value of the constraint thereof;
selecting a target fitness function as follows:
in order to ensure that the energy consumption of ballistic flight is the lowest, an optimization target is selected, the reserved energy of a projectile is the largest, namely the consumed energy is the smallest, under the condition that constraint conditions are met, and a target fitness function is selected as follows:
Figure FDA0003186474100000021
wherein J is a fitness function,
Figure FDA0003186474100000022
when alpha is a control amount, J is made to be maximum Y (t)0)、V(t0) Initial values for altitude and speed, respectively;
step four, optimizing the trajectory by using an adaptive genetic algorithm
The method for optimizing the trajectory by using the adaptive genetic algorithm comprises two parts, namely trajectory parameter and fitness calculation and a genetic algorithm, wherein the trajectory parameter and fitness calculation comprises the following specific steps:
4.11 adopting a quintic polynomial interpolation method to continuously process the initial generation population, and approximating the discrete control variable to a continuous control function;
4.12 substituting the control function into a system state differential equation system, namely an equation (1), and calculating the motion parameter of the projectile body at each moment by using a Runge-Kutta method;
4.13 calculate the objective fitness function according to the solved motion parameters, and perform penalty function processing on the fitness of the ballistic individual who does not meet the constraint conditions set in step three, wherein the fitness is positive, and the individual who does not meet the constraint conditions can be effectively eliminated by reducing, so as to prevent the occurrence of wrong optimal solution, and the penalty function expression is as follows:
Jc＝cJ (7)
formula (III) J, JcRespectively representing the fitness before and after the penalty function is applied, wherein c is a proportionality coefficient;
the specific execution steps of the genetic algorithm are as follows:
4.21 Generation of New Generation populations: determining the size N of the initial generation population, and determining the number M of parents and the initial conditions of evolution algebra;
4.22 selection of the parent: selecting M/2 pairs of matrixes from N individuals of the initial generation population by using a selection operator according to the discrete controlled quantity individuals obtained by calculation and the fitness function values corresponding to the discrete controlled quantity individuals;
4.23 intersection operation: for the selected M/2 pairs of parents, according to the cross probability PcApplication ofPerforming cross operation by using a decimal cross operator to generate M intermediate individuals;
4.24 mutation operation: for M intermediate individuals, according to the cross probability PmPerforming mutation operation by using a decimal mutation operator to generate M candidate individuals;
4.25 calculating and selecting the fitness of the M candidate individuals to generate a new generation parent, and executing the evolution process again until a termination condition is met;
4.26 calculating and outputting the optimal individual and the corresponding motion parameter and trajectory.
2. The method of claim 1, wherein the method comprises the following steps: in the second step, the specific calculation method is as follows:
Figure FDA0003186474100000031
in the formula, alpha and alphav、αaThe angle of attack, the rate of change of the angle of attack, and the magnitude of acceleration of the change of the angle of attack are respectively; t is t0Is the initial time, t is the flight time; b0、b1、b2、b3、b4、b5The intermediate function is calculated specifically as follows:
Figure FDA0003186474100000032
t in formula (9)0、t1、α0、α1、
Figure FDA0003186474100000033
Respectively representing the states of the flight time, the attack angle, the change rate of the attack angle and the change acceleration of the attack angle at the starting point and the ending point of each section, wherein T is the time interval between two adjacent control points, and d is the difference value of the starting point time and the ending point time of the attack angle;
for U ═ alpha1,α2,...αr,αr+1]And calculating all adjacent two points by using the method to obtain the final smooth continuous control function alpha (t).
3. The low energy ballistic optimization method based on system state constraints according to claim 1 or 2, characterized in that: in the second step, a fourth-order Runge-Kutta method is adopted to solve a differential equation set through alpha (t) and an initial state, and parameter information of the position, the speed and the trajectory inclination angle of the trajectory is obtained; the specific solving method is as follows:
one differential equation and its initial values are known as follows:
y'＝f(t,y),y(t0)＝y0 (10)
where y is a function with respect to t, f (t, y) is a function containing t and y, and y (t0)＝y0Is shown at an initial time t0An initial value of (d);
let tnDenotes the time at which n time intervals h have elapsed after the initial time, ynIs shown at tnThe value of time y;
y after the next time interval hn+1The values are solved as follows:
Figure FDA0003186474100000041
wherein y isn+1From the value y at the present momentnAdd h times an estimated slope, which is a weighted average of the following four slopes:
Figure FDA0003186474100000042
wherein k is1Is the slope at the beginning of time; k is a radical of2Is the slope of the midpoint of the time segment, and the slope k is adopted by the Euler method1To calculate y at point tnA value of + h/2; k is a radical of3Is made of k2To calculate the midpoint slope of the y value; k is a radical of4Is the end of the time periodSlope of (1), using k3The y value is calculated.
4. The method of claim 3, wherein the method comprises the following steps: in the second step: and selecting a proper step length h in the calculation, calculating the y value at any moment after the initial time, and directly selecting the value at a certain moment if the value at the certain moment needs to be used.
5. The method of claim 1, wherein the method comprises the following steps: in the second step 4.22, adopting a fitness roulette method, and selecting a better individual for multiple times according to the selected probability of the individual; wherein each individual is selected at a time with a probability of:
Figure FDA0003186474100000051
wherein N is the total number of individuals, i is the serial number of the individuals, JiThe fitness function value of the ith individual.
6. The method of claim 5, wherein the method comprises the following steps: selecting N times of rotation for the turntable surface area of the roulette method according to fitness probability distribution, and leaving the selected individuals, wherein the specific selection method comprises the following implementation flows:
1) generating a group of random numbers with the length of N and between 0 and 1 and sequencing the random numbers from small to large;
2) the probabilities of each individual are summed and sorted, which is equivalent to sequential accumulation;
3) and sequentially summing and sequencing the random number and the individual probabilities, comparing the results, if the random number is smaller than the summation probability, selecting the individual corresponding to the summation probability, namely, each individual corresponds to a probability interval, and if the random number is in the interval, selecting the individual.
7. The method of claim 1, wherein the method comprises the following steps: in step two 4.23, the cross probability is processed in a self-adaptive way, so as to reduce invalid cross, and the cross probability is:
Figure FDA0003186474100000052
wherein, Pc1Is a lower bound of the cross probability, Pc2Is the upper bound of the crossover probability; j. the design is a squaremFor larger fitness values in two matrices crossing each other, JavgAverage of all maternal fitness values, JmaxIs the maximum of all maternal fitness values.
8. The method of claim 1, wherein the method comprises the following steps: in step two 4.24, the mutation probability is adaptively processed to reduce the invalid mutation, and the mutation probability is:
Figure FDA0003186474100000053
wherein, Pm1Is a lower bound of the cross probability, Pm2Is the upper bound of the crossover probability; j. the design is a squaremIs the greater fitness of two parents crossing each other, JavgIs the average of all maternal fitness, JmaxIs the maximum value of all maternal fitness.
CN202110863258.4A 2021-07-29 2021-07-29 Low-energy-consumption trajectory optimization method based on system state constraint Active CN113742643B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110863258.4A CN113742643B (en) 2021-07-29 2021-07-29 Low-energy-consumption trajectory optimization method based on system state constraint

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110863258.4A CN113742643B (en) 2021-07-29 2021-07-29 Low-energy-consumption trajectory optimization method based on system state constraint

Publications (2)

Publication Number Publication Date
CN113742643A true CN113742643A (en) 2021-12-03
CN113742643B CN113742643B (en) 2023-06-23

Family

ID=78729399

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110863258.4A Active CN113742643B (en) 2021-07-29 2021-07-29 Low-energy-consumption trajectory optimization method based on system state constraint

Country Status (1)

Country Link
CN (1) CN113742643B (en)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP3211621A1 (en) * 2016-02-29 2017-08-30 The Boeing Company Method and electronic device for optimizing a flight trajectory of an aircraft subject to time of arrival control constraints
CN109506517A (en) * 2018-11-21 2019-03-22 中国人民解放军空军工程大学 A kind of midcourse guidance Method of Trajectory Optimization of belt restraining
CN111442697A (en) * 2020-02-07 2020-07-24 北京航空航天大学 Over-emphasis guidance method and trajectory shaping guidance method based on pseudo-spectrum correction
CN112115544A (en) * 2020-08-10 2020-12-22 南京理工大学 Rocket track optimization method based on improved genetic algorithm

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP3211621A1 (en) * 2016-02-29 2017-08-30 The Boeing Company Method and electronic device for optimizing a flight trajectory of an aircraft subject to time of arrival control constraints
CN109506517A (en) * 2018-11-21 2019-03-22 中国人民解放军空军工程大学 A kind of midcourse guidance Method of Trajectory Optimization of belt restraining
CN111442697A (en) * 2020-02-07 2020-07-24 北京航空航天大学 Over-emphasis guidance method and trajectory shaping guidance method based on pseudo-spectrum correction
CN112115544A (en) * 2020-08-10 2020-12-22 南京理工大学 Rocket track optimization method based on improved genetic algorithm

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
刘运鹏;李伶;: "高超声速导弹高空再入减速段弹道优化设计", 航天控制 *

Also Published As

Publication number Publication date
CN113742643B (en) 2023-06-23

Similar Documents

Publication Publication Date Title
CN110308649B (en) PID parameter optimization method based on PSO-SOA fusion algorithm and applied to industrial process control
CN113325706B (en) Hypersonic aircraft reentry trajectory optimization method based on improved control parameterization
CN111773722B (en) Method for generating maneuver strategy set for avoiding fighter plane in simulation environment
CN108416421A (en) The dynamic Algorithm of Firepower Allocation of bat algorithm is improved based on DDE
CN114063644B (en) Unmanned fighter plane air combat autonomous decision-making method based on pigeon flock reverse countermeasure learning
CN104252132A (en) Adaptive genetic algorithm-based interplanetary orbit control optimization method
CN110046471A (en) Based on the radiator optimization method for improving PSO Neural Network algorithm
CN113741509B (en) Hypersonic gliding aircraft hold-down section energy management method
CN113742643A (en) Low-energy-consumption trajectory optimization method based on system state constraint
Dong et al. Gliding motion optimization for a biomimetic gliding robotic fish
CN114815878B (en) Hypersonic aircraft collaborative guidance method based on real-time optimization and deep learning
CN112464557B (en) Flying wing unmanned aerial vehicle redundant control surface control method based on improved hybrid multi-target PSO
CN107831781B (en) Method and system for controlling movement of robotic fish
CN115685764A (en) Task self-adaptive anti-interference tracking control method and system for variable-span aircraft
CN112379693B (en) Intelligent parallel Gaussian pseudo-spectrum aircraft reentry track optimization method
CN109754108A (en) Unit Economic load distribution method based on fluctuating acceleration coefficient Chaos-Particle Swarm Optimization
CN111077896B (en) Liquid-filled flexible spacecraft parameter optimization method based on improved layering algorithm
CN114580138B (en) Bessel Newton-based missile multi-constraint terminal guidance law design method
CN113311871B (en) Guidance law optimization method and system for jump-glide missile
CN113110576A (en) Self-adaptive fixed time convergence cooperative guidance method capable of realizing continuous switching
CN111897218A (en) Multi-physical-field coupling calculation self-adaptive step length method based on fuzzy control theory
CN116663239B (en) Missile escape area attack distance calculation method based on golden section method
CN114970881B (en) Offline reinforcement learning method and device based on convex hull constraint
CN118331305A (en) Guidance rocket glider maneuvering trajectory self-adaptive planning method considering terminal speed constraint
CN116050259A (en) Trajectory correction method

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