Probability uncertainty closed-loop robust trajectory planning method considering saturation of executing mechanism

Abstract

The invention relates to a probability uncertainty closed-loop robust trajectory planning method considering execution mechanism saturation, belongs to the technical field of aircraft control, and aims to reduce the sensitivity degree of a guided ammunition scheme trajectory to uncertain factors and give full play to the limited control ability while avoiding execution mechanism saturation as much as possible. Firstly, a control quantity of a deterministic energy optimal trajectory is acquired as an initial value of a subsequent algorithm; secondly, modeling is carried out on main uncertain factors, and an uncertainty quantification module is built based on a PCE method; then, designing a corresponding constraint condition and a target function, establishing an open-loop robust programming model, weighing and considering algorithm precision and calculation time consumption, designing a closed-loop tracking PID controller, and coupling with the open-loop programming model; and finally, carrying out closed-loop robust planning by adopting a direct targeting method. The method fully considers the matching of the trajectory and the guide control system, and overcomes the problems of insufficient description precision of a nonlinear model and sensitivity to an initial value in a traditional method.

Description

Technical Field

[0001] The invention belongs to the field of aircraft control, and in particular relates to a probabilistic uncertain closed-loop robust trajectory planning method considering actuator saturation. Background Art

[0002] Among the many key technologies in the research of guided missiles, trajectory planning has always been a core and hot issue. Limited by cost and the harsh mechanical environment during the launch process, guided munitions are usually small in size and have small aerodynamic control surfaces. There are special requirements for the navigation and control devices that can be used on the projectiles and the overall structure layout. In addition, the projectile is usually unpowered in the controlled trajectory section. Therefore, guided munitions have the characteristics of a small lift-to-drag ratio, limited maneuverability and anti-interference capabilities. Under such conditions, the rationality of the trajectory design is particularly critical.

[0003] At present, most of the research on trajectory design discusses the problem of deterministic optimal control, without considering the influence of uncertain factors, so the designed trajectory is too ideal. In the actual flight of the projectile, there are uncertainties in the dynamic model, aerodynamic and environmental parameters, and the state at the starting point, while the planning design variables, performance indicators, and constraints are closely related to these parameters. The planning results obtained based on the deterministic model may lead to design variables exceeding the range, performance indicators deviating, and constraints breaking the boundaries. This problem can be alleviated through robust planning.

[0004] In the robust planning process, it is crucial to obtain the propagation law of random interference along the trajectory under the action of the current control quantity. The Linear Analysis of Covariance (LAC) method can obtain statistical results similar to Monte Carlo Simulation (MCS) by integrating the error propagation equation once. Due to its high efficiency, it has been widely used in robust planning problems. However, LAC requires linearization of the model, which will lead to large approximation errors when facing problems with a high degree of nonlinearity; at the same time, robust planning based on the LAC method is more sensitive to initial values, because it is shaped based on the reference trajectory, and the robust optimal trajectory is often not within the neighborhood of the reference trajectory, which may cause the algorithm to fall into a local optimum and terminate early.

[0005] The currently available technical achievements are all open-loop robust planning methods, which ignore the coupling between the proposed trajectory and the guidance control system. In actual environments, the deviations caused by various uncertainties are usually eliminated by the guidance control system. Therefore, the matching between the proposed trajectory and the guidance control system has an important impact on the final mission execution effect. For UAVs with sufficient maneuverability, open-loop planning, a processing method with performance loss, is not a big problem. However, for guided munitions with limited rudder control capabilities and unpowered flight, it is very urgent to reduce losses as much as possible and fully tap their maneuverability. Therefore, it is necessary to carry out closed-loop robust trajectory planning that takes into account the feedback of the guidance control system. So far, the discussion on such issues is still insufficient. Summary of the invention

[0006] The technical problems to be solved by the present invention are:

[0007] In order to avoid the shortcomings of the prior art, the present invention provides a probabilistic uncertain closed-loop robust trajectory planning method taking into account the saturation of the actuator, aiming to reduce the sensitivity of the proposed trajectory to uncertain factors and give full play to its limited control capability while avoiding the saturation of the actuator as much as possible.

[0008] In order to solve the above technical problems, the technical solution adopted by the present invention is:

[0009] A probabilistic uncertain closed-loop robust trajectory planning method considering actuator saturation, characterized by comprising the following steps:

[0010] S1: For the simplified three-degree-of-freedom dynamic model of guided munitions, the influence of random factors is ignored, and deterministic trajectory planning is performed with optimal control energy as the objective function. The optimal control quantity is obtained as the initial value of the open-loop robust planning model;

[0011] S2: Consider the main uncertainties in the projectile flight process, and model the deviation from the nominal value as obeying a specific probability distribution to obtain an uncertainty factor model; introduce the uncertainties into the dynamics model to obtain a random dynamics system;

[0012] S3: Based on the polynomial chaos expansion PCE method, the uncertainty of the random dynamic system is quantified to obtain the statistical characteristics of the uncertain factors under the action of the specified control quantity, including:

[0013] S3.1: According to the specific statement of the robust programming problem, if the problem has a certain terminal time, the stochastic dynamics system is retained; if the problem has a free terminal time, the substitution method is used to replace the independent variable from time t with other state variables in the stochastic dynamics system that are monotonically changing and bounded, or to define auxiliary variables;

