CN111783232B - An adaptive optimization method for the return segment of a recoverable rocket based on cluster analysis - Google Patents

Abstract

The invention discloses an adaptive optimization method for the return segment of a recoverable rocket based on cluster analysis, comprising: establishing a return trajectory optimization problem model of a sub-stage return segment; solving the state transition matrix of terminal constraint variables relative to design variables; The data elements in the state transition matrix are divided into sensitivity by clustering algorithm, and the data elements with higher sensitivity are screened out to form a state transition matrix M; according to the set linearity criterion, linear constraint variables, linear The design variables and the corresponding linear state transition matrix are selected from M to screen the nonlinear constraint variables and nonlinear design variables; for the linear parameters and nonlinear parameters, different optimization methods are used to obtain the optimal design variable values. The invention can effectively deal with the complex landing constraint problem, has high precision, convergence and calculation efficiency, reduces the dimension and complexity of the parameter space of the optimization problem in the return segment, and has good adaptability and expansibility.

Description

An adaptive optimization method for the return segment of a recoverable rocket based on cluster analysis

technical field

The invention belongs to the field of spacecraft flight control, in particular to an adaptive optimization method for the return segment of a recoverable rocket based on cluster analysis.

Background technique

In recent years, with the advent of the era of commercial spaceflight, the vertical take-off and landing technology of retrievable launch vehicles has become a research hotspot in the field of international spaceflight. At present, vertical take-off and landing technology has become a feasible and most promising technology for recoverable launch vehicles, which greatly reduces the cost of spacecraft launch missions. As the core key technology of vertical take-off and landing, the optimal trajectory design of the first sub-stage return section of the entire launch vehicle (hereinafter referred to as the return section) is a systematic and holistic trajectory optimization problem. There is a mutual coupling and cross-linking relationship between them, and the requirements for landing accuracy are very high.

The return segment is an important basis for the orbital design of the entire probe and the overall design of the rocket. It includes multiple flight segments. The reentry speed is high, the return process is complicated, the range and speed span are large, the flight environment changes significantly, and various factors are strong disturbance factors. . The trajectory optimization of the return segment must not only ensure high-precision and stable landing at the terminal, but also comprehensively consider the dynamic pressure, overload, heat flux density, engine performance, and terminal landing position, speed, attitude and other constraints of each flight segment. It is necessary to carry out related research on the nonlinear problems of dimensional constraints and design variables. Such systematic and holistic trajectory optimization problems need to be considered and satisfied simultaneously under a unified model, and should not be dealt with in a single segment. In conclusion, multi-stage synchronous planning is a key direction for future research. In recent years, many scholars at home and abroad have proposed their own optimization design methods for the optimization of the return segment of recoverable launch vehicles, but most of them have one of the following problems:

(1) The model is more idealized and does not conform to the actual flight conditions of the entire return trajectory of the first sub-stage of the launch vehicle;

(2) The research mainly focuses on the single flight segment of the vertical landing segment;

(3) The flight time of the first sub-stage is short, the design variables and constraint variables are considered less, and the adjustment ability is limited, which cannot meet the needs of returning to the optimal trajectory of multiple flight segments;

(4) For different return tasks or different flight segments, it is necessary to reconstruct the optimization problem combined with the relevant algorithms, the derivation process is more complicated, and the calculation is more complicated;

(5) The self-adaptive optimization of the return trajectory and the automatic allocation of the algorithm cannot be realized.

In high-dimensional parameter space (parameter space composed of complex constraints and design variables) and high dynamics, traditional optimization methods can no longer meet the requirements in terms of accuracy, efficiency, optimality, adaptability, etc., and it is easy to appear non-convergence. Therefore, it is necessary to develop a more accurate and more applicable adaptive trajectory optimization method.

SUMMARY OF THE INVENTION

The purpose of the present invention is to provide a method based on cluster analysis for the return section of a recoverable rocket to solve the problem that the traditional optimization method cannot satisfy the optimization problem of the return section of the recoverable rocket from the perspectives of accuracy, efficiency, optimality, adaptability, etc. The adaptive optimization method can effectively deal with complex landing constraints, has high accuracy, convergence and computational efficiency, reduces the dimension and complexity of the parameter space of the optimization problem in the return segment, and has good adaptability and scalability.

For solving the above-mentioned technical problems, the technical scheme adopted in the present invention is:

An adaptive optimization method for the return segment of a recoverable rocket based on cluster analysis, comprising:

Step 1, establish a return trajectory optimization problem model of a sub-level return segment; the optimization problem model includes the flight segment division information, design variable information and constraint variable information of the return segment; wherein the constraint variable information includes design variable constraint information, process Constraint variable information and end constraint variable information;

It also features the following steps:

Step 2, solve the state transition matrix of the terminal constraint variables relative to the design variables;

Step 3: Use a clustering algorithm to perform sensitivity division on the data elements in the state transition matrix, and identify and divide data elements with higher sensitivity from the state transition matrix according to design requirements to form a state transition after preliminary dimensionality reduction. matrix M;

Step 4, according to the set linearity criterion, filter linear constraint variables, linear design variables and corresponding linear state transition matrices from M, and filter nonlinear constraint variables and nonlinear design variables from M;

Step 5: For linear constraint variables, linear design variables and corresponding linear state transition matrices, use a linear search algorithm to obtain optimal design variable values; for nonlinear constraint variables and nonlinear design variables, use nonlinear The direct optimization method in the optimization method obtains the optimal design variable values.

As a preferred way, the step 1 includes the following steps:

Step 11: Divide the return segment into four flight segments: the attitude adjustment segment, the power deceleration segment, the atmospheric reentry segment and the vertical landing segment;

Step 12, establishing a three-degree-of-freedom dynamic equation of a sub-level return trajectory;

Step 13: Determine the objective function, design variable information and constraint variable information of the first-level trajectory optimization.

As a preferred way, the three-degree-of-freedom dynamic equation established in the step 12 is:

分别为x、y、z、V、θ、σ、m的变化率；R为地球半径；r为地心距；P为一子级上发动机的额定推力；ε_n为变推力因子且0≤ε_n≤1，n＝0时ε_n为动力减速段的变推力因子，n＝1时ε_n为垂直着陆段的变推力因子；I_sp为发动机比冲；g₀为海平面重力加速度；g_r＝-μ/r²，μ为地球引力常数；α为一子级的攻角；β为一子级的侧滑角；X_q、Y_q、Z_q分别为一子级在速度坐标系x向、y向、z向下的气动力分量。Among them, x, y, z are the coordinates of a sub-stage in the x, y, and z directions of the emission coordinate system; V is the velocity value of a sub-stage; θ is the velocity inclination of a sub-stage; σ is a sub-stage yaw angle of ; m is the current quality of a sub-class;

are the rate of change of x, y, z, V, θ, σ, m respectively; R is the radius of the earth; r is the distance from the center of the earth; P is the rated thrust of the engine on a sub-stage; ε _n is the variable thrust factor and 0≤ ε _n ≤ 1, when n=0, ε _n is the variable thrust factor of the power deceleration stage, and when n=1, ε _n is the variable thrust factor of the vertical landing segment; I _sp is the engine specific impulse; g ₀ is the sea-level gravitational acceleration; g _r =-μ/r ² , μ is the gravitational constant of the earth; α is the angle of attack of a child; β is the sideslip angle of a child; X _q , Y _q , Z _q are the velocity coordinates of a child respectively It is the aerodynamic component in the x-direction, y-direction, and z-downward.

As a preferred way, the objective function determined in the step 13 is represented by the following formula: minJ=-m(t _f ), wherein J is the objective function, t _f is the landing time, and m(t _f ) is a sub the mass of the stage at the moment of landing.

As a preferred way, the design variable information determined in the step 13 is:

X ^s =[t ₁ t ₂ t ₃ t ₄ α ₁ α ₂ β ₁ β ₂ ε ₁ ε ₂ ] ^T ,

Among them, X ^s is a matrix composed of design variables, t ₁ is the taxiing time of a sub-stage in the attitude adjustment section, t ₂ is the flight time of a sub-stage in the power deceleration section, and t ₃ is a sub-stage in the atmospheric re-entry section. t ₄ is the flight time of a sub-stage in the vertical landing section, α ₁ is the attack angle of a sub-stage in the power deceleration section, α ₂ is the attack angle of a sub-stage in the vertical landing section, and β ₁ is a The sideslip angle of the child in the power deceleration section, β ₂ is the sideslip angle of a child in the vertical landing section, ε ₁ is the engine thrust factor of a child in the power deceleration section, and ε ₂ is the vertical landing of a child segment engine thrust factor.

As a preferred manner, the design variable constraint information includes attack angle constraint information, sideslip angle constraint information, and engine thrust constraint information.

As a preferred manner, the process constraint variable information includes one or more of dynamic pressure limitation conditions, overload limitation conditions, engine shutdown and ignition time interval conditions of two connected power stages, and terminal landing mass limitation conditions of a sub-stage.

As a preferred manner, the terminal constraint variable information includes: within the set accuracy range, the position of the first sub-stage when landing reaches the specified position, the speed of the first sub-stage landing reaches the specified speed, and the first sub-stage landing time reaches the specified speed, and the first sub-stage landing The posture is vertical.

As a preferred manner, in the step 3, K-means clustering algorithm is used to divide the sensitivity of the data elements in the state transition matrix.

As a preferred manner, when the K-means clustering algorithm is executed, the cluster center is initialized according to the data experience of the first-level return task model.

Compared with the prior art, the present invention can effectively deal with the complex landing constraint problem, has high precision, convergence and calculation efficiency, reduces the dimension and complexity of the parameter space of the optimization problem in the return segment, and has good self-adaptation. sex and expandability.

Description of drawings

Figure 1 is a flight sectional view of the return segment.

Figure 2 is a flow chart of the present invention.

Figure 3 is a graph of the optimization iterative process.

Figure 4 is a cross-sectional comparison diagram of the main parameters.

Detailed ways

The overall technical idea of the method of the present invention is as follows:

First, a problem model with multiple flight segments, multiple constraints, and multiple variables is established. Then, an adaptive optimization algorithm is proposed through the analysis of parameter characteristics: the high-dimensional parameter space state transition matrix of the terminal constraint variables relative to the design variables is solved, which describes the sensitivity of the terminal constraint variables relative to the design variables perturbation; The class algorithm divides the sensitivity of the data elements in the state transition matrix, and filters out linear constraint variables, linear design variables and linear state transition matrix according to the set linearity criterion, and filters out nonlinear constraint variables and design variables; Different optimization algorithms are automatically assigned according to the linear and nonlinear characteristics of different parameters.

The present invention is described in detail below in conjunction with the accompanying drawings:

1 Return to the establishment of the multi-flight segment model

The division model of one-stage return multi-flight segment, the dynamic model of return trajectory and the optimization problem model under multi-variable and multi-constraint are established.

1.1 One-child return multi-flight segment division

The trajectory optimization problem studied in the present invention is only for the recovery process of the first sub-stage of the launch vehicle, so the ascending section is not considered. A sub-level return recovery includes two ways of returning to the original field and not returning to the original field. The present invention studies the recovery method that does not return to the original field. A typical sub-stage that does not return to the original field usually undergoes multiple flight segments such as the attitude adjustment segment, the high-altitude re-entry segment, the power deceleration segment, the atmospheric re-entry segment, and the vertical landing segment during the return process after completing the ascent segment flight. For different flight missions, the return segment division is also different. Considering that the attitude adjustment section and the taxiing section are basically outside the dense atmosphere, the aerodynamic influence can be ignored, and the attitude adjustment section is also unpowered taxiing and relies on the reverse thrust control system for attitude adjustment, so these two flight sections are only affected by gravity. . The return segment division model of the present invention is simplified, and the return segment is divided into four flight segments: an attitude adjustment segment, a power deceleration segment, an atmospheric reentry segment and a vertical landing segment. The flight profile of the return segment is shown in Figure 1.

