A tight formation path planning method for aerial refueling

Abstract

The present invention discloses a close formation path planning method for aerial refueling, belonging to the technical field of aerial refueling path planning for unmanned aerial vehicles. Based on wind tunnel test data of the rendezvous phase and the formation-entering phase, the present invention obtains the key points of the receiving aircraft's path in the two phases of aerial refueling, defines the vortex negative impact area as an obstacle body, builds a three-dimensional obstacle model for path planning, and obtains a path planning map; based on the key points of the path and the set constraints, uses a key point heuristic algorithm to perform initial planning in the map; establishes a mathematical model of the fixed-wing unmanned aerial vehicle and analyzes its differential flatness characteristics; uses the differential flatness characteristics to establish a cost function for the path optimization process, and uses a flat output minimum mutation algorithm in combination with the constraints to optimize the path obtained by the initial planning, thereby obtaining the final aerial refueling path. While ensuring safety, the present invention takes into account the impact of vortices on the receiving aircraft and obtains a smooth and feasible path.

Description

Technical Field

[0001] The present invention belongs to the technical field of aerial refueling path planning for unmanned aerial vehicles (UAVs), and in particular relates to a close formation path planning method for aerial refueling. Background Art

[0002] Aerial refueling for drones (UAVs) refers to the autonomous refueling mission between a tanker and a receiving aircraft, without pilot intervention. This involves autonomously completing steps including rendezvous, pre-docking, docking, refueling, and formation recovery, thus automating the aerial refueling function for unmanned aircraft. Aerial refueling can significantly improve fighter aircraft's endurance and payload capacity, increasing the range and combat radius of combat aircraft and extending their flight time. Path planning for aerial refueling involves planning a path from a starting point to a destination while satisfying certain constraints, optimizing specified performance indicators. Constraints primarily include spatial coordination constraints, temporal coordination constraints, and UAV constraints. Performance indicators include path length, path smoothness, path safety, and mission completion time. When a tanker and receiving aircraft fly in close formation during aerial refueling operations, the tanker generates wingtip vortices, which significantly impact the receiving aircraft. Furthermore, the receiving aircraft must maintain a precise trajectory relative to the tanker to complete the refueling operation. Therefore, real-time and safe trajectory planning during aerial refueling is crucial.

[0003] Currently, many methods have been developed to address the path planning problem for drones. However, most traditional path planning methods and intelligent path optimization algorithms suffer from high computational complexity and are prone to falling into local optimal solutions. Furthermore, most fail to consider the influence of eddy currents and only plan paths in a two-dimensional environment without considering the influence of eddy currents. Furthermore, most studies only consider spatial constraints when planning paths, ignoring temporal constraints. This makes it impossible to achieve UAV refueling path planning that satisfies continuous temporal and spatial variations. Summary of the Invention

[0004] To address the above issues, the present invention provides a tight formation path planning method for aerial refueling, which can be used for aerial refueling path planning in a vortex environment. While ensuring safety, the present invention takes into account the impact of vortices on the receiving aircraft and obtains a smooth and feasible path.

[0005] Based on wind tunnel test data from the rendezvous and formation phases, this method calculates the impact of vortices on the receiving aircraft at different locations, deriving key points along the receiving aircraft's path during both phases of aerial refueling. Areas negatively impacted by vortices are defined as obstacles, and a three-dimensional obstacle model is constructed for path planning. This model then generates an aerial refueling path planning map. Based on the determined key points of the path, a keypoint heuristic algorithm is used to perform initial planning within the map, subject to pre-defined constraints. A mathematical model of a fixed-wing unmanned aerial vehicle (UAV) is established, and its differential flatness properties are analyzed, allowing the UAV's position and its derivatives to be used to describe other state variables and control variables of the UAV system. This differential flatness property is used to establish a cost function for the path optimization process. A flat output minimum catastrophe algorithm is then used, combined with the constraints, to optimize the initially planned path, resulting in the final aerial refueling path. This method, while ensuring safety, accounts for the impact of vortices on the receiving aircraft, resulting in a smooth, feasible path.

[0006] The technical solution adopted by the present invention is: a close formation path planning method for aerial refueling, the method comprising the following steps:

[0007] Step 1: Based on wind tunnel test data from the rendezvous and formation phases of aerial refueling, key points of the receiving aircraft's path during the rendezvous and formation phases are obtained according to the changes in the lift coefficient, drag coefficient, roll, pitch, and yaw moment coefficients of the receiving aircraft at different positions in the vortex.

[0008] A three-dimensional obstacle model for path planning is built based on the selected obstacle (e.g., the area with negative impact of vortex is defined as the obstacle), and a map for aerial refueling path planning is obtained;

[0009] Step 2: Determine the constraints of the two phases of aerial refueling. Based on the key points of the receiving aircraft's path in the rendezvous phase and the formation phase, use the Key-point Heuristic A* algorithm (KHA*) to perform initial path planning in the aerial refueling path planning map to obtain the nominal planned path.

[0010] Step 3: Establish a fixed-wing UAV model and analyze the differential flatness characteristics of the fixed-wing UAV to describe the state and control variables of the UAV system using the position of the UAV and its derivatives.

[0011] Step 4: Based on the differential flatness characteristics of the fixed-wing UAV, the cost function of the path optimization process is established. Based on the constraints of the two stages of aerial refueling, the Flat-output MinimumSnap algorithm (FMS) is used to optimize the nominal planning path. Optimize and obtain the final aerial refueling path.

[0012] Furthermore, in step 1, the key points of the receiving aircraft's path in the rendezvous phase and the formation phase are determined based on the minimum value of the cost function W;

[0013] Among them, the expression of the cost function W is:

[0014]

[0015] Where Δcmx represents the rolling moment coefficient, Δcmy represents the pitching moment coefficient, Δcmz represents the yaw moment coefficient, K represents the lift-to-drag ratio (lift / drag), and λ i ,i∈[1,2,…,4] represents the weight of the corresponding rolling moment coefficient, pitching moment coefficient, yaw moment coefficient or lift-to-drag ratio.

[0016] Preferably, λ1=λ2=λ3=1, λ4=0.001.

[0017] Furthermore, in step 1, the area where the negative impact of the eddy current exceeds a set threshold is defined as an obstacle.