[0014] S3.2: According to the generalized PCE theory, the WSRSM based on Guass integral node sampling is used to expand the dimension of the stochastic dynamical system and obtain a high-dimensional deterministic dynamical model with PCE coefficients as state variables;

[0015] S3.3: According to the maximum interaction restriction truncation strategy, the high-order cross terms above the second order in the high-dimensional deterministic dynamics model are truncated to obtain a truncated model;

[0016] S3.4: Integrate the truncated model obtained in step S3.3 under the specified control quantity to obtain the variation law of all PCE coefficients with the independent variables, and obtain the statistical characteristics of the state variables in the random dynamic system according to the generalized PCE theory;

[0017] S4: Based on the statistical characteristics of the state variables obtained in step S3, the corresponding constraints and objective functions are designed to establish an open-loop robust programming model; the optimal control problem of the open-loop robust programming model is transformed into a nonlinear programming problem by a direct method to obtain the open-loop robust optimal trajectory of the projectile;

[0018] S5: Design a trajectory tracking algorithm based on the open-loop robust optimal trajectory of step S4, consider the coupling between the proposed trajectory and the guidance control system, and establish a closed-loop robust planning model;

[0019] S6: Deviation of the uncertainty factors of the uncertainty factor model in step S2, using the trajectory tracking algorithm designed in step S5 to track the open-loop robust optimal trajectory obtained in step S4 to obtain the closed-loop control quantity of the projectile, introducing the closed-loop control energy loss into the objective function in a weighted form, and iteratively obtaining the closed-loop robust optimal trajectory of the projectile.

[0020] A further technical solution of the present invention: In step S1, for a simplified three-degree-of-freedom dynamic model of guided ammunition:

[0021]

[0022] In the formula, F x 、F y 、F z are the total drag, total lift and total lateral force acting on the projectile, m is the mass of the projectile, g is the acceleration of gravity, θ and ψ are the ballistic inclination and ballistic deviation, x, y and z are the range, height and lateral deviation of the projectile, ρ is the air density, V is the projectile velocity, q = ρV 2 / 2 is the dynamic pressure, S is the characteristic area of ​​the projectile, C x0 is the zero-lift drag coefficient, k c is the total induced drag coefficient, are the partial derivatives of the total lift coefficient and the total side force coefficient with respect to the quasi-attack angle α and the quasi-sideslip angle β, respectively. For the axisymmetric shape,

[0023] Ignoring the influence of random factors, deterministic trajectory planning is performed with optimal control energy as the objective function:

[0024]

[0025] Where b and p are the corresponding boundary constraints and path constraints, t0 and t f are the initial and final moments respectively;

[0026] The sequential quadratic programming method is used to solve the above optimal control problem and obtain the optimal control quantity U * As the initial value of the robust programming model.

[0027] A further technical solution of the present invention: In step S2, the random factors in the real flight environment are abstracted into aerodynamic parameters, air density and start-up control point states, and are modeled as deviations between real values ​​and nominal values, and are assumed to obey mutually independent normal distributions:

[0028]

[0029] In the formula, They represent the true values ​​of the zero-lift drag coefficient, the derivative of the total lift coefficient, and the derivative of the total lateral force coefficient, respectively. d 、N l 、N l′ is a random number that follows a standard normal distribution N(0,1), σ d , σ l , σ l′ Respectively represent the degree to which the above parameters deviate from the reference value, N ρ is a random number that follows the standard normal distribution N(0,1), N V 、N θ 、N ψ 、N x 、N y 、N z is a random number that follows a standard normal distribution N(0,1), σ V , σ θ , σ ψ , σ x , σ y , σ z Respectively represent the standard deviation of each state at the starting control point, V0, θ0, ψ0, t0, y0, z0 are the initial values ​​of projectile velocity, trajectory inclination, trajectory deviation, time, projectile height, and projectile side deviation, respectively. is the true value of air density taking into account the uncertainty, are the true values ​​of the starting control point states considering uncertain factors.

[0030] A further technical solution of the present invention: In the step S3.1, the terminal flight moment is considered to be free, and the independent variable of the dynamic model is converted into the range x:

[0031]

[0032] In step S3.2, according to the generalized PCE theory, the WSRSM based on Guass integral node sampling is used to expand the dimension of the dynamic model to obtain a high-dimensional deterministic dynamic model with PCE coefficients as state variables:

[0033]

[0034] In the formula, S PCE ∈R 396 is the PCE coefficient of all state variables, P PCE ∈R 8 is the PCE coefficient of the random parameter; H = (Ψ T W s Ψ) -1 Ψ T W∈R 396×792 is the regression coefficient matrix, W s ∈R 792×792 =diag(ω s ) is the weight matrix;

[0035] In step S3.3, according to the maximum interaction limit truncation strategy:

[0036]

[0037] In the formula, λ = [λ1, λ2,..., λ d ] is a multi-index, A d，p is a set of multi-indexes, and ||λ||0 represents the rank of the multi-index;

[0038] Truncating the high-order cross terms above the second order in the full-order PCE model to reduce the amount of model calculation;

[0039] In step S3.4, the truncated model obtained in step S3.3 is integrated under the specified control amount to obtain the variation law of all PCE coefficients with the independent variables, and the statistical characteristics of the state variables in the original random dynamics model are obtained according to the generalized PCE theory:

[0040]

