A model predictive static programming terminal guidance method based on control barrier function

Abstract

Aiming at the problem of missile terminal guidance with no-fly zone and impact angle constraints, a model predictive static programming guidance law based on control barrier function is proposed. The proposed algorithm inherits the high efficiency of the traditional model predictive static programming while ensuring the no-fly zone constraints. The main reasons are as follows: first, the projection operator under the linear inequality constraints has an explicit expression, avoiding additional optimization solving; second, due to the forward invariance of the control barrier function, the no-fly zone constraints do not need to be imposed at all discrete nodes, but only at one or a few discrete nodes, reducing the scale of solving. The simulation results show that the proposed algorithm can meet the no-fly zone and impact angle constraints while ensuring high computational efficiency. Subsequent research will focus on trajectory optimization and guidance problems with multiple no-fly zone constraints, even dynamic no-fly zone constraints.

Description

Technical Field

[0001] This invention relates to a terminal guidance method based on control obstacle function and model predictive static programming, which includes no-fly zone and landing angle constraints, and belongs to the field of guidance technology. Background Technology

[0002] With the increasing complexity of the combat environment and the growing sophistication of enemy air defense and anti-missile systems, guidance systems, in addition to meeting basic attack accuracy requirements, often also aim to satisfy attack angle constraints to increase damage effects, attack time constraints to achieve coordinated combat missions, or no-fly zone constraints to improve terminal penetration capabilities. Domestic and international scholars have conducted in-depth research on terminal guidance problems under one or more of these constraints.

[0003] Besides classical proportional guidance and its variations, predictive control, sliding mode variable structure control, Lyapunov theory, and feedback linearization are widely used in guidance law design. Among these, Model Predictive Static Programming (MPSP) has seen significant development. Combining model predictive control and static programming, it can efficiently solve a class of finite-time two-point boundary value problems with terminal equality constraints, and has been successfully applied in many guidance problems such as terminal guidance, ascent guidance, reentry guidance, and powered descent guidance. The solution process of MPSP can be described as follows: first, the guidance problem to be solved is equivalently constructed as a nonlinear least squares problem, then transformed into a series of linear least squares problems, solved successively until convergence. MPSP has high computational efficiency, mainly because: firstly, the problem to be solved is a feasibility problem, which is easier to solve than an optimization problem; secondly, the transformed linear least squares problem, as a simplest type of quadratic convex optimization problem, has an explicit analytical solution and does not depend on any optimization solver.

[0004] However, traditional MPSP methods can only handle terminal equality constraints such as landing angle constraints, but cannot handle process inequality constraints such as no-fly zone constraints. Bhitre et al. used the slack variable method to transform inequality constraints into equality constraints, constructing a high-dimensional dynamic system without process constraints for solution; however, slack variables increase the solution size and reduce convergence performance. Hong et al. constructed a convex programming problem with process constraints within the MPSP framework, using the primordial-dual interior point method to solve iteratively until convergence; however, the constructed convex programming problem no longer has an explicit analytical solution, which increases the solution complexity to some extent.

[0005] Therefore, it is necessary to propose an MPSP method that can handle process constraints and has high computational performance to solve the missile terminal guidance problem considering no-fly zones and angle of impact constraints. Summary of the Invention

[0006] To address the terminal guidance problem that simultaneously considers no-fly zone and angle-of-landing constraints, this invention proposes a MPSP terminal guidance method based on the control obstacle function (CBF). This method, building upon traditional MPSP which can only handle angle-of-landing constraints, utilizes CBF to transform no-fly zone constraints into inequality constraints of control variables. Then, using a projection method, it embeds these control inequality constraints into the traditional MPSP solution framework, resulting in an MPSP terminal guidance method capable of handling both no-fly zone and angle-of-landing constraints.

[0007] The specific technical solution of the present invention is as follows:

[0008] A model prediction static programming final guidance method based on a control barrier function includes the following steps:

[0009] Step 1: Establish the guidance coordinate system and the dynamic equations of a three-degree-of-freedom particle;

[0010] Step 1-1: The origin O of the guidance coordinate system is selected as the projection of the missile's center of mass onto the horizontal plane at the moment of shift change. The Ox axis points towards the target direction, the Oy axis points upward, and the Oz axis is determined according to the right-hand rule.

[0011] Step 1-2: Establish the three-degree-of-freedom particle dynamics equations for the missile as follows:

[0012]

[0013] In the formula: [x,y,z] T Let V be the missile's position vector, V be the velocity magnitude, γ be the trajectory inclination angle, ψ be the trajectory deflection angle, D be the drag magnitude, m be the missile's mass, g be the acceleration due to gravity, and a be the velocity vector. y and a z These are the normal and lateral command accelerations, respectively.

[0014] Steps 1-3: The formula for calculating the resistance D is:

[0015]

[0016] In the formula: Where S is the dynamic pressure, C is the reference area, and S is the dynamic pressure. D The drag coefficient is C. Based on the characteristics of the polar curve, the drag coefficient C... D The calculation formula can be expressed as:

[0017]

[0018] In the formula: C D0 is the zero-lift drag coefficient, and K is the induced drag factor; both are functions of the speed magnitude.

[0019] Steps 1-4: Terminal constraints are expressed as follows:

[0020]

[0021] In the formula: The position vector of the target. Let be the desired landing angle. Assuming the no-fly zone is an infinitely tall cylinder, the no-fly zone constraints are expressed as:

[0022]

[0023] In the formula: x NFZ and z NFZ Located in the center of the no-fly zone, R NFZ This is the radius of the no-fly zone.

[0024] Steps 1-5: Combining formulas (1)-(5), the missile terminal guidance problem can be described as the following optimal control problem P0:

[0025]

[0026] In the formula: x=[x,y,z,V,γ,ψ] T Let u be the state vector, and u = [a y ,a z ] T Let f be the control vector, and f be the compact form of the right-hand side function in the dynamic equations. Given the specified initial state, s is the compact form of the no-fly zone constraint function, h is the compact form of the terminal constraint function, and t... f This refers to the terminal moment of a missile's flight.

[0027] Step 2: Set the independent variable to the actual time t∈[t0,t... f The conversion formula for virtual time τ∈[0,1] is as follows:

[0028] t = t0 + (t f -t0)τ (7)

[0029] The optimal control problem P0 can be rewritten as:

[0030]

[0031] Step 3: The optimal control problem P1 is discretized using the explicit Euler method.

[0032] Select N+1 discrete nodes, uniformly set to 0 = τ0 < τ1 < ... < τ N =1. The state vector x and control vector u are discretized at the above nodes, and the dynamic differential constraints can be transformed into finite-dimensional algebraic constraints, as follows:

[0033] x i+1 =x i +Δτ i g(x i ,u i ), i = 0, 1, ..., N-1 (9)

[0034] In the formula: x i =x(τ) i Let u be the state vector at the i-th node. i =u(τ) i ) represents the control vector at the i-th node. Δτ i Let the distance between the i-th node and the (i+1)-th node be denoted as . Similarly, other constraints are also converted into algebraic constraints at discrete nodes.

[0035] Therefore, the continuous-time optimal control problem P1 is transformed into a finite-dimensional nonlinear programming (NLP) problem, as follows:

[0036]

[0037] Step 4: Next, we will further transform the NLP problem P2 into a constrained nonlinear least squares problem.

[0038] By simultaneously establishing the dynamic constraints and initial conditions in the NLP problem, intermediate state nodes are gradually eliminated, ensuring that the state vector at any given time is solely the control variable {u0, u1, ..., u...}. N-1} and terminal time t f The function. Define the discrete-form state trajectory as... and control history as The terminal constraints are rephrased as follows:

[0039] H(U,t f )=h(x N )=h(u0,u1,…,u N-1 ,t f )=0 (11)