1.2 Establishment of the dynamic model

When studying the movement of the center of mass in the return process of the first sub-stage of the launch vehicle, without considering the control and adjustment process of its attitude and the rotation of the earth, the three-degree-of-freedom dynamic equation of the return trajectory of the first sub-stage is deduced, which is expressed by the following differential equations :

分别为x、y、z、V、θ、σ、m的变化率；R为地球半径；r为地心距；P为一子级上发动机的额定推力；ε_n为返回段第n次有动力飞行段(这里的有动力指的是发动机提供的推力，动力减速段是第一次有动力的飞行段，垂直着陆段是第二次有动力的飞行段，另外两个飞行段无动力)且0≤ε_n≤1，n＝0时ε_n为动力减速段的变推力因子，n＝1时ε_n为垂直着陆段的变推力因子；I_sp为发动机比冲；g₀为海平面重力加速度；g_r＝-μ/r²，μ为地球引力常数；α为一子级的攻角；β为一子级的侧滑角；X_q、Y_q、Z_q分别为一子级在速度坐标系x向、y向、z向下的气动力分量。Among them, x, y, z are the coordinates of a sub-stage in the x, y, and z directions of the emission coordinate system; V is the velocity value of a sub-stage; θ is the velocity inclination of a sub-stage; σ is a sub-stage yaw angle of ; m is the current quality of a sub-class;

are the rate of change of x, y, z, V, θ, σ, m respectively; R is the radius of the earth; _r is the distance from the center of the earth; P is the rated thrust of the engine on a sub-stage; Powered flight segment (powered here refers to the thrust provided by the engine, the powered deceleration segment is the first powered flight segment, the vertical landing segment is the second powered flight segment, and the other two flight segments are unpowered) And 0≤εn ≤1, when _n =0, εn is the variable thrust factor of the power deceleration segment, and when _n =1, εn _is the variable thrust factor of the vertical landing segment; _Isp is the engine specific impulse; g ₀ is the sea level Gravitational acceleration; g _r =-μ/r ² , μ is the gravitational constant of the earth; α is the angle of attack of a sub-level; β is the sideslip angle of a sub-level; X _q , Y _q , Z _q are respectively a sub-level Aerodynamic components in the x, y, and z directions of the velocity coordinate system.

X _q , Y _q , Z _q are defined as follows:

ρ₀为海平面大气密度，h₀为参考高度(为常数)；C_x、C_y、C_z分别为阻力、升力和侧力系数，且C_z＝-C_y，定义如下：In the formula, S _M is the characteristic reference area of the first sub-stage of the launch vehicle (constant); ρ is the local atmospheric density of the earth,

ρ ₀ is the air density at sea level, h ₀ is the reference height (which is a constant); C _x , C _y , and C _z are drag, lift and lateral force coefficients, respectively, and C _z = -C _y , defined as follows:

为常数。In the formula, C _x0 ,

is a constant.

1.3 Optimization problem description

The trajectory optimization of the first-level return segment refers to the goal of minimum fuel consumption, adjusting various design variables to make the aircraft meet the flight process constraints and land softly at the set landing point.

1.3.1 Objective function

In order to make the carrier rocket carry the best carrying capacity and have enough fuel for the guidance and correction maneuver, the performance index of the whole process of the first-stage return segment is the least fuel consumption, that is, the maximum landing mass at the end, so the objective function of trajectory optimization is:

minJ=-m(t _f ) (4)

Among them, J is the objective function, t _f is the landing moment, and m(t _f ) is the mass of a sub-stage at the landing moment.

1.3.2 Design variables

Design variables for each flight segment:

Attitude adjustment section: no power taxiing, relying on the reverse thrust control system to adjust the angle of attack to 0° to meet the ignition attitude requirements of the engine, the design variable is: taxiing time t ₁ .

Power deceleration section: turn on the engine to decelerate, reduce the re-entry velocity, and control the angle of attack and sideslip angle at the same time. The main design variables are: flight time t ₂ , engine thrust factor ε ₁ , angle of attack α ₁ , and sideslip angle β ₁ .

Atmospheric re-entry stage: This stage is designed to fly according to the 0° angle of attack, completely using aerodynamic deceleration, and the design variable is: flight time t ₃ .

Vertical landing stage: the engine is re-ignited to start the engine to decelerate, and the main design variables are: flight time t ₄ , engine thrust factor ε ₂ , angle of attack α ₂ , and sideslip angle β ₂ .

The flight time of each flight segment and the thrust coefficient of the power deceleration segment and the vertical landing segment can be obtained by direct optimization. The angle of attack and sideslip angle of the power deceleration section and the vertical landing section are optimized by designing discrete points, and between the discrete points, the angle of attack and sideslip angle can be obtained by a linear interpolation method. The program expressions of the attack angle and sideslip angle of the power deceleration section are:

The program expressions of the angle of attack and sideslip angle of the vertical landing segment are:

According to the designed flight procedure, the following variables are selected as design variables:

X ^s =[t ₁ t ₂ t ₃ t ₄ α ₁ α ₂ β ₁ β ₂ ε ₁ ε ₂ ] ^T (7)

1.3.3 Constraints

A child will be restricted by various complex constraints during the return process, mainly including design variable constraints, process constraints and precision landing constraints.

(1) Control constraints

Control constraints mainly include angle of attack, sideslip angle and engine thrust constraints.

分别为各个飞行段的攻角、侧滑角和推力的最小和最大值。

是时间变量,

分别为各个飞行段的开始时间和终端时间。

表示在时间区间

具有连续的一阶导数。In the formula, p=1, 2, 3, and 4 represent the attitude adjustment segment, the power descent segment, the atmospheric reentry segment and the vertical landing segment, respectively.