[0041] Where ξ∈R 10 is a d-dimensional cell array, Φ j (ξ) is the value of φ j (ξ k ) The orthogonal polynomial function obtained by taking the tensor product, φj (ξ k ), k=1,2,...,d are the one-dimensional orthogonal basis functions corresponding to the uncertain factors of each dimension, S ij is the j-th PCE coefficient of the ith state variable.

[0042] A further technical solution of the present invention: In step S4, corresponding constraints and objective functions are designed to establish an open-loop robust programming model:

[0043]

[0044] J o =D(y f )+D(z f ) (13)

[0045] In the formula, α max , β max are the maximum values ​​of the quasi-attack angle α and the quasi-sideslip angle β, V fmin is the minimum value of the terminal velocity, y T 、z T are the y-axis and z-axis coordinate values ​​of the target position respectively; the operator E represents the mean, and the operator D represents the variance;

[0046] The direct shooting method is used to transform the optimal control problem into a nonlinear programming problem and obtain the open-loop robust optimal trajectory of the projectile.

[0047] A further technical solution of the present invention: In step S5, a trajectory tracking algorithm is designed:

[0048]

[0049] In the formula, and are the height and lateral deviation at range x on the reference trajectory, is the open-loop control variable, k P , k I , k D Respectively represent the proportional parameter, integral parameter and differential parameter, α c , β c are the control values ​​of quasi-attack angle α and quasi-sideslip angle β respectively.

[0050] A further technical solution of the present invention is as follows: in step S6, the closed-loop control energy loss is introduced into the objective function in a weighted form, and the closed-loop robust optimal trajectory of the projectile is iteratively obtained:

[0051]

[0052] In the formula, ω D is the weight of the terminal position distribution in the closed-loop performance index.

[0053] A computer system, characterized in that it comprises: one or more processors, and a computer-readable storage medium for storing one or more programs, wherein when the one or more programs are executed by the one or more processors, the one or more processors implement the above-mentioned method.

[0054] A computer-readable storage medium, characterized in that it stores computer-executable instructions, which are used to implement the above method when executed.

[0055] A computer program product, characterized by comprising computer executable instructions, wherein the instructions are used to implement the above method when executed.

[0056] The beneficial effects of the present invention are:

[0057] The present invention provides a probabilistic uncertain closed-loop robust trajectory planning method considering actuator saturation, which has the following advantages over the prior art:

[0058] 1. The present invention adopts an accurate and efficient PCE method to quantify uncertainty, overcoming the problems of insufficient accuracy in describing nonlinear models and sensitivity to initial values ​​in the traditional LAC method. At the same time, the method of the present invention performs maximum interaction restriction basis truncation on the full-order PCE model, which greatly reduces the dimension of the problem and the computational time while ensuring the accuracy of the model.

[0059] 2. The present invention adopts the WSRSM transformation dynamics model based on Guass integral node sampling. When the number of sampling points is determined, the Guass integral nodes and corresponding weights are uniquely determined, avoiding the additional randomness introduced by traditional Latin Hypercube Design (LHD) sampling; at the same time, unlike traditional SRSM, WSRSM fully considers the differences in the degree of influence of sample points on the regression process, while improving the accuracy of the PCE model, it does not increase substantial computational overhead.

[0060] 3. The closed-loop robust trajectory planning method proposed in the present invention couples the planning of the proposed trajectory with the tracking control, fully considering the limited steering control capability of the guided munition, and avoiding the saturation of the actuator as much as possible while reducing the sensitivity of the proposed trajectory to uncertain factors. BRIEF DESCRIPTION OF THE DRAWINGS

[0061] The drawings are only for the purpose of illustrating particular embodiments and are not to be considered limiting of the present invention. Like reference symbols denote like components throughout the drawings.

[0062] Figure 1It is a flow chart of the probabilistic uncertain closed-loop robust trajectory planning method considering actuator saturation provided by the present invention.

[0063] Figure 2 It is a flow chart of the probabilistic uncertain closed-loop robust trajectory planning method considering actuator saturation of the present invention. DETAILED DESCRIPTION

[0064] In order to make the purpose, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention. In addition, the technical features involved in the various embodiments of the present invention described below can be combined with each other as long as they do not conflict with each other.

[0065] The present invention provides a probabilistic uncertain closed-loop robust trajectory planning method considering actuator saturation, such as Figure 1 As shown, the following steps are included:

[0066] S1: For the simplified three-degree-of-freedom dynamic model of guided munitions, the influence of random factors is ignored, and deterministic trajectory planning is performed with optimal control energy as the objective function. The optimal control quantity is obtained as the initial value of the subsequent robust planning algorithm.

[0067] S2: Consider the main uncertainties in the projectile flight process, and model their deviations from the nominal value as obeying a specific probability distribution; introduce the uncertainties into the dynamics model to obtain a random dynamics system. For uncertainties that are not suitable for probability description, they are regarded as unmodeled interference and compensated by the guidance and control system.

[0068] S3: Considering the uncertainty factor model established in step S2, based on the Polynomial Chaos Expansion (PCE) method, the uncertainty of the random dynamic system is quantified to obtain its statistical characteristics under the action of the specified control quantity, mainly including:

[0069] S3.1: According to the specific statement of the robust programming problem, if the terminal time of the problem is determined, the stochastic dynamics system is retained; if the terminal time of the problem is free, the substitution method is used to replace the independent variable from time t with other state variables that are monotonically changing and bounded in the stochastic dynamics system or define auxiliary variables.

