Optimal middle guidance trajectory generation method for pneumatically controlled gliding aircraft based on terminal relaxation and voyage convex programming

Abstract

The invention discloses a terminal relaxation and voyage convex programming-based optimal middle guidance trajectory generation method for a pneumatic control gliding aircraft, and relates to the technical field of gliding aircraft guidance. The method comprises the following steps: Step 1, establishing an optimal control problem P0 of a flight range domain of a pneumatic control aircraft; 2, the optimal control problem P0 is converted into a second-order cone problem P2; and Step 3, discretizing the problem P2 from an infinite-dimensional second-order cone problem to a finite-dimensional second-order cone problem, and finally solving the problem P2 through an iterative algorithm, and obtaining an initial reference trajectory generation method. The control quantity can be optimized based on the voyage domain model provided by the invention, and the method does not depend on the monotonicity of the attack angle parabolic surface and the height; the invention further provides a terminal relaxation technology, and the problem that terminal selection is harsh in the optimal problem solving process is solved. The four-order Runge-Kutta method is applied to discretization based on the model, and is more superior to a trapezoidal method and a pseudo-spectral method.

Description

Technical Field

[0001] The present invention relates to the technical field of gliding aircraft guidance, and specifically to an optimal mid-course guidance trajectory generation method for an aerodynamically controlled gliding aircraft based on terminal relaxation and range convex programming. Background Art

[0002] High-speed gliding aircraft pose a major challenge to defense systems due to their high speed, high maneuverability, and large airspace. Mid-course guidance plays an important role in enabling the aircraft to reach the best interception state and successfully intercept high-speed gliding targets in the terminal guidance stage.

[0003] For high-speed gliding targets, mid-course guidance based on dynamic models is widely used. Mid-course guidance trajectory planning can be solved using indirect methods, direct methods, heuristic or intelligent algorithms. The indirect method is based on the maximum principle, transforms the problem into a two-point boundary value problem, and then derives an analytical solution. The advantage of this method is high solution accuracy, but it is difficult to apply to nonlinear models under complex constraints. The direct method discretizes the continuous problem of infinite elements and transforms it into a finite element problem, which is solved using a nonlinear programming algorithm. The advantage of this method is that it does not require the derivation of first-order optimality conditions, but the disadvantage is that the solution time is longer than that of the indirect method. However, with the improvement of computer technology and algorithms, the direct method has made great progress. The pseudospectral method is a type of direct method that uses global polynomials to discretize the problem of infinite elements for solution and can solve complex nonlinear problems with high accuracy. Model predictive static programming also belongs to the direct method. It has relatively simple numerical operations and high computational efficiency, but it is difficult to apply to mid-course guidance with aerodynamic control coupling and severe process constraint problems. Heuristic algorithms are currently a research hotspot, with relatively high solution accuracy, but low computational efficiency for trajectory planning. Intelligent algorithms such as reinforcement learning can solve trajectory planning problems and generate trajectories quickly, but they require pre-training.

[0004] The convex optimization algorithm is also a direct method. Due to its advantages such as high computational efficiency and global optimality, it can provide reliable solutions in a short time. Currently, it is widely applied in directions such as space spacecraft, aerodynamic aerospace (missiles, reentry vehicles), unmanned aerial vehicles, and ground vehicles. The sequential convex programming method can handle highly nonlinear systems and non-convex path constraints by iteratively constructing and solving convex sub-problems until approaching the optimal solution of the original problem. There are certain differences in the solution effects of nonlinear systems in different domains. Existing research shows that solving problems in a distance domain with a decreasing altitude reduces one state dimension of the dynamic equation and has higher solution efficiency. Existing research shows that solving in the time domain relies on the angle-of-attack profile to avoid the coupling problem of dual-channel control (angle of attack and bank angle). Existing research shows that solving in the energy domain reduces the dimension of the velocity state, but it requires predicting the range and terminal velocity, with certain errors. In addition to the domain of the model, different interpolation methods also affect the solution efficiency of the sequential convex programming algorithm. The common ones are mainly divided into three categories: the trapezoidal method, the pseudospectral method, and the adaptive adjustment interpolation method.

[0005] Currently, one of the greatest challenges in sequential convex programming solution lies in the initial guess. Currently, for some models with non-coupled control, this problem has been overcome through means such as variable redefinition, that is, it does not rely on linearization and is not affected by the initial guess. However, for most complex models, it is difficult to avoid this problem. When there is no good initial guess, introducing virtual control can overcome the problem of iterative infeasibility, but this also increases the dimension of one state and the number of iterations. A good initial guess reduces the number of iterations. Existing research uses the indirect method to generate the initial guess trajectory. Existing research uses virtual control to generate the initial trajectory of the first stage as the initial guess trajectory of the second stage to ensure the feasibility of problem solving. Existing research generates the initial guess trajectory using the predictor-corrector method according to the quasi-equilibrium gliding condition. How to balance the generation efficiency and accuracy of the initial guess trajectory is a difficult problem.

[0006] Therefore, we propose an optimal midcourse guidance trajectory generation method for aerodynamically controlled gliding vehicles based on terminal relaxation and range convex programming to solve the problems raised above.

[0007] The above information disclosed in this background technology is only used to increase the understanding of the background technology of the present invention. Therefore, it may include prior art that is not known to those of ordinary skill in the art. Summary of the Invention