are the minimum and maximum values of the angle of attack, sideslip angle and thrust of each flight segment, respectively.

is the time variable,

are the start time and terminal time of each flight segment, respectively.

expressed in the time interval

has a continuous first derivative.

(2) Process constraints

The return process needs to consider the limitations of dynamic pressure and overload.

是速度系下的推力分量；n_max为箭体能够承受的最大过载；q_max为动压上限。In the formula,

is the thrust component under the velocity system; n _max is the maximum overload that the arrow body can bear; q _max is the upper limit of dynamic pressure.

At the same time, considering the working performance of the engine, there must be a certain interval between the ignition times of the two adjacent power sections when the engine is turned off and turned on.

Δt≥Δt _min (10)

In the formula, Δt _min is the lower limit of the time interval between engine shutdown and ignition of two adjacent power stages.

In addition, a sub-stage has undergone multiple flight stages before landing in the vertical landing segment, and the propellant remaining is limited, and it is necessary to ensure that the terminal landing mass is greater than the structural mass

m(t _f ) ^(p=4) ≥m _dry (11)

In the formula, m _dry is the mass of a sub-level structure, that is, the lower limit of the remaining mass during the flight.

The connection conditions between each flight segment are

(3) End precise landing constraints

The landing at the end of the vertical landing segment requires that the position and speed of the first stage must reach the specified target and accuracy, and its attitude must be vertical. In the research of the present invention, geodetic longitude, geodetic latitude and altitude are used as end position constraints, relative velocity is used as velocity constraints, and ballistic inclination is used as attitude constraints.

In the formula, λ ^f , B ^f , H ^f , V ^f , θ ^f are the geodetic longitude, geodetic latitude, altitude, relative velocity, and ballistic inclination of the landing point calculated by ballistic integration; λ _f , B _f , H _f , V _f , θ _f are the geodetic longitude, geodetic latitude, altitude, relative velocity, and ballistic inclination of the target landing site, respectively; ζ _λ , ζ _B , ζ _H , ζ _V , ζ _θ are the corresponding terminal landing accuracy, respectively.

For the convenience of expression, in the following research, the process inequality constraint is denoted as F, and the end precise landing constraint is denoted as G ^f =[λ ^f -λ _f ,B ^f -B _f ,H ^f -H _f ,V ^f -V _f ,θ ^f -θ _f ] ^T .

2 Algorithm design

Considering all kinds of complex constraints and design variables, an adaptive optimization algorithm strategy, which is an algorithm allocation strategy based on cluster analysis, is proposed. The high-dimensional parameter space state transition matrix obtained by the finite difference approximation is further processed, and the high-sensitivity data elements in the matrix and their corresponding design variables are automatically identified and divided by the cluster analysis method. The design variables are divided into linear parameters and nonlinear parameters, and finally different optimization algorithms are assigned to the linear parameters and nonlinear parameters. The algorithm is not only simple in logic, but also reduces the dimension of the high-dimensional parameter space optimization problem and realizes the adaptive allocation of the optimization algorithm. See the following sections for details.

2.1 High-dimensional parameter space state transition matrix solution

变化的主要设计变量

从而通过调整

使得

瞄准预定目标。For optimization problems in high-dimensional parameter spaces, it is necessary to analyze the main design variables that affect the end-constraint variables. Each element in the state transition matrix represents the sensitivity of the constraint variables to the perturbation of the design variables. Therefore, by identifying the parameters of the elements in the state transition matrix, the influencing constraint variables can be obtained.

Changed key design variables

so that by adjusting

make

Aim for the intended target.

等式约束

(本发明优化问题模型中n＝10,m＝5)，它们之间存在非线性函数关系To make the problem general, note the design variables

equality constraints

(n=10, m=5 in the optimization problem model of the present invention), there is a nonlinear functional relationship between them

G ^f =f(X ^s ,t) (14)

The state transition matrix between the end constraints and the design variables is

表示X^s中对应的设计变量；

为末端约束G^f对设计变量X^s的状态转移矩阵，描述了约束变量对设计变量摄动的敏感度。目前求解状态转移矩阵有两种常用的方法，一种是计算出状态转移矩阵的显示表达式，在轨迹积分中同步求解；另一种是采用数值差分方法近似求解。一子级的返回段动力学模型比较复杂，末端约束和设计变量之间没有显示的函数关系，因此得到不到状态关系矩阵的解析形式。本发明在计算中采用有限差分近似公式求解

中的各个元素：In the formula,

represents the corresponding design variable in X ^s ;

For the state transition matrix of the end constraint G ^f to the design variable X ^s , the sensitivity of the constraint variable to the perturbation of the design variable is described. At present, there are two commonly used methods for solving the state transition matrix, one is to calculate the display expression of the state transition matrix and solve it simultaneously in the trajectory integration; the other is to use the numerical difference method to approximate the solution. The dynamic model of the return segment of the first sub-stage is relatively complex, and there is no explicit functional relationship between the end constraints and the design variables, so the analytical form of the state relationship matrix cannot be obtained. The present invention adopts finite difference approximation formula to solve

Elements in:

为

的偏差，偏差过大会影响收敛性和计算精度，偏差过小会引起计算误差，根据数值仿真经验，通常取0.1％-1％设计变量区间的小偏差(本发明取1％)。由此，根据数值仿真计算可求出状态转移矩阵。When using the formula (16) to calculate,

for

If the deviation is too large, it will affect the convergence and calculation accuracy, and if the deviation is too small, it will cause calculation errors. Thus, the state transition matrix can be obtained according to the numerical simulation calculation.

2.2 Cluster Analysis of Data Elements of State Transition Matrix