[0070] S3.2: According to the generalized PCE theory, the weighted stochastic response surface method (WSRSM) based on Guass integral node sampling is used to expand the dimension of the random dynamic system and obtain a high-dimensional deterministic dynamic model with PCE coefficients as state variables.

[0071] S3.3: According to the maximum interaction limit truncation strategy, high-order cross terms with smaller effects (above the second order) in the full-order PCE model are truncated to reduce the amount of model calculation.

[0072] S3.4: Integrate the truncated model obtained in step S3.3 under the specified control quantity to obtain the variation law of all PCE coefficients with the independent variables, and obtain the statistical characteristics of the state variables in the random dynamic system based on the generalized PCE theory.

[0073] S4: Based on step S3, design corresponding constraints and objective functions in combination with problem requirements, and establish an open-loop robust programming model. Use the direct method to transform the optimal control problem into a nonlinear programming problem and obtain the open-loop robust optimal trajectory of the projectile.

[0074] S5: Weighing accuracy and computational cost, a simple and effective trajectory tracking algorithm is designed. Based on step S4, the coupling between the trajectory and the guidance control system is considered to establish a closed-loop robust planning model.

[0075] S6: The uncertain factors modeled in step S2 are pulled aside to simulate a relatively harsh real environment. The algorithm designed in step S5 is used to track the open-loop robust optimal trajectory obtained in step S4 to obtain the closed-loop control quantity of the projectile. The closed-loop control energy loss is introduced into the objective function in a weighted form, and the closed-loop robust optimal trajectory of the projectile is iteratively obtained.

[0076] In order to enable those skilled in the art to better understand the present invention, the present invention is described in detail below in conjunction with specific embodiments.

[0077] Example 1

[0078] The embodiment of the present invention provides a probabilistic uncertain closed-loop robust trajectory planning method considering actuator saturation, and the specific steps are as follows:

[0079] S1: For a simplified three-degree-of-freedom dynamic model of guided munitions:

[0080]

[0081] In the formula, F x 、F y 、F zare the total drag, total lift and total lateral force acting on the projectile, m is the mass of the projectile, g is the acceleration of gravity, θ and ψ are the ballistic inclination and ballistic deviation, x, y and z are the range, height and lateral deviation of the projectile, ρ is the air density, V is the projectile velocity, q = ρV 2 / 2 is the dynamic pressure, S is the characteristic area of ​​the projectile, C x0 is the zero-lift drag coefficient, k c is the total induced drag coefficient, are the partial derivatives of the total lift coefficient and the total side force coefficient with respect to the quasi-attack angle α and the quasi-sideslip angle β, respectively. For the axisymmetric shape,

[0082] Ignoring the influence of random factors, deterministic trajectory planning is performed with optimal control energy as the objective function:

[0083] Let the combination of state variables be X = [V, θ, ψ, x, y, z] T , the combination of control variables is U = [α, β] T , then formula (17) can be expressed as:

[0084]

[0085] Where t is time.

[0086] The sequential quadratic programming method is used to solve the optimal control problem (19) and the optimal control quantity U is obtained. * As the initial value of the subsequent robust planning algorithm.

[0087]

[0088] Where b and p are the corresponding boundary constraints and path constraints, t0 and t f are the initial and final moments respectively.

[0089] S2: Consider the main uncertainties in the projectile flight process and model their deviations from the nominal value as obeying a specific probability distribution. Uncertainties that are not suitable for probability description are regarded as unmodeled interference and compensated by the guidance and control system.

[0090] S2.1: Aerodynamic parameter deviation

[0091] Taking into account the differences in the machining errors of the projectile and the ballistic environment, it is assumed that the deviations of the aerodynamic parameters obey independent normal distributions:

[0092]

[0093] In the formula, They represent the true values ​​of the zero-lift drag coefficient, the derivative of the total lift coefficient, and the derivative of the total lateral force coefficient, respectively.d 、N l 、N l′ is a random number that follows a standard normal distribution N(0,1), σ d , σ l , σ l′ They respectively indicate the degree to which the above parameters deviate from the reference values.

[0094] S2.2: Meteorological environment deviation

[0095] The difference between real weather and standard weather is mainly reflected in the deviation of wind field, temperature and atmospheric density. Among them, the influence of wind field on trajectory calculation is not suitable for description by probability, and can be regarded as an unmodeled random force, which is compensated by the guidance and control system; the temperature mainly affects the calculation of the speed of sound, and then affects the flight Mach number, and the change of Mach number will lead to the change of aerodynamic parameters. Therefore, the deviation of temperature can be attributed to the deviation of aerodynamic parameters. In summary, the present invention models the deviation of meteorological environment as the deviation of atmospheric density. The statistical characteristics of the deviation of real atmospheric density from the standard value are related to the altitude y, and the ratio of its standard deviation to the standard atmospheric density can be approximated as:

[0096]

[0097] Therefore, the true atmospheric density can be expressed as

[0098]

[0099] Where N ρ is a random number that follows a standard normal distribution N(0,1).

[0100] S2.3: Start point state deviation