[0018] Furthermore, in step 2, the constraints include but are not limited to: path length constraint, yaw angle constraint, pitch angle constraint, relative distance constraint between the two aircraft, etc.

[0019] Furthermore, step 2 includes:

[0020] The aerial refueling path planning map is divided into multiple three-dimensional elements according to a certain step size, and obstacle blocks and non-obstacle blocks are marked. Each three-dimensional element represents a path node (referred to as node);

[0021] Then, the KHA* algorithm is used to trigger from the set starting point, gradually passing through all the key points. Each node found between two key points must satisfy the minimum movement cost function, and each path node on the path is found in turn, and then connected with the key points to form the nominal planning path.

[0022] The movement cost function includes the movement cost from the starting point to the path node and the estimated movement cost from the path node to the end point.

[0023] Furthermore, step 3 includes:

[0024] Define coordinate system G i ,G b and G s They represent the inertial reference frame, the body reference frame and the stable reference frame respectively;

[0025] Among them, the inertial reference system G iThe origin is located at a point on the ground, x i The axis lies in the ground plane and points in a specific direction, y i Axis and x i The axes are in the same plane and perpendicular to x i Axis, z i The axis points vertically downward, toward the center of the Earth;

[0026] Body reference frame G b The origin is located at the center of mass of the drone, x b The axis is in the symmetry plane of the UAV, parallel to the fuselage axis of the UAV or the average aerodynamic chord of the wing, and points to the direction of the UAV's nose. b The axis is perpendicular to the symmetry plane of the drone and points to the right side of the fuselage. b with x b The axis is vertical and points downwards from the drone;

[0027] Stable reference frame G s The origin is located at the center of mass of the drone, x s The axis coincides with the projection of the drone's airspeed in the symmetry plane, y s Axis and body coordinate axis y b Overlap, z s In the symmetry plane of the aircraft and x s Vertical and pointing down the fuselage;

[0028] Define the regular basis vectors in the three-dimensional coordinate system as e1, e2 and e3, and establish the fixed-wing UAV model as:

[0029]

[0030] Among them, O, v, m, ω are the position, velocity, mass and angular velocity of the receiving aircraft respectively, g is the acceleration of gravity, ω c is the specified angular rate, α is the angle of attack, and D is the velocity of the drone along the negative x b The drag force of the UAV along the x-axis, T is b The thrust generated by the axis, L is the thrust generated by the UAV along the negative z axis. b The lift generated by the shaft; Denotes the stable coordinate system G s To the inertial coordinate system G i The stability-inertia selection matrix, Represents the coordinate system G b The position O in the coordinate system G i In the expression The first derivative of is the coordinate system G b The velocity in the coordinate system G i In the expression, for The first derivative of In the coordinate system G i Measurement coordinate system G b The angular velocity in the coordinate system G b In the statement, Represents the coordinate system G b To coordinate system G i The body-inertia selection matrix, for The first derivative of .

[0031] Furthermore, the drag in the fixed-wing drone model is set to:

[0032]

[0033] Where ρ is the air density, V a is the airspeed of the drone, S is the plane area of ​​the wing, c is the chord length of the wing, and is the aerodynamic coefficient.

[0034] Furthermore, in step 3, define the airspeed of the drone when there is no wind And define the state x and control input u of the drone system as:

[0035]

[0036] The state x and control input u are described in terms of the flat output z and its high-order derivatives as:

[0037]

[0038] Furthermore, in step 4, the cost function of the path optimization process is established based on the differential flatness characteristics of the fixed-wing UAV, including:

[0039] According to the differential flatness characteristics of fixed-wing UAV, based on the current coordinate position x of the UAV in space T (t),y T (t),z T (t) The nominal planned path corresponding to time t is described as:

[0040]

[0041] Among them, the superscript T represents the transposition, and the subscript T is used to represent the current;

[0042] The path is divided into several segments based on the key points of the receiving aircraft path, and t is defined. i (i=1,2,…,j) represents the key point number, t i represents the time corresponding to the i-th key point, j represents the total number of key points, t0 represents the starting time, then Described as:

[0043]

[0044] Minimize the fourth-order derivative of the drone's position to obtain a smooth path for each segment (i.e., a path segment, which is the section of road between two path nodes):

[0045]

[0046] Among them, the drone position r=[x T (t),y T (t),z T (t)] T .

[0047] Furthermore, step 4 also includes time allocation optimization:

[0048] Let the time distribution vector τ=[T1,T2,…,T j ], where T i (i=1,2,…,j) represents the time length allocated to the i-th segment;

[0049] An initial value is set for the time allocation vector τ, and then the gradient descent method is used to iteratively optimize τ until the preset iterative convergence condition is reached, and the final time allocation vector τ is obtained.

[0050] Furthermore, in step 4, the flat output minimum mutation algorithm is used to calculate the nominal planning path When optimizing, set the safety corridor constraint and define the distance from the point on the path at time t to the corresponding segmented straight line segment as:

[0051] d i (t)=(r T (t)-r i )-((r T (t)-r i )·s i )s i (9)

[0052] Among them, r T (t) is the path at time t, r i is the i-th key point, s i For along r i Point to r i+1 The unit vector of

[0053] Define the width of the safety corridor as ε i , requiring that the distance that the path deviates from the straight line segment does not exceed the width of the safety corridor, then the safety corridor constraint is described as:

[0054] |d i (t)|≤ε i ,t i ≤t≤t i+1 (10)

[0055] All key points of the path are taken from the straight line path. Therefore, the path at the key points must satisfy the safety corridor constraint.

[0056] The denser the key points, the closer the path will be to a straight line. Therefore, by introducing intermediate points as new key points, the path can be made closer to a straight line, avoiding collisions with obstacles caused by deviations from the straight line.

[0057] By introducing intermediate points in the segments that do not meet the constraints and re-optimizing the path, after several iterations, the path can be ensured to be within the safe corridor. The path after introducing the intermediate points must meet the following requirements:

[0058]

[0059] Among them, n c To introduce the number of intermediate points, the final path generated after optimization is obtained.

[0060] The technical solution provided by the present invention brings at least the following beneficial effects:

[0061] (1) For UAV aerial refueling, the present invention can plan a safe, smooth and optimal path under the influence of vortex.