[0040] Therefore, the NLP problem P2 is transformed into the following constrained nonlinear least squares problem:

[0041]

[0042] Where: control history U and terminal time t f As the variable to be solved.

[0043] Step 5: Without considering the no-fly zone constraint, the traditional MPSP transforms the nonlinear least squares problem P3 into a series of linear least squares problems for iterative solution.

[0044] Step 5-1: From the initial conjecture Initially, the iteration format is represented as follows:

[0045]

[0046] In the formula: k = 0, 1, ... is the iteration sequence, and the iteration step size σ k ∈[0,1], and The terminal constraint function H is defined with respect to the variables U and t, respectively. f The Jacobian matrix.

[0047] Step 5-2: Due to the variable to be solved in the formula The dimension of the least squares problem is much larger than the number of constraint functions H, resulting in infinitely many solutions. Among these, the minimum norm solution is the most worthy of study, i.e., finding the solution with the smallest iteration step size, such that the difference between two iteration solutions is minimized. The minimum norm least squares problem is constructed as follows:

[0048]

[0049] The explicit analytical solution to problem P4 can be obtained using the Lagrange multiplier method, as shown below:

[0050]

[0051] Step 5-3: Calculate the control sequence U for the next iteration using the formula. k+1 and terminal time Then, the system dynamic equations are integrated to recover the state trajectory X for the next iteration. k+1 Repeat the above steps until the convergence condition is met.

[0052] Step 6: Using the concept of CBF, transform the no-fly zone constraint into an inequality constraint of the control variable to be solved.

[0053] Step 6-1: The no-fly zone constraints are re-represented as follows:

[0054]

[0055] Based on the dynamic equations, the time derivative of the no-fly zone constraint function can be obtained as follows:

[0056]

[0057] Step 6-2: Since the formula does not explicitly contain control variables, inequality constraints for control variables cannot be established. The following is given... and An approximate calculation formula is provided to avoid cases where the control variables are not explicitly included. According to the dynamic equation, the integral of the position variable x(t) is:

[0058]

[0059] By changing the order of integration, the formula can be rewritten as:

[0060]

[0061] Therefore, the position variable x i The time derivative is approximately calculated as follows:

[0062]

[0063] It can be seen that the approximate calculation formula includes the control variable 'a' to be solved. y and a z .

[0064] Step 6-3: Same as above, the position variable z can be obtained. i The approximate time derivative. Combining formulas (17) and (20), we can obtain t. i No-fly zone constraint function s at time t i The approximate time derivative of is expressed in compact form as follows:

[0065]

[0066] Where: coefficient A j B j Both C and are functions of state variables, expressed as follows:

[0067]

[0068] Using the forward invariance of CBF, the inequality constraints for the control variables are constructed as follows:

[0069]

[0070] Because of the forward invariance of CBF, it is not necessary to impose no-fly zone constraints at all nodes, but only to impose inequality constraints at one or a few nodes.

[0071] Step 7: Using projection technology, the control inequality constraints are embedded into the solution framework of the traditional MPSP, that is, embedded in step 5, to obtain an MPSP algorithm that can handle no-fly zone constraints.

[0072] Step 7-1: In each iteration of MPSP, the state parameters in the formula are replaced by the state variables from the previous iteration. This yields the linear inequality constraints on the control variables to be solved, as follows:

[0073]

[0074] Where: coefficient and Ck It is calculated using the state variables from the previous iteration using formula (22).

[0075] Step 7-2: After obtaining the iterative method and iteration step size using the traditional MPSP, project the unconstrained iterative solution onto the set defined by the linear inequality constraints. The iterative format is as follows:

[0076]

[0077] In the formula: For the projection operator on linear inequality constraints, the iteration direction is... It is also calculated using the formula in step 5.

[0078] Step 7-3: Obtain the control variable U for the next iteration from the formula. k+1 and terminal time Then, the dynamic equations are integrated to recover the state trajectory X for the next iteration. k+1 Repeat the above steps until the convergence condition is met.