For different return tasks, changes in flight segment tasks or design variables will lead to changes in the state transition matrix, so it is necessary to use a method that can automatically identify matrix elements to avoid manual participation. Cluster analysis is an important data analysis method in statistics and has been widely used in pattern recognition, machine learning, data analysis and image processing. , the whole division is an unsupervised process, and finally the data similarity within the same cluster is as high as possible and the data between clusters is as independent as possible. The invention adopts the classical K-means clustering algorithm to perform clustering processing on the data elements in the state transition matrix. The following briefly introduces the clustering process of matrix data elements by the K-means clustering algorithm in combination with the state transition matrix in the high-dimensional parameter space.

中的元素进行min-max标准化，对原始元素数据进行线性变换，使结果落于[0,1]区间，记标准化后的矩阵为A，则A中对应元素为(1) First convert the state transition matrix

The elements in the min-max are normalized, and the original element data is linearly transformed so that the result falls in the [0,1] interval, and the standardized matrix is A, then the corresponding element in A is

中，记标准化矩阵A＝[A₁,...,A_i,...,A_m]^T，可知每个A_i包含n个数据元素，即A_i＝(A_i1,...,A_ij,...,A_in)，需要将数据元素分为p类。选择p个初始聚类中心，本发明研究中，将数据元素分为两类，即高敏感度参数和低敏感度参数，设初始迭代次数为b＝0，初始化需选择两个聚类中心

(2) In the parameter space

, denote the standardized matrix A=[A ₁ ,...,A _i ,...,A _m ] ^T , it can be known that each A _i contains n data elements, that is, A _i =(A _i1 ,..., A _ij ,...,A _in ), the data elements need to be divided into p classes. Select p initial cluster centers. In the research of the present invention, the data elements are divided into two categories, namely high-sensitivity parameters and low-sensitivity parameters, and the initial iteration number is set to b=0, and two cluster centers need to be selected for initialization

在A_i中把与聚类中心

距离近的数据元素归到相应的簇中，即若

就把A_ij划分到簇

中。(3) For any data element in A _i , find its Euclidean distance to p cluster centers:

Put and cluster center in A _i

Data elements with close distances are grouped into the corresponding clusters, that is, if

Divide _Aij into clusters

middle.

其中A_ij为类

中的数据元素，N_p为类

中数据元素的个数。(4) Calculate the new cluster center according to the data elements in each cluster

where _Aij is the class

The data elements in , where N _p is the class

The number of data elements in .

时表示算法收敛，即聚类中心在迭代更新后数值保持不变，ε为迭代停止阈值。否则令b＝b+1，转到步骤(3)，继续迭代。(5) When

When it means that the algorithm converges, that is, the value of the cluster center remains unchanged after iterative update, and ε is the iteration stop threshold. Otherwise, let b=b+1, go to step (3), and continue to iterate.

The K-means clustering algorithm is simple, efficient, and easy to implement. The simulation shows that the convergence value can be iterated for 5 to 10 times. However, studies have shown that the selection of cluster centers will affect the convergence and efficiency of clustering, and the random selection of initial cluster centers will also bring greater uncertainty to the algorithm. Therefore, when the state transition matrix is clustered, the reasonable selection of the cluster center can improve the accuracy of the parameter sensitivity clustering results. Based on the guidance of prior information, the invention initializes the clustering center according to the data experience of a large number of first-level return task models, improves the clustering effect in the early stage of iteration, guides the accurate clustering of data elements, and at the same time reasonably imposes constraints in the clustering process to improve the clustering performance.

范围内影响约束变量的主要设计变量，记为

实现了问题模型的初步降维。记初步降维后的状态转移矩阵为M，则So far, the high-sensitivity data elements in the state transition matrix are automatically screened out through cluster analysis, so as to obtain the results in the local interval.

The main design variables in the range that affect the constraint variables, denoted as

The initial dimensionality reduction of the problem model is achieved. Denote the state transition matrix after initial dimensionality reduction as M, then

表示

中对应的设计变量；p_n为聚类后主要设计变量的个数，且p_n＜n，

In the formula,

express

The corresponding design variables in ; p _n is the number of main design variables after clustering, and p _n <n,

作为输入，利用打靶法计算设计变量在设计区间内变化时输出约束变量的变化值，对应的变化矩阵定义为差分矩阵，记为N。用矩阵N中的元素除以M对应位置的元素，得到的矩阵记为C，且C_ik＝N_ik/M_ik，其中每个元素的含义是设计变量在设计区间和局部区间内变化时，末端约束变量变化的倍数。本发明研究中偏差

取1％设计变量设计区间，即局部偏差到设计区间扩大了100倍。定义线性度指标l(本发明研究中，线性度指标取l＝[0.9,1.1]×100)，若由聚类分析得到的高敏感度数据元素所在的位置在C中对应的主要元素的值在线性指标范围之内，说明该主要数据元素对应的设计变量和约束变量之间存在全局的线性关系，由此可筛选出设计变量和约束变量中存在的线性参数及对应的线性状态转移矩阵(记为Q)，从而实现了线性参数与非线性参数的解耦，同时也对优化问题进行再次降维。Next, the design variables will be

As the input, the target method is used to calculate the change value of the output constraint variable when the design variable changes within the design interval, and the corresponding change matrix is defined as the difference matrix, denoted as N. Divide the element in matrix N by the element at the corresponding position of M, and the obtained matrix is denoted as C, and C _ik =N _ik /M _ik , where the meaning of each element is that when the design variable changes in the design interval and the local interval, The fold of change in the end constraint variable. Bias in the present study

Taking 1% of the design variable design interval, that is, the local deviation to the design interval is expanded by 100 times. Define the linearity index l (in the study of the present invention, the linearity index takes l=[0.9,1.1]×100), if the position of the high-sensitivity data element obtained by the cluster analysis is the value of the corresponding main element in C Within the range of the linear index, it indicates that there is a global linear relationship between the design variables and the constraint variables corresponding to the main data element, so that the linear parameters existing in the design variables and the constraint variables and the corresponding linear state transition matrix ( Denoted as Q), which realizes the decoupling of linear parameters and nonlinear parameters, and also reduces the dimension of the optimization problem again.

2.3 Optimization Problem Reconstruction and Adaptive Allocation of Optimization Algorithms