[0062] (2) The present invention uses two-layer planning of front-end and back-end. The front-end uses the KHA* algorithm combined with the aerodynamic analysis results to plan an initial path. The back-end uses the FMS algorithm, adds safety corridor constraints, and optimizes the time and smoothness of path planning to obtain the optimal path, ensuring safety while improving planning efficiency.

[0063] (3) The present invention is simple and easy to use, and is convenient for engineering implementation. BRIEF DESCRIPTION OF THE DRAWINGS

[0064] In order to more clearly illustrate the technical solutions in the embodiments of the present invention, the following briefly introduces the drawings required for use in the description of the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0065] Figure 1 The distribution diagram of Δcmx, Δcmy, Δcmz and K of the receiving aircraft when the span is 1.5 times, and the black dots indicate the corresponding extreme values.

[0066] Figure 2 The distribution diagram of Δcmx, Δcmy, Δcmz and K of the receiving aircraft when the span is 3 times, and the black dots indicate the corresponding extreme values.

[0067] Figure 3 The distribution diagram of Δcmx, Δcmy, Δcmz and K of the receiving aircraft when the extension is 5 times, and the black dots indicate the corresponding extreme values.

[0068] Figure 4 The distribution diagram of Δcmx, Δcmy, Δcmz and K of the receiving aircraft when the extension is 10 times, and the black dots indicate the corresponding extreme values.

[0069] Figure 5 The distribution diagram of Δcmx, Δcmy, Δcmz and K of the receiving aircraft when the extension is 15 times, and the black dots indicate the corresponding extreme values.

[0070] Figure 6 Schematic diagram of the vortex distribution generated by the tanker wing.

[0071] Figure 7 To customize the 3D environment map based on the eddy distribution, the eddies are modeled as cylinders in the figure.

[0072] Figure 8 This is the force analysis diagram of the fixed-wing UAV.

[0073] Figure 9 This is the initial path diagram planned in the vortex map using the KHA* algorithm. The path is represented by a solid line in the figure, and the key points are represented by “*” points in the figure.

[0074] Figure 10 This is the path map after the initial path in the vortex map is optimized by the FMS algorithm.

[0075] Figure 11 This is the collision diagram between the path and the obstacle when there is no safety corridor constraint.

[0076] Figure 12 This is the path diagram after adding the safety corridor constraint.

[0077] Figure 13 This is a graph showing the change in pitch angle during path planning.

[0078] Figure 14 This is a graph showing the change in roll angle during path planning. DETAILED DESCRIPTION

[0079] In order to make the purpose, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be described in detail and completely in conjunction with the drawings in the implementation of the present invention. Obviously, the described embodiments are only part of the embodiments of the present application, not all of the embodiments. Generally, the components of the embodiments of the present invention described and shown in the drawings can be arranged and designed using different configurations. Therefore, the following detailed description of the embodiments of the present invention provided in the drawings is not intended to limit the scope of the present application for protection, but merely represents selected embodiments of the present invention.

[0080] The present invention proposes a close formation path planning method for aerial refueling, which is used to reduce the adverse effects of tanker vortexes during aerial refueling and ensure the safety, smoothness, and economy of the UAV's path during aerial refueling. The concept of the embodiment of the present invention is as follows: first, the two stages of forming a formation are analyzed, including the rendezvous stage and the formation entry stage; then, based on wind tunnel test data, the influence of vortices on the receiving aircraft at different positions is calculated, thereby obtaining the key points of the receiving aircraft's path in the two stages of aerial refueling, and the vortex negative impact area is defined as an obstacle body, a three-dimensional obstacle model for path planning is constructed, and an aerial refueling path planning map is obtained based on this; different constraints are determined for the situations where vortices affect the rendezvous stage and the formation entry stage, and the key point heuristic A* algorithm (KHA*) is used to perform initial planning in the map based on the obtained path key points; a mathematical model of a fixed-wing UAV is established, and the differential flatness characteristics of the fixed-wing UAV are analyzed, so that the position of the UAV and the various derivatives of the position are used to describe other state quantities and control quantities of the UAV system; the differential flatness characteristics are used to establish a cost function for the path optimization process, and the flat output minimum mutation algorithm (FMS) is used in combination with the constraints to optimize the path obtained by the initial planning to obtain the final aerial refueling path.

[0081] As a possible implementation, the embodiment of the present invention provides a tight formation path planning method for aerial refueling, which specifically includes the following steps:

[0082] Step 1: Analyze wind tunnel test data to identify key points in the aerial refueling process and build a path planning map;

[0083] Based on wind tunnel test data, the aerodynamic data of the receiving aircraft in the two stages of aerial refueling (rendezvous stage and formation stage) were analyzed. The aerodynamic data were used to calculate the impact of vortices on the receiving aircraft at different positions, including the changes in the lift coefficient, drag coefficient, roll, pitch and yaw moment coefficients of the receiving aircraft at different positions in the vortex. The key points of the receiving aircraft's path in the two stages of aerial refueling were obtained, and the areas with greater negative impacts of vortices were defined as obstacles. A three-dimensional model for path planning was constructed and defined as a map for path planning.

[0084] The two stages of aerial refueling are analyzed, and the influence of vortex exists in both stages:

[0085] 1) Meeting stage:

[0086] Before entering the formation, the tanker and the receiving aircraft need to prepare for rendezvous. First, the tanker should arrive and cruise along the holding pattern at the lowest fuel consumption rate. At this time, the receiving aircraft needs to gradually approach from behind the tanker's left wing, always maintaining a vertical separation of a standard altitude (such as at least 305 meters). When the longitudinal separation between the two sides is reduced to the standard separation (such as 1852 meters), the two aircraft should achieve visual contact to verify each other's identity. Afterwards, the receiving aircraft joins the refueling route and climbs steadily, gradually approaching the tanker and reaching the rendezvous phase target point (also known as the rendezvous phase end point).

[0087] 2) Formation phase:

[0088] When the longitudinal distance between the two sides is reduced to within the target distance (such as 500 meters), the receiving aircraft departs from the end of the rendezvous phase, and after confirming that there are no other aircraft behind the tanker, it joins the formation from the left. After reaching the desired position behind the tanker, the formation phase is completed.