[0079] The beneficial effects of this invention are as follows:

[0080] This invention utilizes CBF (Constant-Based Function) technology to establish a direct relationship between no-fly zone constraints and control variables (i.e., inequality constraints regarding control variables). Then, using projection technology, this inequality constraint is embedded into the solution process of MPSP (Multi-Level Programming) without process constraints, achieving guidance law solution with no-fly zone and landing angle constraints. This method inherits the efficiency of traditional model predictive static programming while ensuring the handling of process constraints (no-fly zone constraints). This invention has two advantages: first, it only needs to construct inequality constraints on the CBF time derivative at one or a few discrete nodes, effectively reducing the solution scale, while other methods must consider the no-fly zone constraint at all discrete nodes; second, the constructed CBF time derivative constraint can be regarded as a linear inequality constraint regarding control variables, and its projection operator has an explicit expression, which can effectively improve computational efficiency. Attached Figure Description

[0081] Figure 1 This is a schematic diagram of the three-dimensional motion trajectory in the embodiment;

[0082] Figure 2 This is a schematic diagram of the xOz plane motion trajectory in the embodiment;

[0083] Figure 3 This is the trajectory inclination angle variation curve in the embodiment;

[0084] Figure 4 This is the curve showing the change in normal acceleration in the embodiment;

[0085] Figure 5 This is the lateral acceleration variation curve in the embodiment. Detailed Implementation

[0086] The present invention will be further described below with reference to the accompanying drawings and embodiments.

[0087] This embodiment presents an MPSP terminal guidance method capable of simultaneously handling no-fly zones and angle-of-land constraints, comprising the following steps:

[0088] Step 1: Establish the guidance coordinate system and the dynamic equations of a three-degree-of-freedom particle;

[0089] Step 1-1: The origin O of the guidance coordinate system is selected as the projection of the missile's center of mass onto the horizontal plane at the moment of shift change. The Ox axis points towards the target direction, the Oy axis points upward, and the Oz axis is determined according to the right-hand rule.

[0090] Step 1-2: Establish the three-degree-of-freedom particle dynamics equations for the missile as follows:

[0091]

[0092]

[0093]

[0094]

[0095]

[0096]

[0097] In the formula: [x,y,z] T Let V be the missile's position vector, V be the velocity magnitude, γ be the trajectory inclination angle, ψ be the trajectory deflection angle, D be the drag magnitude, m be the missile's mass, g be the acceleration due to gravity, and a be the velocity vector. y and a z These are the normal and lateral command accelerations, respectively.

[0098] Steps 1-3: The formula for calculating the resistance D is:

[0099]

[0100] In the formula: Where S is the dynamic pressure, C is the reference area, and S is the dynamic pressure. D The drag coefficient is C. Based on the characteristics of the polar curve, the drag coefficient C... D The calculation formula can be expressed as:

[0101]

[0102] In the formula: C D0is the zero-lift drag coefficient, and K is the induced drag factor; both are functions of the speed magnitude.

[0103] Steps 1-4: The terminal constraints can be expressed as:

[0104]

[0105] In the formula: The position vector of the target. Let be the desired landing angle. Assuming the no-fly zone is an infinitely tall cylinder, the no-fly zone constraints can be expressed as:

[0106]

[0107] In the formula: x NFZ and z NFZ Located in the center of the no-fly zone, R NFZ This is the radius of the no-fly zone.

[0108] Steps 1-5: The missile terminal guidance problem can be described as the following optimal control problem:

[0109] P0:min J=||h(x(t f ))|| 2

[0110]

[0111]

[0112] s(x(t))≥0

[0113] In the formula: x=[x,y,z,V,γ,ψ] T Let u be the state vector, and u = [a y ,a z ] T Let f be the control vector, and f be the compact form of the right-hand side function in the dynamic equations. Given the specified initial state, s is the compact form of the no-fly zone constraint function, h is the compact form of the terminal constraint function, and t... f This refers to the terminal moment of a missile's flight.