[0101] After the projectile is launched, under the influence of various uncertain factors, the state value at the starting control point will inevitably deviate from the design value of the trajectory. The present invention assumes that the state deviation obeys a mutually independent normal distribution, and the design value of the plan at the starting control point is X0 = [V0, θ0, ψ0, x0, y0, z0] T , then the real state at the starting control point can be expressed as:

[0102]

[0103] Where N V 、N θ 、N ψ 、N x 、N y 、N z is a random number that follows a standard normal distribution N(0,1), σ V , σ θ , σ ψ , σ x, σ y , σ z Respectively represent the standard deviation of each state at the starting control point.

[0104] S3: Considering the uncertainty factor model established in step S2, based on the Polynomial Chaos Expansion (PCE) method, the uncertainty of the random dynamic system is quantified to obtain its statistical characteristics under the action of the specified control quantity, including:

[0105] S3.1: According to the specific statement of the robust programming problem, if the terminal time of the problem is determined, the original dynamic model is retained; if the terminal time of the problem is free, the substitution method is used to replace the independent variable from time t to other state variables in the system that are monotonically changing and bounded, or to define auxiliary variables.

[0106] For general considerations, it is assumed that the terminal attack time of the trajectory is free and unknown, and the position of the target to be hit is fixed. During the flight, the range x of the projectile is monotonically increasing and bounded. Therefore, the dynamic model equation (17) can be transformed into:

[0107]

[0108] The combination of system state variables after transformation is S = [V, θ, ψ, t, y, z] T , the combination of control variables is still U = [α, β] T , then formula (24) can be expressed as:

[0109]

[0110] Similarly, formula (23) can be transformed into:

[0111]

[0112] S3.2: According to the generalized PCE theory, the weighted stochastic response surface method (WSRSM) based on Guass integral node sampling is used to expand the dimension of the stochastic ballistic dynamics equations to obtain a high-dimensional deterministic dynamic model with PCE coefficients as state variables.

[0113] The main uncertainties in the projectile flight process can be attributed to two categories: parameter uncertainty and initial state uncertainty. P =4, the dimension d of the initial state is uncertain S =6, so the uncertainty factor dimension of the entire random dynamic system is d = 10. The random parameter is P = [σ d ,σ l,σ l′ ,σ ρ ] T And introduce it into the dynamic system shown in formula (25), then the stochastic dynamic model can be expressed as:

[0114]

[0115] The non-intrusive PCE method based on Guass integral node sampling and WSRSM is used to transform Equation (27) into a high-dimensional deterministic dynamic model with PCE coefficients as state variables. Taking the PCE order p = 2, the number of PCE coefficients corresponding to each state variable after expansion is:

[0116]

[0117] The number of sampling points used for regression is:

[0118] n s =2(P+1)=132 (29)

[0119] According to the generalized chaos polynomial theory, each state variable in equation (27) can be expanded into a polynomial form:

[0120]

[0121] Where ξ∈R 10 is a d-dimensional cell array, where each element is the n corresponding to the uncertainty factor s Gaussian integral sampling point combination. The one-dimensional orthogonal basis function corresponding to each dimension of uncertainty is φ j (ξ k ),k=1,2,...,d, then Φ j (ξ) is the value of φ j (ξ k ) is the orthogonal polynomial function obtained by taking the tensor product. S ij is the jth PCE coefficient of the ith state variable, and its total number is:

[0122] n PCE =d S (P+1)=396 (31)

[0123] For random parameters that follow a normal distribution, a 1-order, 1-dimensional orthogonal polynomial model can accurately describe its uncertainty. Therefore, each random parameter can be expanded into the form of a sum of two terms:

[0124] P i =p i0 Φ0(ξ P )+p i1 Φ1(ξ P),i=1,2,...,d P (32)

[0125] In the formula, ξ P ∈R 4 is d P dimensional cell array, where each element is the n-dimensional cell array corresponding to the uncertainty factor. s Gaussian integral sampling point combination.

[0126] At this point, the random system shown in formula (27) is transformed into a high-dimensional deterministic system:

[0127]

[0128] In the formula, S PCE ∈R 396 is the PCE coefficient of all state variables, P PCE ∈R 8 is the PCE coefficient of the random parameter, F∈R 792 To n s Substitute the sample points into the true response value of the original random system formula (27), Ψ∈R 792×396 is the coefficient matrix of the transformed deterministic system, which is described in detail as follows:

[0129] Let Ω(ξ i )=[Φ0(ξ i ),Φ1(ξ i ),...Φ P (ξ i )],i=1,2,...,n s , define the intermediate matrix:

[0130]

[0131] In the formula, the symbol represents the Kronecker product, then:

[0132]

[0133] Taking into account the different importance of each sample point, the weight vector is recorded as:

[0134]

[0135] In the formula, ω is n s The weight corresponding to each sample point.

[0136] At this point, according to WSRSM, the overdetermined system of equations shown in equation (17) can be solved by weighted least quadratic regression:

[0137]

[0138] In the formula, H=(Ψ T W s Ψ) -1 Ψ T W∈R 396×792 is the regression coefficient matrix, W s ∈R 792×792 =diag(ω s ) is the weight matrix.

[0139] S3.3: According to the maximum interaction limit truncation strategy, high-order cross terms with smaller effects (above the second order) in the full-order PCE model are truncated to reduce the amount of model calculation.