[0089] Based on the wind tunnel test data, the aerodynamic analysis of the receiving aircraft in the two stages of aerial refueling (rendezvous stage and formation stage) was carried out. The wingspan of the receiving aircraft was defined as b. The key points of the path were selected based on the data of five cross-sections with the longitudinal distance X between the two aircraft being 1.5, 3.0, 5.0, 10.0 and 15.0 times the span. The key points are the positions where the lift-to-drag ratio is the largest and the moment coefficient changes the least. The data analysis is as follows: Figures 1 to 5 The data for the formation phase is within five times the longitudinal span, while the data for the rendezvous phase is greater than five times the span. As shown in the table below, K is the lift-to-drag ratio, Δcmx, Δcmy, and Δcmz are the coefficient changes in roll, pitch, and yaw moments, respectively. X, Y, and Z are the relative longitudinal, lateral, and vertical distances between the two aircraft, respectively. The unit of measurement is b = 50m.

[0090] Table 1: Longitudinal spacing is 1.5 times the extension length

[0091]

[0092] Table 2: Longitudinal spacing is 3 times the extension length

[0093]

[0094] Table 3 Longitudinal spacing is 5 times the extension length

[0095]

[0096]

[0097] Table 4 Longitudinal spacing is 10 times the extension length

[0098]

[0099] Table 5 Longitudinal spacing is 15 times the extension length

[0100]

[0101] It can be seen from the above table that when the longitudinal spacing between the two machines is 1.5 times the span, Y=-0.375b and Z=0.14b are both the positions with the maximum gain and the positions with relatively small changes in the moment coefficient; when the longitudinal spacing between the two machines is 3 times the span, the position with the maximum gain is also the position with the smallest change in Δcmx. For Δcmy and Δcmz, the positions with the smallest changes are around z=-0.25b and z=0.043b respectively; when the longitudinal spacing between the two machines is 5 times the span, the position with the maximum gain is z=-0.21b, the positions with the smallest changes in Δcmx and Δcmy are almost the same at z=-0.14b, and the position with the smallest change in Δcmz is around z=-0.04b.

[0102] When the distance between the two aircraft is between 1.5b and 5b, it is found that when Y = -0.375b and Z = -0.12b, the three situations all obtain the maximum gain at this time. At the same time, as the longitudinal spacing increases, the rolling moment coefficient reaches the minimum value at 3 times the extension, the pitch moment gradually decreases, and the yaw moment does not change significantly.

[0103] When the longitudinal distance between the two machines is 10 times the span, the position with the maximum gain is also the position with the minimum disturbance of Δcmy. For Δcmx and Δcmz, the positions with the maximum gain and the minimum disturbance are conflicting, especially for Δcmx, which is completely opposite.

[0104] When the longitudinal distance between the two machines is 15 times the span, the position with the maximum gain is also the position with the minimum disturbance of Δcmy and Δcmx, so the position with the maximum gain can be used directly.

[0105] From this, the cost function W can be established as:

[0106]

[0107] where λ i , i∈[1,2,…,4] represents the weight corresponding to each item. According to the degree of change of the four variables, the parameters are respectively set as λ1=λ2=λ3=1,λ4=0.001, and the minimum cost W min The corresponding position is the path key point.

[0108] According to the aerodynamic data at different cross-sectional positions, the corresponding minimum cost positions at different cross-sectional positions are obtained, as shown in Table 6. The unit in the table is span b.

[0109] Table 6 Key point positions

[0110]

[0111] At the same time, the area with the greatest negative impact of vortex is obtained and modeled as an obstacle. The schematic diagram of the vortex distribution generated by the tanker wing is shown in the figure below. Figure 6 As shown in the figure, we can see the distribution of vortices behind the wing. By analyzing the vortices, we can define the area with the greatest impact as a three-dimensional cylindrical obstacle, and build a three-dimensional model for path planning as shown in the figure below. Figure 7 As shown, this model is defined as a map for path planning.

[0112] Step 2: Determine the constraints for the two phases of aerial refueling and use the KHA* algorithm for initial path planning;

[0113] During the path planning process, it is necessary to consider the path length constraints, yaw angle constraints, pitch angle constraints, and relative distance constraints between the two aircraft. Based on the key path points obtained in step 1, without considering the kinematic and dynamic constraints of the aircraft itself, the KHA* algorithm is used to perform initial path planning in the constructed 3D environment model to obtain an initial path without considering any constraints.

[0114] During path planning, the following constraints need to be considered:

[0115] 1) Path length constraint:

[0116] The path length is used to represent the flight distance of the UAV. Let L max is the maximum flight distance of the UAV. Then, the path length constraint is calculated as

[0117]

[0118] Among them, (x i ,y i ,z i ),i=1,2,...,N,is the coordinate of the i-th path point in the three-dimensional planning environment, and N is the total number of path points.

[0119] 2) Yaw angle constraint:

[0120] Considering the smoothness of the path, the maximum turning angle is introduced to constrain the turning angle, that is, the yaw angle ψ. Assume that vector a i Equal to [x i -x i-1 ,y i -y i-1 ,] T ,||ai || is vector a i The yaw angle constraint of the UAV is

[0121]

[0122] Among them, ψ max Indicates the maximum yaw angle of the drone.

[0123] 3) Pitch angle constraint:

[0124] Similar to the turning angle, the climbing angle is the pitch angle, and the pitch angle θ can be calculated as:

[0125]

[0126] Among them, θ max is the maximum pitch angle, -θ max is the minimum pitch angle.

[0127] 4) Relative distance constraints between two machines:

[0128] The drone must maintain the relative distance after changing the flight attitude or turning, and will not turn frequently. Therefore, it is necessary to limit the path segment between adjacent path points to be greater than the shortest path segment length L min . The path segment constraint can be expressed as

[0129] L min -||p i+1 -p i ||≤0 (16)

[0130] Among them, p i+1 Yes (x i+1 ,y i+1 ,z i+1 ), p i Yes (x i ,y i ,z i ), p i+1 and p i are two adjacent waypoints.

[0131] The KHA* algorithm considers the cost of moving from the starting point to the node and the estimated cost of moving from the node to the end point, as well as the critical path points obtained in step 1, without considering its own constraints. Before path planning, the map must be simplified into a set of quantifiable nodes. The map area is divided into a number of (three-dimensional) elements using a specific step size. To further improve path planning accuracy, the elements are typically smaller, using a smaller step size (which can be set based on the actual scenario to ensure a sufficient number of elements). These elements are referred to as small elements. Each element records the location of its center of mass and is manually labeled as either an obstacle element or a non-obstacle element based on whether its center of mass is within an obstacle. Taking into account the element's volume and the drone's radius, the obstacle is widened so that more elements are recorded as being within the obstacle. This ensures that when the drone's center of mass enters an element marked as non-obstacle, it will not collide with the obstacle due to the influence of the drone's volume and the element's volume.