[0114] Step 2: Set the independent variable to the actual time t∈[t0,t... f The conversion formula for virtual time τ∈[0,1] is as follows:

[0115] t = t0 + (t f -t0)τ

[0116] The optimal control problem P0 can be reconstructed as:

[0117] P1:min J=||h(x(1))|| 2

[0118]

[0119]

[0120] s(x(τ))≥0

[0121] Step 3: Discretize the optimal control problem P1 using the explicit Euler method.

[0122] Select N+1 discrete nodes, uniformly set to 0 = τ0 < τ1 < ... < τ N =1. The state vector x and control vector u are discretized at the above nodes, and the dynamic differential constraints can be transformed into finite-dimensional algebraic constraints, as follows:

[0123] x i+1 =x i +Δτ i g(x i ,u i ), i = 0, 1, ..., N-1

[0124] In the formula: x i =x(τ) i Let u be the state vector at the i-th node. i =u(τ) i Let be the control vector at the i-th node. Similarly, other constraints can be converted into algebraic constraints at discrete nodes.

[0125] Therefore, the continuous-time optimal control problem P1 can be transformed into a finite-dimensional nonlinear programming (NLP) problem, as follows:

[0126] P2:min J=||h(x N )|| 2

[0127] stx i+1 =x i +Δτ i g(x i ,u i ), i = 0, 1, ..., N-1

[0128]

[0129] s i (x i )≥0, i=0,1,…,N

[0130] Step 4: Further transform the NLP problem P2 into a constrained nonlinear least squares problem.

[0131] By combining the dynamic constraints and initial conditions in an NLP problem, intermediate state nodes can be gradually eliminated, so that the state vector at any given time is simply the control variable {u0, u1, ..., u...}. N-1} and terminal time t f The function. Define the discrete-form state trajectory as... and control history as The terminal constraint can be reformulated as:

[0132] H(U,t f )=h(x N )=h(u0,u1,…,u N-1 ,t f ) = 0

[0133] Therefore, the NLP problem P2 can be transformed into the following constrained nonlinear least squares problem:

[0134] P3:min J=||H(U,t f )|| 2

[0135] sts i (U,t f )≥0, i=0,1,…,N

[0136] Where: control history U and terminal time t f As the variable to be solved.

[0137] Step 5: If the no-fly zone constraint is not considered, the traditional MPSP transforms the nonlinear least squares problem P3 into a series of linear least squares problems to be solved iteratively.

[0138] Step 5-1: From the initial conjecture Initially, the iteration format is represented as follows:

[0139]

[0140] In the formula: k = 0, 1, ... is the iteration sequence, and the iteration step size σ k ∈[0,1], and The terminal constraint function H is defined with respect to the variables U and t, respectively. f The Jacobian matrix.

[0141] Step 5-2: Due to the variables to be solved in The dimension of the least squares problem is much larger than the number of constraint functions H, resulting in infinitely many solutions. Among these, the minimum norm solution is the most worthy of study, i.e., finding the solution with the smallest iteration step size, such that the difference between two iteration solutions is minimized. The minimum norm least squares problem is constructed as follows:

[0142]

[0143]

[0144] The explicit analytical solution to problem P4 can be obtained using the Lagrange multiplier method, as shown below:

[0145]

[0146] Step 5-3: Calculate the control sequence U for the next iteration. k+1 and terminal time Then, the system dynamic equations are integrated to recover the state trajectory X for the next iteration. k+1 Repeat the above steps until the convergence condition is met.

[0147] Step 6: Using the concept of CBF, transform the no-fly zone constraint into an inequality constraint of the control variable to be solved.

[0148] Step 6-1: Reconstruct the no-fly zone constraints as follows:

[0149]

[0150] Based on the dynamic equations mentioned above, the time derivative of the no-fly zone constraint function can be obtained as follows:

[0151]