[0140] The number of state variables of the deterministic system after transformation is as high as 396, and the number of Gaussian integral sampling points for regression is also as high as 132. Under this problem scale, the computational time of a single uncertainty propagation is barely acceptable. However, in the subsequent robust planning of the glide trajectory, the iterative process will repeatedly call the uncertainty propagation module, and the convergence speed of the optimization algorithm will decrease significantly when the number of variables is too large, which will lead to a significant reduction in computational efficiency. It is known that with the increase of the dimension d of the uncertainty factor, the total number of terms in the orthogonal polynomial increases exponentially, but in many problems, the terms of the orthogonal polynomial are not equally important. Usually, the important terms in the expansion often only involve a few variables. In most actual engineering systems, the output response is mainly affected by the input variables of each dimension and their low-order cross terms, and the influence of high-order cross terms is relatively small. Based on the above analysis, the present invention adopts the maximum interaction limit truncation strategy to directly remove high-order polynomials, thereby achieving the purpose of reducing the amount of calculation.

[0141] In the full-order PCE model (30), a multivariate orthogonal polynomial is constructed by direct tensor product. When the PCE order is p, it needs to satisfy:

[0142]

[0143] Where λ=[λ1,λ2,...,λ d ] is a multi-index, A d,p is a collection of multiple indices, λ i is the order of the univariate orthogonal polynomial basis corresponding to the j-th dimension of the multivariate orthogonal polynomial, and the operator |·|1 represents the 1-norm.

[0144] The main idea of ​​the maximum interaction restriction truncation is: from the multi-index set A defined by formula (38) d,p Select low rank λ to form a subset A d,p,r , so that λ not only satisfies ||λ||1≤p, but also satisfies ||λ||0≤r.

[0145]

[0146] Where ||λ||0 represents the rank of the multi-index.

[0147] For the present invention, d=10, p=2. Therefore, according to formula (23), r=1 is taken. After removing the high-order cross terms, the number of PCE coefficients P+1 of each dimensional random variable is reduced from 66 in the full-order to 11. Therefore, after the dimensional expansion transformation, the total number of state variables of the deterministic dynamic system is reduced from 396 to 66, and the number of Gaussian integral sampling points used for regression is reduced from 132 to 22, a reduction of 83.3%.

[0148] S3.4: Integrate the truncated model obtained in step S3.3 under the specified control quantity to obtain the variation law of all PCE coefficients with the independent variables, and obtain the statistical characteristics of the state variables in the original random dynamics model based on the generalized PCE theory.

[0149] By numerically integrating equation (37), we can obtain the variation law of all PCE coefficients with range x. The mean μ and standard deviation σ of each state variable in the original random system can be directly obtained through the PCE coefficient:

[0150]

[0151] S4: Based on step S3, design corresponding constraints and objective functions in combination with problem requirements, and establish an open-loop robust programming model. Use the direct method to transform the optimal control problem into a nonlinear programming problem and obtain the open-loop robust optimal trajectory of the projectile.

[0152] S4.1: Constraints

[0153] After considering the influence of uncertain factors, there are countless (bounded) flight trajectories under the same open-loop control input. Therefore, the focus of open-loop robust planning is not the specific value of the state, but its statistical law, and the statistical moments of state variables should also be considered in the constraints.

[0154] S4.1.1: Boundary constraints

[0155] The state deviations at the starting control points follow mutually independent normal distributions, so the constraint on the initial state mean is:

[0156]

[0157] In order to ensure the final strike accuracy and strike effect, the constraint on the terminal state mean is:

[0158]

[0159] S4.1.2: Process Constraints

[0160] Considering the limited rudder control capability of the projectile, the control quantity constraint during flight is:

[0161]

[0162] S4.2: Performance Indicators

[0163] For the trajectory of the program, the dispersion of the terminal position is usually used as a measure of its sensitivity to uncertainty, while less attention is paid to the dispersion of other terminal states and the state envelope during the projectile flight. Therefore, the performance index of the open-loop robust planning is designed as:

[0164] J o =D(y f )+D(z f ) (44)

[0165] In the above formulas, the operator E represents the mean, and the operator D represents the variance.

[0166] S5: Weighing accuracy and computational cost, a simple and effective trajectory tracking algorithm is designed. Based on step S4, the coupling between the trajectory and the guidance control system is considered to establish a closed-loop robust planning model.

[0167] Considering that the guidance module will be called repeatedly during the iteration of the closed-loop planning algorithm, a sufficiently simple and effective trajectory tracking algorithm is required. The traditional LQR algorithm needs to linearize the model near the reference trajectory. When the deviation is large, a large tracking error will occur. At the same time, the LQR algorithm needs to solve the Riccati equation, which increases the difficulty of solving the algorithm. The present invention uses a PID controller with an analytical expression and is independent of the model to track the trajectory of the solution. Let the height and side deviation at the range x on the reference trajectory be and The control amount is Then the position error between the projectile and the reference trajectory can be expressed as:

[0168]

[0169] The closed-loop control quantity after feedback correction can be expressed as:

[0170]

[0171] In the formula, k P , k I , k D They represent the proportional parameter, integral parameter and differential parameter respectively.

[0172] The dynamic model and related constraints of closed-loop robust planning are the same as those in open-loop planning and will not be described here.

[0173] The closed-loop control energy loss during the projectile flight is recorded as:

[0174]

[0175] The performance index of closed-loop robust planning can be designed as:

[0176] J c =ω D J o +E c (48)