[0132] After the map is divided, the path can be planned. The KHA* algorithm starts from the rendezvous starting point s to the docking destination. Each search evaluates the current search position through the cost function C(n). The cost function C(n) is:

[0133] C(n)=p(n)+q(n) (17)

[0134] Where n represents the nth node, p(n) is the actual cost function from the initial node s to the current node n, and q(n) is the estimated cost function from the current node n to the end node o (i.e., the movement cost from the current node to the end node o, such as the Manhattan distance between the current node and the end node o). Combined with the key points of aerodynamic constraints, they serve as the heuristic function of the A* algorithm.

[0135] Then, the accumulated cost function p(n e ) and q(n e ):

[0136] p(n e )=p(n)+Cnn e (18)

[0137]

[0138] Among them, x ,o y ,o z ,n x ,n y ,n z are the coordinates of the target point and the current point respectively, Cnn e From node n to node n eStarting from the starting point, all key frames and key points constitute the path nodes. Each round of search between adjacent nodes is based on the minimum cost function C(n). Each point is connected in sequence to obtain a nominal planning path from the rendezvous starting point to the docking end point, which is recorded as

[0139] In this embodiment, the processing process of the heuristic A* algorithm includes:

[0140] The input parameters are the initial node s (starting point) and the end node o (target point), and the output result is a nominal planning path connecting the initial node s and the end node o.

[0141] First, initialize the priority queue as a container for storing the nodes to be expanded, and set the actual cost function value of the initial node s to 0, and the actual cost function values ​​of the remaining nodes to infinity.

[0142] When the priority queue is not empty, the loop process begins. In the loop, the node n with the smallest C(n) value in the queue is first removed and marked as expanded. If node n is exactly the target node o, the loop terminates. For all unexpanded neighbor nodes n of node n, e If n is found e The actual cost function value is still infinite, indicating that n e It has not been expanded yet. At this time, based on the actual cost function value of node n (if node n is the initial node s, its value is 0; otherwise it is the value after the latest update), calculate p(n) according to formulas (18) and (19) e ) and q(n e ), that is, the cost function p(n e ) and q(n e ), and n e Join the queue. If a neighbor node n is found e p(n e ) and q(n e ) is greater than the sum of p(n) and Cnn e The sum of n e Has been expanded, now needs to update n e The current actual cost function value is p(n)+Cnn e .

[0143] When all nodes in the queue have been expanded, the sequence of nodes traversed by the algorithm is arranged in reverse order to obtain the nominal planning path.

[0144] Step 3: Establish a mathematical model of the fixed-wing UAV and analyze the differential flatness characteristics of the fixed-wing UAV;

[0145] A mathematical model of fixed-wing UAV is established, and the differential flatness characteristics of the fixed-wing UAV are analyzed. The state and control quantities represented by flat outputs are obtained, that is, the state and control quantities of the UAV system can be described by the position of the UAV and its various position derivatives.

[0146] Define the coordinate system as a right-handed coordinate system, expressed as G a ={i a ,j a ,k a}, where i a ,i a and k a is the three-dimensional (x, y, z) basis vector. The expression In the coordinate system G b Measured in coordinate system G a The vector x in the coordinate system G c In the statement, Indicates that the coordinate system G a The coordinates in are converted to coordinate system G b The rotation matrix of the coordinates in . The regular basis vectors in the (three-dimensional real field) are denoted as e1, e2 and e3. Coordinate system G i ,G b and G s They represent the inertial reference system, the body reference system, and the stable reference system. Among them, the inertial reference system can also be called the ground coordinate system, which is used for orientation indication and navigation control. Its origin is located at a point on the ground, and the x-axis (i.e., a ) is located in the ground plane and points in a specific direction. The y-axis (i.e., j a ) is in the same plane as the x-axis and is perpendicular to the x-axis. Its direction follows the right-hand rule. The z-axis (i.e., k a ) vertically downward, pointing to the center of the earth. The body reference system is used to describe the inherent motion of the drone, with its origin located at the center of mass of the drone, the x-axis (i.e., b ) is in the symmetry plane of the UAV, parallel to the fuselage axis of the UAV or the average aerodynamic chord of the wing, pointing to the direction of the UAV's nose, and the y-axis (i.e., j b ) axis is perpendicular to the symmetry plane of the drone and points to the right side of the fuselage. The z axis (i.e., k b ) is perpendicular to the x-axis and points downwards. The stable reference frame is used to describe the motion characteristics of the drone after a small disturbance. Its origin is at the center of mass of the drone, and the x-axis (i.e. s ) coincides with the projection of the UAV airspeed in the symmetry plane, and the y-axis (i.e., j s ) and the body coordinate axis j b Coincident, z axis (ie k s ) in the symmetry plane of the aircraft and i s Vertically and pointing down the fuselage.

[0147] set up In the generalized vector The skew-symmetric matrix operator is defined above, such as

[0148]

[0149] The mathematical model of fixed-wing UAV is established as follows, and the force analysis of UAV is as follows: Figure 8 As shown:

[0150]

[0151] Where O is the position of the "unmanned aerial system", and the i, b, and s in the subscripts represent the variables defined in the inertial coordinate system, the body coordinate system, and the stable coordinate system, respectively. Denotes the stable coordinate system G s To the inertial coordinate system G i The stability-inertia selection matrix is: v is the velocity of the UAV, g is the acceleration of gravity, m is the mass of the system, and the angular velocity ω of the UAV is: Represents the coordinate system G b The position O in the coordinate system G i In the expression The first derivative of is the coordinate system G b The velocity in the coordinate system G i In the expression, for The first derivative of In the coordinate system G i Measurement coordinate system G b The angular velocity in the coordinate system G b In the statement, Represents the coordinate system G b To coordinate system G i The body-inertia selection matrix, for The first derivative of . That is, the symbol " " on the parameter head represents its first derivative, and the symbol " " represents its second derivative. ω c is the specified angular rate, α is the angle of attack, The drone is along the negative body axis i b resistance, The UAV is along the body axis i b The thrust generated, The UAV is along the negative body axis k b To simplify the model, it is assumed that the deviations in rudder deflection, sideslip angle, and angle of attack α0 are relatively small and the drag caused by the control surfaces can be ignored. Using these assumptions, the drag term can be simplified to