[0008] The purpose of the present invention is to provide an optimal midcourse guidance trajectory generation method for aerodynamically controlled gliding vehicles based on terminal relaxation and range convex programming to solve the problems raised in the above background technology.

[0009] To achieve the above object, the present invention provides the following technical solutions: An optimal mid-course guidance trajectory generation method for an aerodynamic control gliding aircraft based on terminal relaxation and range convex programming, comprising the following steps:

[0010] Step1. Establish the optimal control problem P0 for the range domain of the aerodynamic control aircraft:

[0011] 1) Construct a dynamic model for the range domain;

[0012] 2) Consider boundary constraints and process constraints;

[0013] 3) Formulate the optimal control problem P0;

[0014] Step2. Transform the optimal control problem P0 into a second-order cone problem P2:

[0015] 1) Deal with the non-linear dynamics (20) and process constraints (24) through linearization;

[0016] 2) Relax the control constraints (15) to handle the non-convexity of the constraints;

[0017] 3) Ensure relaxation accuracy so that the optimal solution of the relaxed problem satisfies the control constraints (15);

[0018] Step3. Discretize the problem P2 from an infinite-dimensional second-order cone problem to a finite-dimensional second-order cone problem, and finally use an iterative algorithm to solve the problem P2 and obtain an initial reference trajectory generation method:

[0019] 1) Discretize the problem P2 by the fourth-order Runge-Kutta discretization method RK4;

[0020] 2) Initial guess trajectory generation method;

[0021] 3) Obtain the solution of the original problem P0 through the sequential convex programming method:

[0022] Algorithm 1: Solve the original problem P0

[0023] Input: Initial guess trajectory x (0) , trust region δ x , convergence region ε x , k = 0

[0024] Output: States x and u

[0025]

[0026]

[0027] Preferably, in the Step1, specifically constructing the dynamic model for the range domain is as follows:

[0028] First, assume that the aircraft is in the Earth-fixed coordinate system, where the z-axis and x-axis point east and north respectively, forming a right-handed system with the h-axis; the dynamic model of the aerodynamic control aircraft in the time domain has the following non-dimensional dynamic model

[0029]

[0030] where (h, z, x) represents the position coordinates of the aircraft, scaled according to the Earth's radius r e ; r = 1 + h represents the distance of the aircraft from the Earth's center; V represents the relative velocity of the aircraft with respect to the Earth, scaled according to , and g0 represents the acceleration due to gravity at the Earth's radius; θ represents the flight path angle of the target; ψ represents the flight path deviation angle of the target; σ represents the bank angle; L and D respectively represent the non-dimensional lift and drag of the aircraft

[0031]

[0032] where C L and C D respectively represent the lift coefficient and drag coefficient of the aircraft; S represents the area of the aircraft under force; m represents the mass of the aircraft; ρ is the atmospheric density, which can be expressed as

[0033]

[0034] where ρ0 = 1.225 kg / m 3 , and H = 7254.3 m;

[0035] As Figure 1 shown, after determining the initial state and terminal state of the aircraft, the aircraft trajectory generation problem can be transformed into a problem of decreasing the remaining range of the projection of the initial position of the aircraft in the lateral plane onto the direction of the line connecting the initial position and the target position. Let the projected range of a certain path point (h, z, x) be l, and the initial position and target position of the trajectory be (h0, z0, x0) and (h f , z f , x f ), respectively. Then

[0036]

[0037] At this time, the velocity of the aircraft in the direction of the line connecting the initial position and the target position in the horizontal plane is

[0038] v l = dl / dt = v cosθ cos(ψ p - ψ) (6)

[0039] In the formula, ψ p is the line-of-sight angle between the initial point and the terminal of the aircraft in the lateral plane, expressed as

[0040]

[0041] Substituting Equation (6) into Equation (1), the dynamic model of the aircraft in the time domain can be transformed into that in the range domain

[0042]

[0043] The dynamics (8) is nonlinear with respect to the aerodynamic control variables α and σ; an affine system is constructed using the drag polar and approximated as a linear dynamic system; the drag polar is:

[0044] C D (α,M) = C D0 (M) + K(M)C L (α,M) 2 (9)

[0045] where the zero-lift drag coefficient C D0 and the induced drag factor K can be obtained by interpolating aerodynamic data; the lift and drag coefficients corresponding to the maximum lift-drag ratio can be obtained through Equation (9)

[0046]

[0047] Define a normalization coefficient η

[0048]

[0049] Through Equations (9), (10) and (11), we can obtain

[0050]

[0051] At this time, the dimensionless lift and drag accelerations can be expressed as follows