[0152] Step 6-2: In Step 6-1, the time derivative of the no-fly zone constraint does not explicitly contain control variables, making it impossible to establish inequality constraints for the control variables. The following is given... and An approximate calculation formula is provided to avoid cases where the control variables are not explicitly included. Based on the aforementioned dynamic equations, the position variable x(t) can be integrated as:

[0153]

[0154] By changing the order of integration, the position variable x(t) can be reconstructed as follows:

[0155]

[0156] Therefore, the position variable x i The time derivative can be approximated as follows:

[0157]

[0158] It can be seen that the position variable x i The approximate time derivative contains the control variable a to be solved. y and az .

[0159] Step 6-3: Same as above, the position variable z can be obtained. i The approximate time derivative of t. Therefore, t can be obtained. i No-fly zone constraint function s at time t i The approximate time derivative of is expressed in compact form as follows:

[0160]

[0161] Where: coefficient A j B j Both C and are functions of state variables, expressed as follows:

[0162] A j =[-2(x i -x NFZ sinγ j cosψ j +2(z i -z NFZ sinγ j sinψ j ](t i -t j )

[0163] B j =[2(x i -x NFZ )sinψ j +2(z i -z NFZ cosψ j ](t i -t j )

[0164] C≈2(x i -x NFZ )·V i cosγ i cosψ i -2(z i -z NFZ )·V i cosγ i sinψ i

[0165] Using the forward invariance of CBF, the following inequality constraints for the control variables can be constructed:

[0166]

[0167] Because of the forward invariance of CBF, it is not necessary to impose no-fly zone constraints at all nodes, but only to impose inequality constraints at one or a few nodes.

[0168] Step 7: Using projection technology, the constructed control inequality constraints are embedded into the solution framework of the traditional MPSP, that is, embedded in step 5, to obtain an MPSP algorithm that can handle no-fly zone constraints.

[0169] Step 7-1: In each iteration of MPSP, the state parameters in the control inequality are replaced by the state variables from the previous iteration. This yields the linear inequality constraints on the control variables to be solved, as follows:

[0170]

[0171] Step 7-2: After obtaining the iterative method and iteration step size using the traditional MPSP, project the unconstrained iterative solution onto the set defined by the linear inequality constraints. The iterative format is as follows:

[0172]

[0173] In the formula: For the projection operator on linear inequality constraints, the iteration direction is... It is also obtained from the calculation formula in step 5.

[0174] Step 7-3: Obtain the control variable U for the next iteration k+1 and terminal time Then, the dynamic equations are integrated to recover the state trajectory X for the next iteration. k+1 Repeat the above steps until the convergence condition is met. Specific implementation examples:

[0176] This embodiment performs a missile terminal guidance simulation with no-fly zone constraints and landing angle constraints. All simulations are executed on a desktop computer equipped with an Intel i5-12500H processor with a main frequency of 2.50GHz. All programs are compiled and run in the MATLAB environment.

[0177] The missile has a mass m = 150 kg and a reference area S = 0.0324 m². 2 Zero-lift drag coefficient C D0 =0.0169, induced resistance factor K =0.025.

[0178] The missile's initial position in the guidance coordinate system is x0 = 0m, y0 = 5000m, z0 = 0m; initial velocity is V0 = 300m / s; initial trajectory inclination angle is γ0 = 0°; initial trajectory deflection angle is ψ0 = 0°. The missile's terminal position... Terminal trajectory inclination No-fly zone constraint set to x NFZ =5000m, zNFZ =1000m, R NFZ =2500m.

[0179] The number of discrete nodes in the algorithm is set to N = 50.

[0180] The guidance cycle is set to 0.05s.

[0181] The simulation step size was set to 0.001s.

[0182] The simulation model adds a first-order delay element to the guidance command. The time constant is set to τ = 0.5s.

[0183] The initial conjecture is generated as follows: the initial control variable is set to zero, the initial state variable is generated by integrating the dynamic equation using the initial control variable, and the terminal time is set to the moment when y(t) = 0. This initial conjecture is used in the calculation of the first guidance command in closed-loop guidance; subsequent guidance command calculations all use the result of the previous guidance as the initial conjecture.

[0184] The time t when the no-fly zone constraint is applied i The choice is: using the state variables from the previous iteration, determine the first non-no-fly zone constraint s. k At time t(t)≥0 j The application time satisfies t i ∈[t j -5s,t j That's it.

[0185] Under the aforementioned simulation conditions and parameter settings, guidance calculations were performed using both the MPSP algorithm considering process constraints and the traditional MPSP algorithm. The trajectory, inclination angle, and normal and lateral command acceleration variation curves for the two algorithms are shown below. Figures 1-5 As shown.

[0186] As shown in the figure, both algorithms can satisfy the terminal position and landing angle constraints, and the normal and lateral command accelerations are within reasonable ranges. However, the traditional MPSP algorithm cannot handle no-fly zone constraints, while the algorithm presented in this paper can successfully circumvent no-fly zone constraints.

[0187] From a computational efficiency perspective, the maximum CPU time and average CPU time for guidance command computation using the traditional MPSP algorithm are 8.45ms and 1.39ms, respectively; while the maximum CPU time and average CPU time for guidance command computation using the algorithm proposed in this paper are 11.11ms and 1.86ms, respectively. Comparatively, the algorithm in this paper requires approximately 34% more computation time. Nevertheless, the algorithm in this paper still maintains high computational efficiency and successfully avoids no-fly zone constraints.

Claims

1. A model-based static programming final guidance method for predicting outcomes based on a control barrier function, comprising the following steps: Step 1: Establish the guidance coordinate system and the dynamic equations of a three-degree-of-freedom particle; Step 1-1: Origin of the guidance coordinate system The chosen value is the projection of the missile's center of mass onto the horizontal plane at the moment of shift change. The axis points in the direction of the target. The axis points upwards. The axis is determined according to the right-hand rule; Step 1-2: Establish the three-degree-of-freedom particle dynamics equations for the missile as follows: (1) In the formula: The position vector of the missile. For speed magnitude, For the trajectory inclination angle, For the ballistic deflection angle, For the magnitude of resistance, For missile quality, It is the acceleration due to gravity. and These are the normal and lateral command accelerations, respectively. Steps 1-3: Resistance The calculation formula is: (2) In the formula: For dynamic pressure, For reference area, The drag coefficient; based on the characteristics of the polar curve, the drag coefficient... The calculation formula can be expressed as: (3) In the formula: Zero-lift drag coefficient, Both are induced drag factors and are functions of the velocity magnitude; Steps 1-4: Terminal constraints are expressed as follows: (4) In the formula: The position vector of the target. Let be the desired landing angle; assuming the no-fly zone is an infinitely high cylinder, the no-fly zone constraints are expressed as: (5) In the formula: and Located in the center of the no-fly zone, The radius of the no-fly zone; Steps 1-5: Combining formulas (1)-(5), the missile terminal guidance problem can be described as the following optimal control problem P0: (6) In the formula: For state vectors, For control vectors, This is the compact form of the right-hand side function in the dynamic equation (1). For the specified initial state, This is a compact form of the no-fly zone constraint function. This is a compact form of the terminal constraint function. This refers to the terminal moment of the missile's flight. Step 2: The independent variable is real time. Convert to virtual time The conversion formula is: (7) The optimal control problem P0 can be rewritten as: (8) Step 3: The optimal control problem P1 is discretized using the explicit Euler method. Select A discrete node, uniformly set State vector and control vector Discretizing at the aforementioned nodes, the dynamic differential constraints can be transformed into finite-dimensional algebraic constraints, as follows: (9) In the formula: For the first The state vector at each node For the first Control vectors at each node; For the first The node and the first +1 node distance, similarly, other constraints are also converted into algebraic constraints at discrete nodes; Therefore, the continuous-time optimal control problem P1 is transformed into a finite-dimensional nonlinear programming (NLP) problem, as follows: (10) Step 4: Next, the NLP problem P2 will be further transformed into a constrained nonlinear least squares problem; By simultaneously establishing the dynamic constraints and initial conditions in the NLP problem, intermediate state nodes are gradually eliminated, ensuring that the state vector at any given time is solely the control variable. and terminal time The function; the discrete form of the state trajectory is defined as and control history as The terminal constraints are re-expressed as follows: (11) Therefore, the NLP problem P2 is transformed into the following constrained nonlinear least squares problem: (12) In the formula: control history and terminal time As the variable to be solved; Step 5: Without considering the no-fly zone constraint, the traditional MPSP transforms the nonlinear least squares problem P3 into a series of linear least squares problems for iterative solution; Step 5-1: From the initial conjecture Initially, the iteration format is represented as follows: (13) In the formula: For the iterative sequence, the iteration step size is... , and Terminal constraint functions Regarding the variables to be solved and The Jacobian matrix; Step 5-2: Due to the variable to be solved in formula (13) The dimension is much larger than the constraint function. The number of solutions to the least squares problem is infinite; among them, the least norm solution is the most worthy of study, that is, finding the solution with the smallest iteration step size, such that the difference between the solutions of two iterations is minimized. The least norm least squares problem is constructed as follows: (14) The explicit analytical solution to problem P4 can be obtained using the Lagrange multiplier method, as shown below: (15) Step 5-3: Calculate the control sequence for the next iteration using formulas (13) and (15). and terminal time Then, the system dynamic equations are integrated to recover the state trajectory for the next iteration. Repeat the above steps until the convergence condition is met; Step 6: Using the concept of CBF, transform the no-fly zone constraint into an inequality constraint of the control variable to be solved; Step 6-1: The no-fly zone constraint (5) is re-expressed as: (16) According to the dynamic equation (1), the time derivative of the no-fly zone constraint function can be obtained as follows: (17) Step 6-2: Since formula (17) does not explicitly contain control variables, inequality constraints for control variables cannot be established; the following is given and An approximate calculation formula is used to avoid cases where the control variables are not explicitly included; according to the dynamic equation (26), the position variables... The integral is: (18) By changing the order of integration, formula (46) can be rewritten as: (19) Therefore, the position variable The time derivative is approximately calculated as follows: (20) It can be seen that the approximate calculation formula (48) contains the control variables to be solved. and ; Step 6-3: Same as above, the position variables can be obtained. The approximate time derivative; combining formulas (17) and (20), we can obtain No-fly zone constraint function at time t The approximate time derivative of is expressed in compact form as follows: (21) Where: coefficient , and Both are functions of state variables, represented as follows: (22) Using the forward invariance of CBF, the inequality constraints for the control variables are constructed as follows: (23) Because of the forward invariance of CBF, there is no need to impose no-fly zone constraints at all nodes, but only inequality constraints at one or a few nodes (23); Step 7: Using projection technology, the control inequality constraint (23) is embedded into the solution framework of the traditional MPSP, that is, embedded in step 5, to obtain the MPSP algorithm that can handle the no-fly zone constraint. Step 7-1: In each iteration of MPSP, the state parameters in formula (23) are replaced by the state variables of the previous iteration. This yields the linear inequality constraints on the control variables to be solved, as follows: (24) Where: coefficient , and The formula (22) is calculated using the state variables from the previous iteration; Step 7-2: After obtaining the iterative method and iteration step size using the traditional MPSP, project the unconstrained iterative solution onto the set defined by the linear inequality constraint (24). The iterative format is as follows: (25) In the formula: For the projection operator on the linear inequality constraint (24), the iteration direction is... Similarly, it is calculated using formula (15) in step 5; Step 7-3: Obtain the control variables for the next iteration using formula (25). and terminal time Then, the dynamic equations are integrated to recover the state trajectory for the next iteration. Repeat the above steps until the convergence condition is met.

Images