[0152]

[0153] Where ρ is the air density, is the airspeed of the drone, S is the plane area of ​​the wing, c is the chord length of the wing, and is the aerodynamic coefficient.

[0154] Assuming no wind, the drone's airspeed is

[0155] Assuming that there is no wind during the flight of the UAV, the angle of attack is constant and equal to the trim angle of attack α0, and the sideslip angle is zero, the state and input of the UAV system are defined as follows:

[0156]

[0157] Let the stability-inertia rotation matrix Expressed as

[0158]

[0159] Among them, r j is the jth stable coordinate axis of the fixed wing, expressed in an inertial coordinate system. According to the above two assumptions, the speed of the drone is related to the stable axis i s Alignment, i.e.

[0160]

[0161] in, is the z-axis of the body reference system (z b ) direction vector z b The first derivative of .

[0162] Accordingly, we can get

[0163]

[0164] Taking the dot product of r1 gives

[0165]

[0166] Accordingly, we can get

[0167]

[0168] Since ||r3||=1, we can conclude that

[0169]

[0170] There must be

[0171] r2=r3×r1 (30)

[0172] Then the body-inertia rotation matrix is ​​calculated using the stability-inertia rotation matrix as

[0173]

[0174] in

[0175]

[0176] Rotation Matrix Can be written as flat output Function of , angular velocity vector:

[0177]

[0178] Define the “∨” operator as

[0179]

[0180] The angular velocity can be expressed in terms of the plane output as

[0181]

[0182] The thrust T can be obtained as

[0183]

[0184] Using V a and is defined, and Give available

[0185]

[0186] It can be concluded that all the states x and control inputs u of the system can be described by the flat output z and its derivatives. That is, and

[0187] Then the system is differentially flat, and the flat output can be obtained as:

[0188]

[0189] It can be seen that all states and inputs of the fixed-wing UAV can be defined by the position of the UAV and its derivatives.

[0190] Step 4: Using the differential flatness characteristics of the fixed-wing UAV, a path function for the path optimization process is established, and the FMS algorithm is used to optimize the initial path to obtain the final path;

[0191] Combined with the differential flatness characteristics of the fixed-wing UAV in step 3, it can be obtained that the state and control quantities of the UAV system can be described by the position of the UAV and its various position derivatives. This feature is used to establish the cost function of the path optimization process. Considering the path smoothness, flight time, safety and the constraints in step 2, the FMS algorithm is used to optimize the initial path generated in step 2 to obtain the final path.

[0192] Based on the analysis of the differential flatness characteristics of fixed-wing UAVs, the path of fixed-wing UAVs It can be described by the following formula:

[0193]

[0194] Among them, x T (t),y T (t),z T (t) is the coordinate position of the UAV in space. In order to facilitate the optimization of the path, a set of key points is used to divide the path into several segments, i represents the key point number, t i represents the time corresponding to the i-th key point, and j represents the total number of key points.

[0195]

[0196] In the process of path planning for fixed-wing UAVs, in order to make the path of the UAV as smooth as possible, the fourth-order derivative of the position must be minimized to obtain a smooth path. Specifically for a certain segment i, that is:

[0197]

[0198] where r = [x(t), y(t), z(t)] T is the coordinate position of the drone in space.

[0199] The three parameters of the flattened output are decoupled from each other. Therefore, the three parameters (three-dimensional coordinates) of the flattened output can be optimized separately to achieve flattening optimization.

[0200] Taking parameter x as an example, we perform path optimization. We use an n-order polynomial with time t as the independent variable to describe the path of a certain segment i of a certain parameter of the drone, namely:

[0201] p i (t) = p i0 +p i1 t+p i2 t 2 +…+p in t n ,t i-1 ≤t≤t i (42)

[0202] Define the polynomial coefficient vector P i as follows:

[0203] P i =[p i0 ,p i1 ,p i2 ,...,p in ] T (43)

[0204] Then the polynomial path can be expressed as:

[0205] p i (t)=[0,t,t 2 ,…,t n ]P,t i-1 ≤t≤t i (44)

[0206] P is the polynomial coefficient matrix. If the position of the drone is known, p(t), the velocity, acceleration, jerk, and snap of the fixed-wing drone can be obtained by taking derivatives of each order (jerk is the third-order derivative of position with respect to time, and snap is the fourth-order derivative of position with respect to time). The overall path is divided into j segments according to key points, which can be written as:

[0207]

[0208] P i =[p i0 ,p i1 ,p i2 ,…,p in ] T ,i=1,2,…,j (46)

[0209] You can Rewritten in matrix form in the direction of the parameters currently being solved:

[0210]

[0211] Among them, the matrix Matrix Q i Block Q in the lower right corner i,right The elements satisfy: Where r' is the row index, c' is the column index, and r'≥5,c'≤n+1.

[0212] According to the matrix of the entire path, the form of the entire path can be obtained:

[0213]

[0214] Among them, the matrix

[0215] The optimization problem can be written as follows:

[0216] minP T QP (49)

[0217] stAP=b (50)

[0218] Vector b stores the values ​​of various parameters at the keyframes, which correspond one-to-one to the constraints in matrix A. AP = b is the optimization constraint. Taking the i-th segment as an example, the constraints include:

[0219] 1) The drone is at the end of the segment t i At this moment, its position p i , speed v i , acceleration a i , jerk i 、snap i is a specific value, forming 5 equality constraints;

[0220] [1,t i ,t i 2 ,…,t i n ]P i =p i (51)

[0221] [0,1,2t i ,3t i 2 ,4t i 3 …,nt i n ]P i =v i (52)

[0222] [0,0,2,6t i ,12t i 2 ,…,n(n-1)t i n-2 ]P i =a i (53)

[0223]

[0224] 2) The parameters of the UAV should not change suddenly at the segment boundaries. That is, the difference between the starting time of the i-th segment and the end time of the i-1-th segment should be 0, forming five equality constraints, called continuity constraints.