[0177] In the formula, ω D is the weight of the terminal position distribution in the closed-loop performance index.

[0178] From equation (48), we can see that closed-loop robust trajectory planning is a multi-objective optimal control problem. Considering J o With E c The order of magnitude of the difference is large. In order to balance the planning emphasis, the present invention takes the deterministic energy optimal trajectory without considering random interference as the benchmark trajectory and non-dimensionalizes the performance index. The terminal position dispersion of the benchmark trajectory after uncertainty propagation is Control energy loss is The performance index is converted into:

[0179]

[0180] S6: The uncertain factors modeled in step S2 are pulled aside to simulate a relatively harsh real environment. The algorithm designed in step S5 is used to track the open-loop robust optimal trajectory obtained in step S4 to obtain the closed-loop control quantity of the projectile. The closed-loop control energy loss is introduced into the objective function in a weighted form, and the closed-loop robust optimal trajectory of the projectile is iteratively obtained.

[0181] like Figure 2 As shown, considering the situation where all uncertain factors are pulled to the half limit (the ±1.5σ value of each random parameter is taken in the present invention), tracking the robust optimal trajectory can obtain the closed-loop feedback of the control amount, and then obtain the control energy loss during the projectile flight (according to simulation experience and measured data, the situation where all uncertain factors are pulled to the half limit is already quite harsh. At this time, trajectory tracking is relatively difficult, and the control energy consumption is also large. In the planning process, such a conservative processing method can basically ensure the smooth progress of closed-loop guidance during the actual flight of the projectile). By introducing the control energy loss as a penalty function into the planning performance index, the coupling and matching characteristics of the solution trajectory and the guidance system can be considered in the iterative process, and a closed-loop robust optimal trajectory that takes into account both robustness and feasibility can be planned.

[0182] The above description is only a specific implementation mode of the present invention, but the protection scope of the present invention is not limited thereto. Any technician familiar with the technical field can easily think of various equivalent modifications or substitutions within the technical scope disclosed in the present invention, and these modifications or substitutions should be included in the protection scope of the present invention.

Claims

1. A probabilistic uncertain closed-loop robust trajectory planning method considering actuator saturation, characterized in that: The following steps are involved: S1: For the simplified three-degree-of-freedom dynamic model of guided munitions, the influence of random factors is ignored, and deterministic trajectory planning is performed with optimal control energy as the objective function. The optimal control quantity is obtained as the initial value of the open-loop robust planning model; S2: Consider the main uncertainties in the projectile flight process and model their deviations from the nominal value as obeying a specific probability distribution; introduce the uncertainties into the dynamic model to obtain a random dynamic system; S3: Based on the polynomial chaos expansion PCE method, the uncertainty of the random dynamic system is quantified to obtain the statistical characteristics of the uncertain factors under the action of the specified control quantity, including: S3.1: According to the specific statement of the robust programming problem, if the problem has a certain terminal time, the stochastic dynamics system is retained; if the problem has a free terminal time, the substitution method is used to replace the independent variable from time t with other state variables in the stochastic dynamics system that are monotonically changing and bounded, or to define auxiliary variables; S3.2: According to the generalized PCE theory, the WSRSM based on Guass integral node sampling is used to expand the dimension of the stochastic dynamical system and obtain a high-dimensional deterministic dynamical model with PCE coefficients as state variables; S3.3: According to the maximum interaction restriction truncation strategy, the high-order cross terms above the second order in the high-dimensional deterministic dynamics model are truncated to obtain a truncated model; S3.4: Integrate the truncated model obtained in step S3.3 under the specified control quantity to obtain the variation law of all PCE coefficients with the independent variables, and obtain the statistical characteristics of the state variables in the random dynamic system according to the generalized PCE theory; S4: Based on the statistical characteristics of the state variables obtained in step S3, the corresponding constraints and objective functions are designed to establish an open-loop robust programming model; the optimal control problem of the open-loop robust programming model is transformed into a nonlinear programming problem by a direct method to obtain the open-loop robust optimal trajectory of the projectile; S5: Design a trajectory tracking algorithm based on the open-loop robust optimal trajectory of step S4, consider the coupling between the proposed trajectory and the guidance control system, and establish a closed-loop robust planning model; S6: Deviation of the uncertainty factors of the uncertainty factor model in step S2, using the trajectory tracking algorithm designed in step S5 to track the open-loop robust optimal trajectory obtained in step S4 to obtain the closed-loop control quantity of the projectile, introducing the closed-loop control energy loss into the objective function in a weighted form, and iteratively obtaining the closed-loop robust optimal trajectory of the projectile.

2. According to claim 1, a probabilistic uncertain closed-loop robust trajectory planning method considering actuator saturation is characterized in that: In step S1, for the simplified three-degree-of-freedom dynamics model of guided munitions: In the formula, F x 、F y 、F z are the total drag, total lift and total lateral force acting on the projectile, m is the mass of the projectile, g is the acceleration of gravity, θ and ψ are the ballistic inclination and ballistic deviation, x, y and z are the range, height and lateral deviation of the projectile, ρ is the air density, V is the projectile velocity, q = ρV 2 / 2 is the dynamic pressure, S is the characteristic area of ​​the projectile, C x0 is the zero-lift drag coefficient, k c is the total induced drag coefficient, are the partial derivatives of the total lift coefficient and the total side force coefficient with respect to the quasi-attack angle α and the quasi-sideslip angle β, respectively. For the axisymmetric shape, Ignoring the influence of random factors, deterministic trajectory planning is performed with optimal control energy as the objective function: Where b and p are the corresponding boundary constraints and path constraints, t0 and t f are the initial and final moments respectively; The sequential quadratic programming method is used to solve the above optimal control problem and obtain the optimal control quantity U * As the initial value of the robust programming model.