其中，

为线性参数设计变量，

为降维后的非线性参数优化设计变量。记等式约束

其中，

为线性参数约束变量，

为降维后的非线性参数末端约束变量。对于线性参数，有线性搜索算法，如牛顿迭代、拟牛顿迭代、微分修正方法等。值得注意的是，线性搜索算法在优化中可以增强解的鲁棒性。本发明采用微分修正法进行迭代求解，它本质上是一种牛顿迭代逼近期望值算法，即迭代的打靶法，通过不断调整设计变量使得约束变量最终收敛于期望值，具有足够的适应性和鲁棒性。通过敏感度和线性度的分析可以得到

相对于X_1s的线性状态转移矩阵Q，则Through the above cluster analysis, the parameters are automatically classified for sensitivity and linearity, which realizes the parameter decoupling between the end constraint variables and the design variables, and reduces the complexity and dimension of the optimization problem. Then according to the linear and nonlinear characteristics of the parameters, the corresponding optimization algorithm can be assigned to the adaptation. design variables

in,

Design variables for linear parameters,

Optimizing design variables for nonlinear parameters after dimensionality reduction. note equality constraints

in,

Constrain the variables for linear parameters,

is the end constraint variable of the nonlinear parameter after dimension reduction. For linear parameters, there are linear search algorithms, such as Newton iteration, quasi-Newton iteration, and differential correction methods. It is worth noting that the linear search algorithm can enhance the robustness of the solution in optimization. The invention adopts the differential correction method for iterative solution, which is essentially a Newton iterative approximation to the expected value algorithm, that is, the iterative shooting method. By continuously adjusting the design variables, the constraint variables finally converge to the expected value, which has sufficient adaptability and robustness. . Through the analysis of sensitivity and linearity, we can get

With respect to the linear state transition matrix Q of X _1s , then

和迭代设计变量

自由度不等的情况，并增加迭代因子μ增加迭代的精度和收敛性，则For generality, a least squares strategy is used to solve the constraint variables

and iterative design variables

In the case of unequal degrees of freedom, and increasing the iteration factor μ increases the iteration accuracy and convergence, then

In the formula, μ is the iteration factor, and 0<μ≤1. It is worth noting that, in theory, parameters with a strict linear relationship can converge to the expected value in one iteration. The linear parameters defined in the present invention are measured by the linearity index l, which is not strictly linear. Simulation shows that the constraint variables converge to the expected value. Accuracy requires 5-10 iterations.

可采用直接优化方法，如序列二次规划法(SQP)、内点法、伪谱法等方法进行优化，从而使得非线性末端约束变量

和不等式约束

收敛于要求的精度，目标函数

可以通过对

微分迭代修正和对

优化搜索得到最优指标。注意到使用微分修正算法和直接优化方法需提供初值解，本发明通过遗传算法(GA)对原优化问题优化得到初值解。For nonlinear design variables after dimensionality reduction

Direct optimization methods, such as sequential quadratic programming (SQP), interior point method, pseudospectral method, etc., can be used for optimization, so that the nonlinear end-constraint variable

and inequality constraints

converges to the required accuracy, the objective function

through the

Differential iterative correction and pair

Optimize the search to get the best index. It is noted that the use of the differential correction algorithm and the direct optimization method needs to provide the initial value solution, and the present invention optimizes the original optimization problem through the Genetic Algorithm (GA) to obtain the initial value solution.

2.4 Algorithm Summary

The overall flow of the present invention is shown in Figure 2, which mainly includes three aspects: state transition matrix solution, cluster analysis, problem reconstruction and algorithm allocation; the three are progressive, state transition matrix is the basis of the algorithm, cluster analysis It is the core link, and problem reconstruction and algorithm allocation are the key parts.

The adaptive optimization algorithm starts from the solution of the state transition matrix, identifies and extracts high-sensitivity data elements in the matrix through the cluster analysis algorithm, and then further performs linear analysis combined with the difference matrix to accurately identify the linear and nonlinear parameters and analyze the linearity. State transition matrix, and finally assign different optimization algorithms according to different parameter characteristics. It is worth noting that once the division of flight segments, the number of design variables and the interval of design variables are determined, the linear state transition matrix is a constant value matrix, which does not need to be updated in the process of differential correction calculation, which greatly improves the operation efficiency. In a word, the biggest advantage of this algorithm is that it does not require manual participation in the whole process, and the algorithm is simple and efficient, and can accurately realize automatic identification and classification of parameters and automatic allocation of optimization algorithms.

3 Simulation verification

然后将状态转移矩阵min-max标准化得到标准化矩阵A，如表1所示。对A中数据元素进行聚类分析，得到高敏感度元素，如A中的粗斜元素所示，对应于状态转移矩阵中的元素如表1中

粗斜元素所示。通过约束变量对设计变量的敏感度分析得到影响约束变量的主要设计变量，即粗斜元素所对应的的设计变量。Calculate the state transition matrix according to the finite difference approximation formula

Then normalize the state transition matrix min-max to get the normalized matrix A, as shown in Table 1. Cluster analysis is performed on the data elements in A, and high-sensitivity elements are obtained, as shown by the thick oblique elements in A, corresponding to the elements in the state transition matrix as shown in Table 1

Bold italic elements are shown. Through the sensitivity analysis of constraint variables to design variables, the main design variables that affect the constraint variables, that is, the design variables corresponding to the rough slope elements, are obtained.

It can be found that within the range of local deviation, the main design variable that affects the change of the terminal landing longitude λ ^f is the sideslip angle β ₁ of the power descent section, and the main design variable that affects the change of the terminal landing latitude B ^f is the flight time t of the attitude adjustment section. _1. The main design variables that affect the change of the terminal landing height H ^f are the flight time t ₁ and t ₃ of the attitude adjustment section and the atmospheric re-entry section, and the main design variable affecting the change of the terminal landing speed V ^f is the flight time t of the vertical landing section ₄ and thrust factor ε ₂ , the main design variable that affects the change of terminal landing speed inclination angle θ ^f is the flight time t ₄ of the vertical landing segment, where t ₁ and t ₄ are coordinated control variables, and t ₁ is the main control of B and H at the same time quantity, t ₄ is the main controlling quantity of both V and θ.