[0225] [1,ti-1 ,t i-1 2 ,…,t i-1 n ]P i -[1,t i-1 ,t i-1 2 ,…,t i-1 n ]P i-1 =0 (56)

[0226] [0,1,2t i-1 ,3t i-1 2 ,4t i-1 3 …,nt i-1 n ]P i -[0,1,2t i-1 ,3t i-1 2 ,4t i-1 3 …,nt i-1 n ]P i-1 =0(57)

[0227] [0,0,2,6t i-1 ,12t i-1 2 ,…,n(n-1)t i-1 n-2 ]P i -[0,0,2,6t i-1 ,12t i-1 2 ,…,n(n-1)t i-1 n-2 ]P i-1 =0 (58)

[0228]

[0229] These 10 equality constraints are written in matrix form as AP=b, which means that the polynomial coefficients of the target path can be solved to obtain a path described by the polynomial.

[0230] Optimizing the time allocation in the planning process can prevent the paths from overlapping. The time allocation vector τ = [T1, T2, ..., T j ], where T i , (i=1,2,…,j) represents the time length assigned to the i-th segment. Then the time t corresponding to each key point is i There are the following relationships:

[0231]

[0232] First, an initial value for τ is set, and then gradient descent is used to iteratively optimize τ. There are several ways to set this initial value: ① Assign each segment the same time length as the starting value. ② Allocate initial time values ​​proportionally to the length of the straight paths within each segment. ③ Taking into account acceleration and deceleration, allocate time based on the length of each segment path to reflect the actual flight speed changes.

[0233] Taking the time length of each segment as the independent variable, define the cost function J:

[0234]

[0235] Find the gradient of J Optimize τ using the gradient descent algorithm:

[0236]

[0237] Among them, β represents the step size of each iteration, τ k represents the time allocation vector for the k-th iteration.

[0238] Repeat the above process iteratively until the cost function meets certain indicators, such as the cost function is less than a given value, or the change in the cost function between two iterations is less than a certain proportion, then the optimal segmentation time optimization is considered to be completed.

[0239] In order to prevent collisions during the planning process, it is necessary to add safety corridor constraints, as shown in formulas (10) and (11).

[0240] The steps are verified through examples. After data analysis, the key points of the receiving aircraft's path during aerial refueling and a three-dimensional model for path planning are obtained. The KHA* algorithm is used to perform initial planning in the map to obtain an initial path. Analysis of the differential flatness characteristics of fixed-wing UAVs shows that the state and control quantities of the UAV system can be described by the UAV's position and its various position derivatives. This characteristic is used to establish a cost function for the path optimization process. The FMS algorithm is used to optimize the path obtained by the initial planning in combination with constraints to obtain the final path.

[0241] Example

[0242] In the MATLAB simulation, the KHA* algorithm is used for unconstrained initial path planning. The map is modeled according to the existence of eddy currents. The map is as follows: Figure 7 The basic parameters are shown in Table 7, and the obstacle is the cylinder in the figure.

[0243] Table 7 Basic parameters of planning map

[0244]

[0245] Input the map and the starting and ending point information into the KHA* algorithm module to obtain an initial path as follows Figure 8 As shown, each point on the path records the three-dimensional coordinates of the point. Figure 8 As can be seen, the path generated using KHA* has large turns at corners and is not smooth enough. The path is represented by a solid line in the figure. The starting point is used as the first keypoint, and the next point is searched sequentially. The line connecting the keypoint and the next target point is checked to see if it collides with an obstacle. If not, the next point is checked. If so, the previous point is used as the new keypoint, and the same method is used to find the next keypoint. Keypoints are represented by "*" dots in the figure. The coordinates of the six keypoints are shown in Table 8:

[0246] Table 8 Key point numbers and coordinates

[0247]

[0248] After obtaining the key points, a polyline path was generated. However, the polyline path was not smooth and had large turning angles, necessitating FMS optimization. The squared norm of the fourth-order derivative of position was used as the cost function. Constraints at different stages were also considered, with the constraint parameters set as shown in Table 9. Based on the constraints being satisfied, optimization was performed with the goal of minimizing the cost function. This resulted in a segmented path, each described by a polynomial. After the path was generated, each segment was checked for collisions with obstacles. A certain number of points were sampled along the path at a certain sampling rate, and the position of each point was calculated. Based on the previously loaded map information, it was determined whether the point was within an obstacle element. If any point on a path segment was within an obstacle element, a safe corridor constraint was marked for that segment and added. After adding the constraint, the path was regenerated and collision detection was repeated. If collisions still occurred, a narrower safe corridor width was added until no collisions occurred. To avoid path overlap, the path was time-optimized to obtain the final path.

[0249] Table 9 Smoothing optimization constraint parameter settings

[0250]

[0251] Optimize the generated initial path and get the optimized path as follows Figure 10 As shown, it can be seen Figure 10 and Figure 9 It is smoother than in corners; Figure 11 、 Figure 12This is a detailed diagram of adding a safety corridor constraint when a path segment collides with an obstacle. It can be seen that when the safety corridor constraint is not added ( Figure 11 ) The planned path collides with the obstacle. The collision path segment is between the second key point and the third key point. After adding the safety corridor constraint optimization, the collision is avoided. Figure 12 shown. Figure 13 、 Figure 14 They are the changes of pitch angle and roll angle during the path planning process. It can be seen from the figure that the maximum values ​​of pitch angle and roll angle do not exceed the given constraints and are both within the constraint range.

[0252] Experimental results show that by optimizing the initial path while satisfying basic constraints and adding safety corridor constraints, a safe and smooth optimal path can be obtained.

[0253] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent replacements for some of the technical features therein. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the various embodiments of the present invention.

[0254] The above are only some embodiments of the present invention. For those skilled in the art, several modifications and improvements can be made without departing from the inventive concept of the present invention, which all fall within the scope of protection of the present invention.

Claims