[0052] L = L * η, D = D * [1 + η 2 / 2 (13)

[0053] where

[0054] Next, define new control variables

[0055] u1 = ηcosσ, u2 = ηsinσ, u3 = η 2 (14)

[0056] The new control variables satisfy

[0057]

[0058] Selecting the affine control (15) can limit the upper and lower bounds of the bank angle σ

[0059] σ min ≤σ≤σ max (16)

[0060] Due to the limitations of the angle of attack and bank angle of the aircraft, there are certain constraints on the control variable u; generally, the bank angle control σ of an aerodynamic control aircraft can satisfy σ max =-σ min , Equation (16) can be transformed into

[0061] -u1tanσ max ≤u2≤u1tanσ max (17)

[0062] Assume that η is a non - negative value, with its lower limit being 0 and upper limit being Then the value range of u3 is

[0063]

[0064] wherein

[0065] According to Equations (15), (17), and (18), the control constraints can be constructed as

[0066]

[0067] Next, Equation (8) can be transformed into the non - linear dynamics of an affine system

[0068]

[0069] where x = [h, z, x, v, θ, ψ] T , u = [u1, u2, u3] T ,

[0070]

[0071] Relative to the time domain, in the range domain, the matrix coefficient B(x) before the control variable can satisfy That is, B(x)u has no quadratic term with respect to x.

[0072] Preferably, in the said Step1, considering the boundary constraints and process constraints specifically are:[[]]

[0073] Assume the initial boundary condition x0 of the aircraft, and the initial state is x(l0), where the initial projected range l0 = 0, then the initial state constraint is

[0074] x(l0) = x0 (21)

[0075] Assume the terminal boundary condition x of the aircraftf , the terminal projection range is \(l\) f The state of which is \(x(l\) f ), then the terminal state constraint is

[0076] x(l f ) = \(x\) f (22)

[0077] The overload constraint, heat flux density and dynamic pressure constraint in the process constraint are

[0078]

[0079] For the convenience of expression, the process constraint (23) can be expressed as

[0080]

[0081] Among them, represents the upper bound of the \(j\)-th process constraint.

[0082] Preferably, in the said Step1, the optimal control problem is specifically expressed as:

[0083] Using the nonlinear dynamics equation (20), control constraint (19), boundary constraints (21), (22) and process constraint (24), after adding the objective function, a nonlinear optimal problem with strong equality constraints can be obtained; due to the limited power of the aerodynamic control gliding vehicle, it is difficult to accurately reach the terminal state, and this problem is avoided through the terminal relaxation method; specifically, the terminal state constraint (22) of the equation is transformed as follows

[0084]

[0085] Among them, \(c_1\) and \(c_2\) respectively represent the coefficients before the distance term and the angle term, which are used to adjust the weights of the corresponding parameters; \(\gamma=[\gamma\) h , \(\gamma\) z , \(\gamma\) x , \(\gamma\) θ , \(\gamma\) ψ is the relaxation variable corresponding to the parameter, which does not include the velocity state; the terminal state constraint (25) can be abbreviated as

[0086] |x(l f ) - \(x\) f | \(\leq c\gamma\) (26)

[0087] The optimal midcourse guidance takes the relaxation variable with the penalty term added and the maximum terminal velocity as the objective function, and the form is as follows

[0088] \(J_0 = -\kappa\) v \(v\) f + \(\kappa\) γ \(\sum\) γ \(\gamma\) (27)

[0089] In summary, the mid-course guidance problem of the pneumatic control gliding vehicle can be transformed into the optimal control problem P0, which is in the following form

[0090] P0: min J0

[0091] subject to control constraints (19), dynamic equations (20), initial constraints (21), process constraints (24), and terminal constraints (26);

[0092] In addition, the optimal control problem P0 is currently non-convex, and its non-convexity mainly comes from three sources, namely control constraints (19), non-linear dynamics (20), and process constraints (24).

[0093] Preferably, in the Step2, the specific method of linearizing the non-linear dynamics (20) and process constraints (24) is as follows:

[0094] Assume that the solution in the k-th iteration is {x (k) , u (k)}, where x (k) = [h (k) , z (k) , x (k) , v (k) , θ (k) , ψ (k) T and u (k) = [u1 (k) , u2 (k) , u3 (k) ; It has been described in Step1 Therefore, for the dynamic equation (20), linearizing with respect to {x (k) , u (k)} gives

[0095]

[0096] where, g(x (k) ) = F(x (k) ) - A(x (k) )x (k) ,

[0097] And

[0098]

[0099] In A(x (k) ),

[0100]

[0101] ​

[0102] The second non-convexity is the process constraint (24), which, when linearized with respect to (h (k) , V (k) , u3 (k) ), gives

[0103]

[0104] where

[0105]

[0106] For ease of expression, the process constraint can be described as

[0107] L Uj ≤ 0, j = 1, 2, 3. (30)

[0108] To ensure the effectiveness of the linearization, a trust region constraint

[0109] |x - x (k) | ≤ δ (31)

[0110] where δ ∈ R 5 is a constant vector related to the corresponding state.

[0111] Preferably, in the said Step2, the relaxation of the control constraint (15) is specifically:

[0112] Relax the equality constraint (15) to the inequality constraint (32);

[0113]

[0114] The relaxed control set is as shown in equation (33);

[0115]

[0116] To more clearly describe the relaxation of the control constraint, a schematic diagram is drawn as Figure 2 shown; the control set of the original control constraint (19) is as Figure 2 (a), (c) shown, and the control set is the arc surface in the rotating paraboloid with the angle in [-σ max , σ max , and the maximum radius is It can be seen from this that the set is clearly non-convex; the relaxed control constraint is as Figure 2 (b), (d) shown, and the control set is the space in the rotating paraboloid with the angle in [-σ max , σ max , and the maximum radius is

[0117] Now, after the convexification of problem P0, a second-order cone problem P1 can be obtained:

[0118] P1: min J0

[0119] subject to the initial constraint (21), the terminal constraint (26), the dynamic equation (28), the process constraint (30), the trust-region constraint (31), and the control constraint (33);

[0120] Theorem 1: If {x * (l), u * (l), γ} is the optimal solution of problem P1 and satisfies then {x * (l), u * (l), γ} is also the optimal solution of problem P0;

[0121] Proof: Assume that the optimal objective function values of problems P0 and P1 are and The only difference between problems P0 and P1 lies in the difference between the control constraints (15) and (32); since the control constraint (32) includes the control constraint (15), problem P0 is a feasible subset of problem P1, that is However, the optimal solution of problem P1 always satisfies So this optimal solution is always feasible for problem P0, that is Therefore, The optimal solution of problem P1 is also optimal for problem P0. Similarly, the optimal solution of problem P0 is also optimal for problem P1.

[0122] However, when the optimal solution of the relaxed problem P1 does not satisfy the control constraint (15), this problem is not necessarily the optimal solution of the original problem P0; therefore, it is crucial that the optimal solution of the relaxed problem satisfies the control constraint (15).

[0123] Preferably, in Step 2, ensuring the relaxation accuracy specifically means:

[0124] Directly solving the optimal problem according to the control constraint (33) will inevitably lead to inaccurate relaxation. Therefore, a regularization technique in the objective function is applied to avoid the problem of inaccurate relaxation [Liu2016_Clcd], that is, an integral term of the ballistic deflection angle is added to J0, as shown in the following problem P2

[0125]

[0126] subject to the initial constraint (21), the terminal constraint (26), the dynamic equation (28), the process constraint (30), the trust-region constraint (31), and the control constraint (33);

[0127] Finally, the non-convex problem P0 is transformed into a second-order cone problem P2, and the optimal solution of problem P2 is used to approximately solve the optimal solution of the non-convex original problem P0.

[0128] Preferably, in the above Step3, the discretization of problem P2 by the fourth-order Runge-Kutta discretization method RK4 is specifically as follows:

[0129] The dynamic equation (28) in the second-order cone problem P2 is different from the general dynamic equation in form. The following depends on the reference trajectory {x (k) , u (k)} of the k-th iteration, and the form of RK4 interpolation is as follows

[0130] x i = x i-1 + (k1 + 2k2 + 2k3 + k4)Δl / 6, i = 1, 2,..., N (34)

[0131]

[0132] where Next, substituting equations (35) - (38) into equation (34) and simplifying, we can obtain

[0133]

[0134] where, H i-1 , G i-1 , G i and d have relatively complex specific expression forms, but can be directly obtained; the following defines as the optimization variable of the problem, then equation (39) can be expressed in matrix form

[0135] Ms = D (40)

[0136] where, the matrices M and D can be obtained from equation (39); in addition, the initial constraint (21) can also be incorporated into equation (40); discretizing other constraints in problem P2, problem P2 can be transformed into problem P3 as follows

[0137] P3:

[0138] subject to Ms = D (42)

[0139] |x N - x p | ≤ cγ (43)

[0140] L Ui ≤ 0 (44)

[0141] (u 1,i ) 2 +(u 2,i ) 2 ≤u 3,i (45)

[0142]

[0143] -u 1,i tanσ max ≤u 2,i ≤u 1,i tanσ max (47)

[0144] |x i -x i (k) |≤δ (48)

[0145] Where, i = 0, 1, …, N.

[0146] When there is no process constraint (44), problem P3 becomes problem P4.

[0147] Preferably, in the said Step3, the method for generating the initial guess trajectory is specifically as follows:[[]]

[0148] Select the normalization coefficient and the bank angle σ i as control variables to generate a trajectory group, where σ i is selected proportionally in [σ min , σ max ;

[0149] Select the trajectories l1 and l2 closest to the terminal position in the trajectory group, the distance differences are Δl1 and Δl2, and the bank angles are σ1 and σ2. According to Equation (49), obtain the bank angle σ3;

[0150] σ3 = (σ2Δl1 + σ1Δl2) / (Δl1 + Δl2) (49)

[0151] Select and the bank angle σ3 as control variables to generate the initial reference trajectory.

[0152] Compared with the prior art, the beneficial effects of the present invention are as follows: Based on the voyage domain model proposed by the present invention, the control variables can be optimized respectively, without relying on the monotonicity of the angle of attack parabola and height; the present invention also proposes a terminal relaxation technique to overcome the problem of harsh terminal selection in solving the optimal problem; the model of the present invention is discretized by the fourth-order Runge-Kutta method, which is more superior than the trapezoidal method and the pseudospectral method.

[0153] The above summary is only for the purpose of the specification and is not intended to be limiting in any way. In addition to the illustrative aspects, embodiments, and features described above, further aspects, embodiments, and features of the present invention will be readily apparent by reference to the drawings and the following detailed description. Description of the Drawings

[0154] Figure 1 is the trajectory in the lateral distance domain;

[0155] Figure 2 is a schematic diagram of the original control set and the relaxed control set;

[0156] Figure 3 is the algorithm convergence graph with and without process constraints;

[0157] Figure 4 is the convergence trajectory of different methods;

[0158] Figure 5 is the simulation graph under different terminal conditions;

[0159] Figure 6 is the simulation graph of different interpolation methods and the number of discrete points. Detailed Embodiments

[0160] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0161] I. Numerical Simulation

[0162] A pneumatic control aircraft reentry glide trajectory optimization problem is provided to verify the effectiveness of the present invention. The aircraft parameters are as follows: m = 900 kg, S = 0.4839 m / s 2 , λ ∈ [0, 4.4016], σ ∈ [-60°, 60°]. The path constraint n max = 5g0, q max = 150 kPa. The coefficients κ v = 0.01, κ γ = 100, κ ψ = 0.1.

[0163] For the numerical simulation, the trust region constraint radius and the iterative convergence radius are respectively

[0164]

[0165] Initial Conditions of the Aircraft

[0166] x0 = [60 km, 0 km, 0 km, 2500 m / s, 0°, 0°]

[0167] Terminal Conditions of the Aircraft

[0168] x f = [27 km, 25 km, 400 km, 2500 m / s, 0°, 10°]

[0169] The numerical simulation uses an Intel Core I7 - 10510 2.30 GHZ processor and is programmed with the ECOS solver.

[0170] 1. Verify the Convergence Mode of the Optimal Midcourse Guidance Method Proposed in the Present Invention

[0171] Next, the convergence of the optimal midcourse guidance method proposed in the present invention is verified, and the influence of process constraints (with or without) on the optimal midcourse guidance trajectory is considered. The simulation discretization method uses RK4, and the number of discrete points N = 40.

[0172] The proposed algorithm converges after 4 iterations without process constraints, and the trajectory of each iteration is as shown in Figure 3 (a), (b). It converges after 5 iterations with process constraints, and the trajectory of each iteration is as shown in Figure 3 (c), (d). Figure 3 (e), (f) are respectively the terminal position error and angle error of the optimized trajectory.

[0173] It should be noted that the initial guess trajectory is generated by interpolating a fixed angle of attack and there is a certain gap from the final converged trajectory. Without process constraints, the first iteration still has a large difference from the final converged trajectory because the algorithm fails to find a solution close to the terminal in the trust region based on the initial reference trajectory. The subsequent trajectories converge rapidly, and the third and fourth iterations almost completely overlap; this shows that the proposed algorithm can find a feasible solution based on an inaccurate initial reference trajectory within non - harsh constraints and a fixed trust region, rather than optimization failure. With process constraints, no process constraints are added in the first iteration to ensure that a feasible solution can be generated based on the initial reference trajectory, and the trajectory converges when process constraints are added subsequently, and the fourth and fifth iterations are basically overlapped. In terms of the terminal error of the converged trajectory, the terminal error changes similarly with and without process constraints and is at the same order of magnitude. The convergence performance of the final trajectory is shown in Table 1, and the running time is the average of 20 simulations.

[0174] Table 1 Results of Algorithm Operation

[0175] Process constraints Run time(ms) Error on position(m) Error on angle(°) No 838.3 1.7e-04 5.4e-11 Yes 1267.2 5.6e-05 8.0e-12

[0176] 2. Compare the performance of three interpolation methods, namely RK4, TM (trapezoidal method), FRPM (flipped Radau pseudospectral method), and the GPOPS-II optimization method.

[0177] Compare the solution performance of the sequential convex programming method using the RK4, TM, and FRPM interpolation methods respectively, and the GPOPS method with fixed terminals for solving the optimal midcourse guidance trajectory, that is, compare the solution performance of Problem P2. As can be seen from Figure 4 (a) and (b), the trajectories of different methods are relatively similar, and all reach the expected target position at the expected angle. Figure 4 (c) and (d) show the trajectory inclination angle and trajectory deflection angle of the convergent trajectory. Although there are certain differences in the ballistic inclination angle and ballistic deflection angle of different methods, the trajectory distributions are basically the same. Figure 4 (e) and (f) show that the process constraints and control quantities of the overload, heat flux density, and dynamic pressure of the trajectory meet the constraints, where the black dotted line is the maximum value of the constraint. Generally speaking, the simulation results of RK4 and TM are relatively similar, which is related to the equal spacing of discrete points; FRPM and GPOPS use non-equidistant discrete points, which are somewhat different from the former two.

[0178] The solution performance of different methods for Problem P3 is shown in Table 2, and the running time is the average value of 20 simulations. The terminal errors of the RK4, TM, and FRPM interpolation methods for solving Problem P3 are relatively close and can be ignored. Among the three methods, the single iteration of FRPM takes the longest time, and the single iteration running times of RK4 and TM are close. The running time of the GPOPS method is significantly the longest. Since the terminal relaxation constraint is not adopted, its optimized terminal error is 0, but when the terminal constraint conditions are changed, the GPOPS method may not be able to solve the trajectory.

[0179] Table 2 Comparison of the performance of different methods

[0180] Method Run time(ms) Number of iterations Error on position(m) Error on angle(°) RK4 768.0 5 5.6exp-05 8.0exp-12 TM 574.8 4 9.6exp-05 6.0exp-12 FPRM 2245.8 6 5.8exp-04 1.4exp-11 GPOPS 4814.638 7 0 0

[0181] 3. Conduct simulations with equal spacing of the terminal position to verify the effectiveness of the optimal midcourse guidance method proposed in the present invention.

[0182] Next, taking the RK4 interpolation method as an example, solve Problem P2 under different terminal position constraints to verify the robustness of the method proposed in the present invention. Adopt equidistant selection of terminal position constraints, with the value in the z-axis direction being 20:0.5:30 and the value in the h-axis direction being 24:0.5:30. As Figure 5 (b) shows, the number of discrete points N = 40, and the other simulation parameters remain unchanged.

[0183] Figure 5The simulation shows that under the selected terminal position conditions, the method proposed in the present invention can generate trajectories without infeasible solutions, indicating that the method proposed in the present invention is stable and can generate optimal midcourse guidance trajectories without the need to modify parameters according to specific terminal constraints. Figure 5 In (c) and (d), it is shown that the optimization errors of different terminals are within a certain range, and the optimization errors of position and angle will slightly increase at the boundary values. When the target point is far from the initial reference trajectory, a new initial reference trajectory can be selected as the reference trajectory. In engineering, preset initial reference trajectories can be selected at equal intervals in space or the optimized trajectory can be used as the initial reference trajectory to generate a trajectory reaching the ideal terminal state.

[0184] 4. Analyze the influence of three interpolation methods, RK4, TM, FRPM, and the number of discrete points on solving the optimal midcourse guidance problem

[0185] Analyze the influence of three interpolation methods, RK4, TM, FRPM, and the number of discrete points on solving problem P2, as a reference for future research by oneself and others. The terminal error is the difference between the trajectory terminal and the set target terminal. The integral trajectory terminal refers to the trajectory terminal obtained by interpolating the control quantity of the optimized trajectory at intervals of 100 m on the x-axis and integrating the dynamic equation (28) according to the interpolated control commands.

[0186] Figure 6 (a) and (b) show the position error and angle error of the optimized trajectory. When the number of discrete points is less than 32, the optimization error of TM is slightly worse than that of RK4, and when it is greater than 32, the optimization error is slightly better than that of RK4. The optimization error of FRPM is generally slightly worse than that of TM and RK4 methods, but the optimization error can be ignored for the volume of the aircraft. Figure 6 (c) and (d) show the single-iteration running time and the number of iterations for solving problem P2 using different interpolation methods. Generally speaking, the single-iteration running time increases with the increase in the number of discrete points, and there is no obvious relationship between the number of iterations and the number of discrete points. The running time and the number of iterations of RK4 and TM are close. When the number of discrete points increases, the running time of TM is slightly greater than that of RK4. Compared with RK4 and TM, the running time of FRPM increases significantly with the number of discrete points, the number of iterations is generally larger, and the running time is significantly higher than that of RK4 and TM.

[0187] II. Conclusion

[0188] Aiming at the problems of dual-channel coupling of pneumatic control aircraft and the harsh selection of the trajectory terminal in midcourse guidance, the present invention proposes an optimal midcourse guidance trajectory generation method with terminal relaxation in the distance domain. Based on the dynamic model in the distance domain, it can overcome the dual-channel coupling problem of pneumatic control aircraft without relying on the monotonicity of the angle of attack profile and altitude. Based on terminal relaxation, it overcomes the problems of limited pneumatic maneuverability and the harsh selection of the trajectory terminal. Under reasonable assumptions, it is proved based on the maximum principle that the convexified second-order cone model is equivalent to the non-convex original model.

[0189] Simulations show that the method proposed in the present invention can overcome the problems of dual-channel coupling and the harsh selection of the trajectory terminal in midcourse guidance. Moreover, under the same guessed trajectory and parameter constraints, the trajectory terminal state can be arbitrarily selected within a certain range. In addition, the discrete method of the fourth-order Runge-Kutta method of the model in the present invention has an efficiency similar to that of the trapezoidal method, and the optimization error is smaller when the number of discrete points is small. As the number of discrete points increases, the solution time-consuming of the pseudospectral method is significantly higher than that of the fourth-order Runge-Kutta method discrete and the trapezoidal method.

Claims

1. A method for generating an optimal mid-course guidance trajectory for an aerodynamically controlled glider vehicle based on terminal relaxation and range convex programming, characterized in that: The steps include: Step 1. Establish the optimal control problem P0 of the aerodynamic control aircraft in the range domain: 1) Construct a dynamic model of the range domain; 2) Consider boundary constraints and process constraints; 3) P0 formulation of the optimal control problem; Step 2, transform the optimal control problem P0 into the second-order cone problem P2: 1) Handling nonlinear dynamics (20) and process constraints (24) through linearization; 2) Relax the control constraint (15) to deal with the non-convexity of the constraint; 3) Ensure the accuracy of relaxation so that the optimal solution of the problem after relaxation satisfies the control constraints (15); Step 3: Discretize problem P2 from an infinite-dimensional second-order cone problem to a finite-dimensional second-order cone problem. Finally, solve problem P2 through an iterative algorithm and obtain the initial reference trajectory generation method: 1) Discretize problem P2 using the fourth-order Runge-Kutta discretization method RK4; 2) Initial guess trajectory generation method; 3) Obtain the solution to the original problem P0 through the sequential convex programming method: Algorithm 1: Solve the original problem P0 Input: Initial guess trajectory x (0) , trust region δ x , the convergence region ε x , k = 0 Output: states x and u 2. The method for generating the optimal mid-course guidance trajectory of an aerodynamically controlled gliding vehicle based on terminal relaxation and range convex programming according to claim 1 is characterized in that: In Step 1, the construction of the range domain dynamics model is specifically as follows: First, assume that the aircraft is in the earth-fixed coordinate system, with the z-axis and x-axis pointing to the east and north respectively, forming a right-handed system with the h-axis; the dynamic model of the aerodynamically controlled aircraft is in the time domain, and its dimensionless dynamic model is as follows Among them, (h, z, x) represents the position coordinates of the aircraft, according to the radius of the earth r e Scaling; r = 1 + h represents the distance from the center of the earth of the aircraft; V represents the relative speed of the aircraft to the earth, according to Scaling, g0 represents the gravitational acceleration at the radius of the earth; θ represents the target's track inclination; ψ represents the target's track deviation; σ represents the roll angle; L, D represent the dimensionless lift and drag of the aircraft respectively; Among them, C L and C D They represent the lift and drag coefficients of the aircraft respectively; S represents the force area of ​​the interceptor; m represents the mass of the aircraft; ρ is the atmospheric density, which can be expressed as Where, ρ0 = 1.225 kg / m 3 , H = 7254.3m; After the initial state and terminal state of the aircraft are determined, the aircraft trajectory generation problem can be transformed into the problem of decreasing the remaining distance of the projection of the line connecting the initial position and the target position of the aircraft in the transverse plane. Suppose the projection distance of a path point (h, z, x) is l, and the initial position and target position of the trajectory are (h0, z0, x0) and (h f ,z f ,x f ),but At this time, the speed of the aircraft in the horizontal plane in the direction of the line connecting the initial position and the target position is v l =dl / dt=vcosθcos(ψ p -ψ) (6) In the formula, ψ p is the line of sight angle between the initial point of the aircraft and the terminal in the horizontal plane, expressed as Substituting equation (6) into equation (1), the dynamic model of the aircraft in the time domain can be transformed into the dynamic model in the range domain. The dynamics (8) is nonlinear for the aerodynamic control variables α and σ. The drag polar line is used to construct an affine system and approximate it as a linear dynamic system. The drag polar line is: C D (α,M)=C D0 (M)+K(M)C L (α,M) 2 (9) Where, zero lift drag coefficient C D0 and induced drag factor K can be obtained by interpolating aerodynamic data; the lift and drag coefficients corresponding to the maximum lift-to-drag ratio can be obtained by formula (9): Define a normalization coefficient η Through equations (9), (10) and (11), we can get At this time, the dimensionless lift and drag acceleration can be expressed as follows In the formula, Define new control variables below u1=ηcosσ,u2=ηsinσ,u3=η 2 (14) The new control variable satisfies Selecting affine control (15) can limit the upper and lower limits of the roll angle σ s min ≤σ≤σ max (16) Since the aircraft is limited by the angle of attack and the roll angle, the control variable u has certain constraints; the roll angle control σ of the general aerodynamic control aircraft can satisfy σ max =-σ min , equation (16) can be transformed into -u1tanσ max ≤u2≤u1tanσ max (17) Assume that η is a non-negative value with a lower limit of 0 and an upper limit of Then the value range of u3 is In the formula, According to equations (15), (17), and (18), the control constraints can be constructed as Equation (8) can be transformed into the nonlinear dynamics of an affine system as follows: where \(x = [h,z,x,v,\theta,\psi]\) T , \(u = [u_1,u_2,u_3]\) T , Compared with the time domain, the matrix coefficient B(x) before the control quantity can be realized in the range domain to satisfy That is, B(x)u has no quadratic term with respect to x.