Table 1 State transition matrix and its corresponding normalization matrix

和初步降维后的状态转移矩阵M，将

作为输入，利用打靶法计算

在设计变量全区间范围内变化时输出约束变量的变化值，得到全局差分矩阵N，矩阵M和N如表2所示。用N除以M中的对应元素得到矩阵C，可以发现，当设计变量范围变化100倍时，λ^f,B^f,H^f对应的粗斜元素的数值在线性度范围之内，而V^f,θ^f对应的粗斜元素不在线性度范围之内。因此，筛选得到线性设计变量t₁,t₃,β₁和线性约束变量λ^f,B^f,H^f以及线性状态转移矩阵Q。The above cluster analysis obtained the main design variables that affect the end constraint variables

and the state transition matrix M after initial dimensionality reduction, the

As input, the target method is used to calculate

When the design variables are changed in the whole interval, the change values of the constraint variables are output, and the global difference matrix N is obtained. The matrices M and N are shown in Table 2. Divide N by the corresponding element in M to get the matrix C. It can be found that when the range of design variables changes by 100 times, the values of the rough slope elements corresponding to λ ^f , B ^f , and H ^f are within the linearity range, while V ^f , θ ^f corresponds to the thick oblique element that is not within the linearity range. Therefore, the linear design variables t ₁ , t ₃ , β ₁ and the linear constraint variables λ ^f , B ^f , H ^f and the linear state transition matrix Q are obtained by screening.

Table 2 Preliminary dimensionality reduction state transition matrix M, difference matrix N, matrix C

The linear state transition matrix Q is expressed as

Therefore, the relationship between the correction amount of λ ^f -λ _f , B ^f -B _f , H ^f -H _f and the adjustment amount of t ₁ , t ₃ , β ₁ is as follows

Through parameter sensitivity and linearity analysis, the linear partial derivative matrix between t ₁ , t ₃ , β ₁ and λ ^f -λ _f , B ^f -B _f , H ^f -H _f is obtained, that is, the differential correction matrix, through Iteratively corrects t ₁ , t ₃ , β ₁ so that λ ^f , B ^f , H ^f converge to the expected values λ _f , B _f , H _f , and the design variables and constraint variables are reduced accordingly when the nonlinear optimization algorithm is used for optimization three dimensions. In optimization problem reconstruction,

The self-adaptive optimization algorithm (abbreviated as AOA) designed by the present invention is compared and verified with the SQP algorithm in the prior art. The end landing position accuracy of the two methods is shown in Table 3. It can be seen that the accuracy of the terminal landing position (including longitude, latitude and altitude) iteratively corrected by AOA is higher than that obtained by the SQP method. The terminal landing speed and speed inclination obtained by the two methods both meet the accuracy requirements, and the effects are equivalent. It shows that the AOA method has good convergence and accuracy.

Table 3 End-to-end precision landing accuracy comparison

It can be seen from Figure 3 that the optimization iterations of the SQP algorithm are 43 times, and the optimization iterations of the AOA algorithm are 16 times, which is only 37.2% of the SQP algorithm. At the same time, the terminal landing mass obtained by AOA is 22582kg, and the terminal landing mass obtained by SQP is 22136kg. Therefore, the AOA method is superior to the SQP algorithm in terms of computational efficiency and optimality.

Figures 4(a)-(h) give the profiles of the main parameters. Figure 4(a)-(b) shows the profile curves of overload and dynamic pressure. It can be seen that the results obtained by the two methods satisfy all constraints, and Fig. 4(c)-(g) gives the profile curves of the main parameters of the terminal landing. In the vertical landing segment, the longitude and latitude are almost unchanged, and the first stage reaches the area above the target landing point, and mainly adjusts the altitude, speed and speed inclination. Figure 4(g) shows the mass profile comparison of the two methods. The slope of the curve of the mass profile represents the mass-second consumption of the engine. During the attitude and atmospheric reentry phases, the mass remains unchanged due to engine shutdown. In the power descent section and the vertical landing section, the mass curves of the two methods are almost parallel, indicating that the thrust factors of the optimized engines are approximately equal. However, the engine startup time of the AOA method in the two power stages is 80.2s and 35.4s respectively, while the engine startup time of the SQP method in the two power stages is 87.6s and 31.4s, respectively. The fuel mass consumed by the AOA method is 12918kg, The fuel consumed by the SQP method was 13264 kg. The entire return process of the AOA method took 398.5s, and the entire return process of the SQP method took 396.2s. For such a free-time optimization problem, the AOA method ensures the optimality of engine ignition timing and fuel consumption.

To sum up, the present invention has the following advantages:

(1) A systematic and holistic trajectory optimization problem model for the entire return segment and multiple flight segments of the first sub-stage of the launch vehicle is established, which is considered and satisfied simultaneously under a unified model.

(2) On the basis of the high-dimensional parameter space state transition matrix, based on the K-means clustering algorithm, the accurate identification and positioning of the highly sensitive data elements in the state transition matrix is successfully realized, and the linearity and nonlinearity are filtered out by linearity analysis. parameters and linear partial derivative matrices; finally, different optimization algorithms are assigned to the linear and nonlinear characteristics of different design variables and constraint variables, which have good accuracy, convergence, computational efficiency and optimality.

(3) In the process of cluster analysis, the main design variables that affect the change of constraint variables can be obtained, which helps to provide a reference for the adjustment and design of the nominal trajectory, and reduces the complexity and dimension of the optimization problem through sensitivity and linearity analysis. , which realizes the decoupling of linear parameters and nonlinear parameters.