3. According to claim 2, a probabilistic uncertain closed-loop robust trajectory planning method considering actuator saturation is characterized in that: In step S2, the random factors in the real flight environment are abstracted into aerodynamic parameters, air density and start-up control point states, and are modeled as deviations between real values ​​and nominal values, and are assumed to obey mutually independent normal distributions: In the formula, They represent the true values ​​of the zero-lift drag coefficient, the derivative of the total lift coefficient, and the derivative of the total lateral force coefficient, respectively. d 、N l 、N l′ is a random number that follows a standard normal distribution N(0,1), σ d , σ l , σ l′ Respectively represent the degree to which the above parameters deviate from the reference value, N ρ is a random number that follows the standard normal distribution N(0,1), N V 、N θ 、N ψ 、N x 、N y 、N z is a random number that follows a standard normal distribution N(0,1), σ V , σ θ , σ ψ , σ x , σ y , σ z Respectively represent the standard deviation of each state at the starting control point, V0, θ0, ψ0, t0, y0, z0 are the initial values ​​of projectile velocity, trajectory inclination, trajectory deviation, time, projectile height, and projectile side deviation, respectively. is the true value of air density taking into account the uncertainty, are the true values ​​of the starting control point states considering uncertain factors.

4. According to claim 3, a probabilistic uncertain closed-loop robust trajectory planning method considering actuator saturation is characterized in that: In step S3.1, the terminal flight moment is considered to be free, and the independent variable of the dynamic model is converted into the range x: In step S3.2, according to the generalized PCE theory, the WSRSM based on Guass integral node sampling is used to expand the dimension of the dynamic model to obtain a high-dimensional deterministic dynamic model with PCE coefficients as state variables: In the formula, S PCE ∈R 396 is the PCE coefficient of all state variables, P PCE ∈R 8 is the PCE coefficient of the random parameter; H = (Ψ T W s Ψ) -1 Ψ T W∈R 396×792 is the regression coefficient matrix, W s ∈R 792×792 =diag(ω s ) is the weight matrix; In step S3.3, according to the maximum interaction limit truncation strategy: Where λ=[λ1,λ2,...,λ d ] is a multi-index, A d,p is a set of multi-indexes, and ||λ||0 represents the rank of the multi-index; Truncating the high-order cross terms above the second order in the full-order PCE model to reduce the amount of model calculation; In step S3.4, the truncated model obtained in step S3.3 is integrated under the specified control amount to obtain the variation law of all PCE coefficients with the independent variables, and the statistical characteristics of the state variables in the original random dynamics model are obtained according to the generalized PCE theory: Where ξ∈R 10 is a d-dimensional cell array, Φ j (ξ) is the value of φ j (ξ k ) The orthogonal polynomial function obtained by tensor product, φ j (ξ k ), k=1,2,...,d are the one-dimensional orthogonal basis functions corresponding to the uncertain factors of each dimension, S ij is the j-th PCE coefficient of the ith state variable.

5. According to claim 4, a probabilistic uncertain closed-loop robust trajectory planning method considering actuator saturation is characterized in that: In step S4, corresponding constraints and objective functions are designed to establish an open-loop robust programming model: min J o =D(y f )+D(z f ) (13) In the formula, α max , β max are the maximum values ​​of the quasi-attack angle α and the quasi-sideslip angle β, V fmin is the minimum value of the terminal velocity, y T 、z T are the y-axis and z-axis coordinate values ​​of the target position respectively; the operator E represents the mean, and the operator D represents the variance; The direct shooting method is used to transform the optimal control problem into a nonlinear programming problem and obtain the open-loop robust optimal trajectory of the projectile.

6. The method for probabilistic uncertain closed-loop robust trajectory planning considering actuator saturation according to claim 5, characterized in that: In step S5, a trajectory tracking algorithm is designed: In the formula, and are the height and lateral deviation at range x on the reference trajectory, is the open-loop control variable, k P , k I , k D Respectively represent the proportional parameter, integral parameter and differential parameter, α c , β c are the control values ​​of quasi-attack angle α and quasi-sideslip angle β respectively.

7. The method of probabilistic uncertain closed-loop robust trajectory planning considering actuator saturation according to claim 6, characterized in that: In step S6, the closed-loop control energy loss is introduced into the objective function in a weighted form, and the closed-loop robust optimal trajectory of the projectile is obtained iteratively: In the formula, ω D is the weight of the terminal position distribution in the closed-loop performance index.

8. A computer system, characterized in that include: One or more processors, and a computer-readable storage medium for storing one or more programs, wherein when the one or more programs are executed by the one or more processors, the one or more processors are enabled to implement the method of claim 1.

9. A computer-readable storage medium, characterized in that Computer executable instructions are stored, and when the instructions are executed, they are used to implement the method of claim 1.

10. A computer program product, characterized in that The invention comprises computer executable instructions, which are used to implement the method of claim 1 when being executed.

Images