3. The method for generating the optimal mid-course guidance trajectory of an aerodynamically controlled gliding vehicle based on terminal relaxation and range convex programming according to claim 1, characterized in that: In Step 1, the boundary constraints and process constraints are specifically considered as follows: Assuming the initial boundary condition x0 of the aircraft, the initial state is x(l0), where the initial projected range l0 = 0, then the initial state constraint is x(l0)=x0 (21) Assume that the terminal boundary condition of the aircraft is x f , terminal projection range l f The state is x(l f ), then the terminal state constraint is x(l f )=x f (22) The overload constraint, heat flux density and dynamic pressure constraint in the process constraints are: For the convenience of expression, the process constraint (23) can be expressed as in, represents the upper bound of the j-th process constraint.

4. The method for generating the optimal mid-course guidance trajectory of an aerodynamically controlled gliding vehicle based on terminal relaxation and range convex programming according to claim 1, characterized in that: In Step 1, the optimal control problem is specifically expressed as: Using the nonlinear dynamic equation (20), control constraints (19), boundary constraints (21), (22) and process constraints (24), after adding the objective function, a nonlinear optimal problem with strong equality constraints can be obtained. Due to the limited power of the aerodynamically controlled glider, it is difficult to accurately reach the terminal state. This problem can be avoided by the terminal relaxation method. Specifically, the terminal state constraint (22) of the equation is transformed as follows: Where c1 and c2 represent the coefficients before the distance term and angle term, respectively, which are used to adjust the weights of the corresponding parameters; γ = [γ h ,γ z ,γ x ,γ θ ,γ ψ ] is the slack variable of the corresponding parameter, excluding the velocity state; the terminal state constraint (25) can be abbreviated as |x(l f )-x f |≤cγ (26) The optimal mid-course guidance uses the relaxation variable with penalty term and the maximum terminal velocity as the objective function, which is as follows: J0=-k v v f +k γ ∑ γ c (27) In summary, the mid-course guidance problem of aerodynamically controlled gliders can be transformed into an optimal control problem P0, which is as follows: P0:min J0 subject to control constraints (19), dynamic equations (20), initial constraints (21), process constraints (24) and terminal constraints (26); In addition, the optimal control problem P0 is currently non-convex, and its non-convexity mainly comes from three sources, namely control constraints (19), nonlinear dynamics (20) and process constraints (24).