(4) It realizes the automation of parameter identification and classification and the self-adaptation of algorithm allocation, which basically does not require manual participation, has good reliability, accuracy and scalability, and also has reference significance for other types of high-dimensional optimization problems.

The embodiments of the present invention have been described above in conjunction with the accompanying drawings, but the present invention is not limited to the above-mentioned specific embodiments, which are merely illustrative rather than limiting. Under the inspiration of the present invention, without departing from the scope of protection of the spirit of the present invention and the claims, many forms can be made, which all fall within the protection scope of the present invention.

Claims (10)

1. An adaptive optimization method for the return segment of a recoverable rocket based on cluster analysis, comprising: Step 1, establish a return trajectory optimization problem model of a sub-level return segment; the optimization problem model includes the flight segment division information, design variable information and constraint variable information of the return segment; wherein the constraint variable information includes design variable constraint information, process Constraint variable information and terminal constraint variable information; It is characterized in that, also comprises the following steps: Step 2, solve the state transition matrix of the terminal constraint variables relative to the design variables; Step 3: Use a clustering algorithm to perform sensitivity division on the data elements in the state transition matrix, and identify and divide data elements with higher sensitivity from the state transition matrix according to design requirements to form a state transition after preliminary dimensionality reduction. matrix M; Step 4, according to the set linearity criterion, filter linear constraint variables, linear design variables and corresponding linear state transition matrices from M, and filter nonlinear constraint variables and nonlinear design variables from M; Step 5: For linear constraint variables, linear design variables and corresponding linear state transition matrices, use a linear search algorithm to obtain optimal design variable values; for nonlinear constraint variables and nonlinear design variables, use nonlinear The direct optimization method in the optimization method obtains the optimal design variable values. 2. The self-adaptive optimization method for the return section of a recoverable rocket based on cluster analysis as claimed in claim 1, wherein the step 1 comprises the following steps: Step 11: Divide the return segment into four flight segments: the attitude adjustment segment, the power deceleration segment, the atmospheric reentry segment and the vertical landing segment; Step 12, establishing a three-degree-of-freedom dynamic equation of a sub-level return trajectory; Step 13: Determine the objective function, design variable information and constraint variable information of the first-level trajectory optimization. 3. The self-adaptive optimization method for the return section of a recoverable rocket based on cluster analysis as claimed in claim 2, wherein the three-degree-of-freedom dynamic equation established in the step 12 is:

Among them, x, y, z are the coordinates of a sub-stage in the x, y, and z directions of the emission coordinate system; V is the velocity value of a sub-stage; θ is the velocity inclination of a sub-stage; σ is a sub-stage yaw angle of ; m is the current quality of a sub-class;

are the rate of change of x, y, z, V, θ, σ, m respectively; R is the radius of the earth; r is the distance from the center of the earth; P is the rated thrust of the engine on a sub-stage; ε _n is the variable thrust factor and 0≤ ε _n ≤ 1, when n=0, ε _n is the variable thrust factor of the power deceleration stage, and when n=1, ε _n is the variable thrust factor of the vertical landing segment; I _sp is the engine specific impulse; g ₀ is the sea-level gravitational acceleration; g _r =-μ/r ² , μ is the gravitational constant of the earth; α is the angle of attack of a child; β is the sideslip angle of a child; X _q , Y _q , Z _q are the velocity coordinates of a child respectively It is the aerodynamic component in the x-direction, y-direction, and z-downward. 4. The self-adaptive optimization method for the return section of a recoverable rocket based on cluster analysis according to claim 2 or 3, wherein the objective function determined in the step 13 is represented by the following formula: minJ=-m(t _f ), where J is the objective function, t _f is the landing moment, and m(t _f ) is the mass of a sub-stage at the landing moment. 5. The self-adaptive optimization method for the return section of a recoverable rocket based on cluster analysis according to claim 2 or 3, wherein the design variable information determined in the step 13 is: X ^s =[t ₁ t ₂ t ₃ t ₄ α ₁ α ₂ β ₁ β ₂ ε ₁ ε ₂ ] ^T , Among them, X ^s is a matrix composed of design variables, t ₁ is the taxiing time of a sub-stage in the attitude adjustment section, t ₂ is the flight time of a sub-stage in the power deceleration section, and t ₃ is a sub-stage in the atmospheric re-entry section. t ₄ is the flight time of a sub-stage in the vertical landing section, α ₁ is the attack angle of a sub-stage in the power deceleration section, α ₂ is the attack angle of a sub-stage in the vertical landing section, and β ₁ is a The sideslip angle of the child in the power deceleration section, β ₂ is the sideslip angle of a child in the vertical landing section, ε ₁ is the engine thrust factor of a child in the power deceleration section, and ε ₂ is the vertical landing of a child segment engine thrust factor. 6 . The adaptive optimization method for the return segment of a recoverable rocket based on cluster analysis according to claim 2 , wherein the design variable constraint information comprises attack angle constraint information, sideslip angle constraint information, and engine thrust constraints. 7 . information. 7. The cluster analysis-based adaptive optimization method for the return section of a recoverable rocket according to claim 2, wherein the process constraint variable information includes dynamic pressure constraints, overload constraints, and engine closures of two connected power stages. One or more of the power-on ignition time interval conditions and the terminal landing mass limit conditions of a sub-stage. 8. The cluster analysis-based adaptive optimization method for the return segment of a recoverable rocket according to claim 2, wherein the terminal constraint variable information comprises: within a set accuracy range, the The position reaches the specified position, the speed of the first stage landing reaches the specified speed, and the attitude of the first stage landing is vertical. 9. The self-adaptive optimization method for the return segment of a recoverable rocket based on cluster analysis as claimed in claim 1, wherein in the step 3, K-means clustering algorithm is used to convert the data in the state transition matrix. Elements are sensitively divided. 10. The self-adaptive optimization method for the return segment of a recoverable rocket based on cluster analysis as claimed in claim 9, wherein when performing the K-means clustering algorithm, the clustering is initialized according to the data experience of the first-level return mission model center.

Images