1. A close formation path planning method for aerial refueling, characterized in that: The method comprises the following steps: Step 1: Based on wind tunnel test data from the rendezvous and formation phases of aerial refueling, key points of the receiving aircraft's path during the rendezvous and formation phases are obtained according to the changes in the lift coefficient, drag coefficient, roll, pitch, and yaw moment coefficients of the receiving aircraft at different positions in the vortex. A three-dimensional obstacle model for path planning is built based on the selected obstacles to obtain an aerial refueling path planning map. Step 2: Determine the constraints of the two phases of aerial refueling. Based on the key points of the receiving aircraft's path in the rendezvous phase and the formation phase, use the key point heuristic A* algorithm to perform initial path planning in the aerial refueling path planning map to obtain the nominal planned path. Step 3: Establish a fixed-wing UAV model and analyze the differential flatness characteristics of the fixed-wing UAV to describe the state and control variables of the UAV system using the position of the UAV and its derivatives. Step 4: Based on the differential flatness characteristics of the fixed-wing UAV, the cost function of the path optimization process is established. Based on the constraints of the two stages of aerial refueling, the flat output minimum mutation algorithm is used to optimize the nominal planning path. Optimize and obtain the final aerial refueling path.

2. The method according to claim 1, wherein In step 1, the key points of the receiving aircraft's path in the rendezvous phase and the formation phase are determined based on the minimum value of the cost function W; Among them, the expression of the cost function W is: Wherein, Δcmx represents the rolling moment coefficient, Δcmy represents the pitching moment coefficient, Δcmz represents the yaw moment coefficient, K represents the lift-to-drag ratio, and λ1 to λ4 represent the weights of the rolling moment coefficient, pitching moment coefficient, yaw moment coefficient, or lift-to-drag ratio, respectively.

3. The method according to claim 1, wherein In step 2, the constraints include: path length constraint, yaw angle constraint, pitch angle constraint and / or relative distance constraint between the two aircraft.

4. The method according to claim 1 or 3, wherein Step 2 includes: The aerial refueling path planning map is divided into multiple three-dimensional elements according to a certain step size, and obstacle blocks and non-obstacle blocks are marked. Each three-dimensional element represents a path node. Then, the KHA* algorithm is used to trigger from the set starting point, gradually passing through all the key points. Each node found between two key points must satisfy the minimum movement cost function, and each path node on the path is found in turn, and then connected with the key points to form the nominal planning path. The movement cost function includes the movement cost from the starting point to the path node and the estimated movement cost from the path node to the end point.

5. The method according to claim 1, wherein Step 3 includes: Define coordinate system G i ,G b and G s They represent the inertial reference frame, the body reference frame and the stable reference frame respectively; Among them, the inertial reference system G i The origin is located at a point on the ground, x i The axis lies in the ground plane and points in a specific direction, y i Axis and x i The axes are in the same plane and perpendicular to x i Axis, z i The axis points vertically downward, toward the center of the Earth; Body reference frame G b The origin is located at the center of mass of the drone, x b The axis is in the symmetry plane of the UAV, parallel to the fuselage axis of the UAV or the average aerodynamic chord of the wing, and points to the direction of the UAV's nose. b The axis is perpendicular to the symmetry plane of the drone and points to the right side of the fuselage. b with x b The axis is vertical and points downwards from the drone; Stable reference frame G s The origin is located at the center of mass of the drone, x s The axis coincides with the projection of the drone's airspeed in the symmetry plane, y s Axis and body coordinate axis y b Overlap, z s In the symmetry plane of the aircraft and x s Vertical and pointing down the fuselage; Define the regular basis vectors in the three-dimensional coordinate system as e1, e2 and e3, and establish the fixed-wing UAV model as: The fixed-wing UAV model is established as: in, Represents the coordinate system G b The position O in the coordinate system G i In the expression The first derivative of is the coordinate system G b The velocity in the coordinate system G i In the expression, for The first derivative of g and m are the acceleration of gravity and the mass of the receiving aircraft respectively. Denotes the stable coordinate system G s To the inertial coordinate system G i The stability-inertia selection matrix is ​​α, which is the angle of attack, and D is the UAV's motion along the negative x direction. b The drag force of the UAV along the x-axis, T is b The thrust generated by the axis, L is the thrust generated by the UAV along the negative z axis. b The lift generated by the shaft, In the coordinate system G i Measurement coordinate system G b The angular velocity in the coordinate system G b In the statement, Represents the coordinate system G b To coordinate system G i The body-inertia selection matrix, for The first derivative of ω c is the specified angular rate.

6. The method according to claim 5, wherein Set the drag in the fixed-wing drone model to: Where ρ is the air density, V a is the airspeed of the drone, S is the plane area of ​​the wing, c is the chord length of the wing, and is the aerodynamic coefficient.

7. The method according to claim 5 or 6, wherein: In step 3, the state x and control input u of the UAV system are defined as: The state x and control input u are described in terms of the flat output z and its high-order derivatives as:

8. The method according to claim 1, wherein In step 4, the cost function of the path optimization process is established based on the differential flatness characteristics of the fixed-wing UAV, including: According to the differential flatness characteristics of fixed-wing UAV, based on the current coordinate position x of the UAV in space T (t),y T (t),z T (t) The nominal planned path corresponding to time t is described as: Wherein, the superscript T indicates transposition; The path is divided into several segments based on the key points of the receiving aircraft path, and t is defined. i (i=1,2,…,j) represents the key point number, t i Indicates the time corresponding to the i-th key point, j indicates the total number of key points, t0 indicates the starting time, Described as: Minimize the fourth-order derivative of the drone's position to obtain a smooth path for each segment: Among them, the drone position r=[x T (t),y T (t),z T (t)] T .

9. The method according to claim 8, wherein Step 4 also includes time allocation optimization: Let the time distribution vector τ=[T1,T2,…,T j ], where T i (i=1,2,…,j) represents the time length allocated to the i-th segment; An initial value is set for the time allocation vector τ, and then the gradient descent method is used to iteratively optimize τ until the preset iterative convergence condition is reached, and the final time allocation vector τ is obtained.

10. The method according to claim 8, wherein In step 4, the flat output minimum mutation algorithm is used to calculate the nominal planning path When optimizing, set the safety corridor constraint as: |d i (t)|≤ε i ,t i ≤t≤t i+1 Among them, d i (t) represents the distance from the point on the path to the corresponding segmented straight line segment at time t, d i (t)=(r T (t)-r i )-((r T (t)-r i )·s i )s i , r T (t) is the path at time t, r i is the i-th key point, s i For along r i Point to r i+1 The unit vector of ε i Indicates the preset safety corridor width.

Images