5. The method for generating the optimal mid-course guidance trajectory of an aerodynamically controlled gliding vehicle based on terminal relaxation and range convex programming according to claim 1, characterized in that: In Step 2, the nonlinear dynamics (20) and process constraints (24) are processed by linearization as follows: Assume that the solution in the kth iteration is {x (k) ,u (k) }, where x (k) =[h (k) ,z (k) ,x (k) ,v (k) ,θ (k) ,ψ (k) ] T and u (k) =[u1 (k) ,u2 (k) ,u3 (k) ]; Step 1 has been explained So for the dynamic equation (20), about {x (k) ,u (k) Linearization can be obtained Among them, g(x (k) )=F(x (k) )-A(x (k) )x (k) , And In A(x (k) )middle, The second non-convexity is the process constraint (24), which is related to (h (k) ,V (k) ,u3 (k) ) can be linearized to obtain in, For ease of expression, the process constraints can be described as L Uj ≤0,j=1,2,3. (30) To ensure the validity of linearization, add trust region constraints |x-x (k) |≤δ (31) Among them, δ∈R 5 is a constant vector, related to the corresponding state.

6. The method for generating the optimal mid-course guidance trajectory of an aerodynamically controlled gliding vehicle based on terminal relaxation and range convex programming according to claim 1, characterized in that: In Step 2, the relaxed control constraint (15) is specifically: Relax the equality constraint (15) to the inequality constraint (32); The relaxed control set is shown in formula (33); The control set of the original control constraint (19) is the rotation paraboloid with an angle between [-σ max ,σ max ], the maximum radius is The set is non-convex; the control set of the relaxed control constraints is the rotation parabola with an angle between [-σ max ,σ max ] space, the maximum radius is Now after the problem P0 is convexified, we can get a second-order cone problem P1: P1:min J0 subject to initial constraints (21), terminal constraints (26), dynamic equations (28), process constraints (30), trust region constraints (31) and control constraints (33); Theorem 1: If {x * (l),u * (l),γ} is the optimal solution to problem P1 and satisfies Then {x * (l),u * (l),γ} is also the optimal solution to problem P0; Proof: Assume that the optimal objective function values ​​of problem P0 and problem P1 are and The only difference between problem P0 and problem P1 is the difference between control constraints (15) and (32); since control constraint (32) contains control constraint (15), problem P0 is a feasible subset of problem P1, that is, But the optimal solution to problem P1 always satisfies Therefore, the optimal solution is always feasible for problem P0, that is, therefore, The optimal solution to problem P1 is also optimal for problem P0, and similarly, the optimal solution to problem P0 is also optimal for problem P1.

7. The method for generating the optimal mid-course guidance trajectory of an aerodynamically controlled gliding vehicle based on terminal relaxation and range convex programming according to claim 1, characterized in that: In Step 2, the relaxation accuracy is ensured as follows: Solving the optimal problem directly based on the control constraint (33) is inevitably inaccurate. Therefore, a regularization technique in the objective function is applied to avoid the inaccurate relaxation problem [Liu2016_Clcd], that is, adding an integral term of the trajectory deviation angle in J0, as shown in the following problem P2 P2:min subject to initial constraints (21), terminal constraints (26), dynamic equations (28), process constraints (30), trust region constraints (31) and control constraints (33); Finally, the non-convex problem P0 is transformed into the second-order cone problem P2, and the optimal solution of the P2 problem is used to approximate the optimal solution of the original non-convex problem P0.

8. The method for generating the optimal mid-course guidance trajectory of an aerodynamically controlled gliding vehicle based on terminal relaxation and range convex programming according to claim 1, characterized in that: In Step 3, the problem P2 is discretized by the fourth-order Runge-Kutta discretization method RK4 as follows: The dynamic equation (28) in the second-order cone problem P2 is different from the general dynamic equation The form of is different, and the following depends on the reference trajectory {x (k) ,u (k) }, the form of RK4 interpolation is as follows x i =x i-1 +(k1+2k2+2k3+k4)Δl / 6,i=1,2,...,N (34) in, Substituting equations (35) to (38) into equation (34), we can obtain Among them, H i-1 , G i-1 , G i The specific expression of and d is relatively complicated, but can be directly calculated; the following definition As the optimization variable of the problem, equation (39) can be expressed in the form of a matrix: Ms=D (40) Among them, the matrices M and D can be obtained from equation (39); in addition, the initial constraint (21) can also be incorporated into equation (40); the other constraints in problem P2 are also discretized, and problem P2 can be transformed into problem P3 as shown below subject to Ms=D (42) |x N -x p |≤cγ (43) L Ui ≤0 (44) (in 1,i ) 2 +(in 2,i ) 2 in 3,i (45) -u 1,i tanσ max ≤u 2,i ≤u 1,i tanσ max (47) Where,i=0,1,…,N. Without process constraints (44), problem P3 becomes problem P4.

9. The method for generating the optimal mid-course guidance trajectory of an aerodynamically controlled gliding vehicle based on terminal relaxation and range convex programming according to claim 1, characterized in that: In Step 3, the initial guess trajectory generation method is specifically as follows: Select the normalization coefficient and the tilt angle σ i As the control quantity, a trajectory group is generated, where σ i In [σ min ,σ max ] Medium proportion selection; Select the trajectories l1 and l2 closest to the terminal position in the trajectory group, the distance difference is Δl1 and Δl2, the roll angle is σ1 and σ2, and the roll angle σ3 is obtained according to formula (49); σ3=(σ2Δl1+σ1Δl2) / (Δl1+Δl2) (49) choose And the roll angle σ3 are used as control quantities to generate the initial reference trajectory.

Images