A trajectory planning method considering low-altitude wind and UAV energy consumption

Abstract

The present invention relates to a method for trajectory planning of a UAV when taking low-altitude wind and UAV energy consumption into consideration when flying in mountainous areas, and belongs to the field of UAV trajectory planning methods. The energy consumption cost of the UAV is added to the objective function of the trajectory planning of the present invention, and a calculation model of the UAV energy consumption is given. The energy consumption constraint of the UAV is added to the model constraint, and the impact of low-altitude wind on the trajectory planning is considered to correct the trajectory, that is, the passing time of the UAV is corrected according to the vector relationship. In terms of algorithm, the model is solved by improving the A* algorithm. A comparative analysis is performed on the algorithm with that before the improvement, and the trajectory planning model considering low-altitude wind and energy consumption is solved with the improved algorithm. The present invention greatly improves the practicality and accuracy of the research on UAV trajectory planning when considering the actual operating environment.

Description

Technical Field

[0001] The invention relates to a track planning method taking low-altitude wind and energy consumption of unmanned aerial vehicles into consideration, and belongs to the field of unmanned aerial vehicle track planning methods. Background Art

[0002] The problem of UAV trajectory planning in low-altitude rescue is increasingly being solved by establishing a trajectory planning model that takes into account the actual operating environment and energy consumption. How to consider the actual operating environment and energy consumption has become a key issue in solving trajectory planning. The purpose of trajectory planning is to fully consider the energy consumption of UAVs and external environmental factors such as wind under the condition of meeting various constraints, so that trajectory planning is more practical. Trajectory planning mainly considers three aspects: obstacle avoidance based on terrain, trajectory correction based on low-altitude wind, and energy consumption as an optimization target.

[0003] (1) Trajectory planning considering obstacle avoidance

[0004] Obstacles are key factors that need to be considered in the trajectory planning of UAVs. Obstacles include static obstacles and dynamic obstacles. The offline trajectory planning problem mainly considers static obstacles, which are mainly mountain obstacles. Therefore, in order to simplify the problem, the offline trajectory planning problem is considered, that is, terrain obstacles are considered for trajectory planning. Offline trajectory planning methods include angle coding group optimization algorithm, A* algorithm, Voronoi method, artificial potential field method, genetic algorithm and other methods to obtain the trajectory of the UAV.

[0005] (2) Considering the impact of low-altitude wind on UAV trajectory planning

[0006] There are many external factors that affect UAVs in actual flight environments, but the impact of low-altitude wind on the flight of UAVs cannot be ignored. At present, the following methods have been used to consider the impact of low-altitude wind on trajectory planning: Zhang Ming et al.

[53] used a data fusion method to predict low-altitude winds and correct aircraft trajectories to address the problem that low-altitude rescue trajectories are easily affected by crosswinds and it is difficult to obtain accurate trajectory planning paths. They also verified the advantages and disadvantages of this method. Lu Yiyu et al.

[54] proposed that considering the influence of meteorological factors is an effective means to improve the accuracy of 4D trajectory prediction. Through the method proposed in this article, it was verified that the predicted trajectory is closer to the actual trajectory. Qu Yaohong et al.

[55] proposed a UAV trajectory planning method that uses combined navigation to estimate wind field information online. Experiments have verified that this method reduces flight time compared to traditional methods. Wu Anping et al.

[56] proposed an online trajectory planning method based on wind field estimation for the trajectory planning problem of small UAVs with angle constraints under wind field conditions. These methods have applied UAV trajectory planning to actual scenarios.

[0007] (3) Consider the impact of energy consumption on trajectory planning

[0008] It is generally believed that the shorter the distance a UAV flies, the lower its energy consumption. However, this evaluation method lacks rationality. The actual energy consumption of a UAV is closely related to the flight status and operating environment of the UAV. In order to make trajectory prediction more energy-efficient, the impact of energy consumption on trajectory planning must be considered. Currently, the following methods have been used to consider the impact of energy consumption on trajectory planning: Abeywickrama et al.

[57] proposed a complete UAV energy consumption model based on empirical research on battery performance. The model quantified the impact of ground power consumption, communication, hovering, vertical and horizontal movement, speed, payload and wind on the energy consumption of the UAV. However, this model requires testing the performance of the actual UAV. It is difficult to obtain an accurate model when conditions do not allow, and it is no longer applicable to other models other than the ones mentioned in the article. Tang et al.

[58] developed a method based on a power transfer model to estimate the real-time state of energy (SOE) and predict the end of discharge (EOD) time of a rotary-wing UAV under dynamic operating conditions. The effectiveness of this method was verified through real experiments on UAVs, but the influence of external factors such as wind was not considered. Liu et al.

[59] established an optimization model for minimizing UAV wireless energy consumption, but this energy consumption model did not take into account the wind factor. Gu et al.

[60] studied a UAV-assisted data dissemination communication system that took into account the data requirements of each IoT device, the limited energy budget of the UAV, the speed of the UAV, and the UAV’s transmission power in each time slot, but failed to consider the impact of wind on the UAV’s power consumption.

[0009] In summary, the research on obstacle avoidance in UAV trajectory planning in China is quite mature. However, in the previous trajectory planning problems, some only considered the influence of wind and the length of the trajectory, without mentioning the specific energy consumption. Although some considered both the energy consumption and the influence of wind, the proposed model is not universal. There are few related studies that consider both energy consumption and wind and propose a more general UAV power consumption model. In order to combine with actual scenarios, it is necessary to consider terrain obstacle avoidance, low-altitude wind and UAV energy consumption at the same time to make trajectory planning more practical. Summary of the invention

[0010] The present invention proposes a trajectory planning method that takes low-altitude wind and UAV energy consumption into consideration. This paper studies how to model the operating environment of the UAV and how to improve the trajectory solution algorithm, and then performs trajectory planning for the UAV. Among them, the initial trajectory of the UAV mainly considers the performance constraints of the UAV, including the maximum range of the UAV, the minimum track segment and the flight altitude limit, etc., and also considers the basic environmental constraints, namely terrain, threats and other factors to establish constraints; then, based on the CFD method, the mountain wind field is simulated, the wind field information is introduced, and the UAV energy consumption is considered to obtain the trajectory, and the model is verified to be better. In terms of the algorithm, considering the safety and operating efficiency of the trajectory, the traditional A* algorithm is improved to make the solution more convenient and rapid, and a safe, energy-saving and reliable trajectory is obtained, which verifies that the improved algorithm is better.

[0011] The present invention adopts the following technical solutions to solve the technical problems:

[0012] The path length, threats, energy consumption, etc. in UAV operation are comprehensively considered as optimization objectives, and the trajectory planning model is established considering the UAV's own performance constraints, flight environment constraints and other constraints. The wind data is simulated by the CFD method, and the impact of low-altitude wind on trajectory planning is considered. The A* algorithm is improved to improve the safety of trajectory planning and the speed of model calculation, so as to obtain a safe and reliable UAV trajectory and complete the UAV trajectory planning work.

[0013] The beneficial effects of the present invention are as follows:

[0014] First: In view of the current situation that the actual energy consumption of UAVs has not been considered in previous trajectory planning, when establishing the UAV trajectory planning model, the impact of UAV energy consumption is added on the basis of the previous considerations of mountain obstacle threats and other threats, making trajectory planning more energy-efficient.

[0015] Second: In view of the current situation that previous trajectory planning did not consider the impact of wind on trajectory planning, a UAV trajectory correction model under low-altitude wind was established by considering the impact of low-altitude wind on UAV trajectory planning. The wind field was simulated based on the CFD method to make it close to the wind field distribution near the mountainous area, providing a basis for solving the low-altitude wind track correction model.

[0016] Third: In view of the current situation that the safety of the trajectory and the slow calculation speed were not considered when solving the trajectory planning in the past, improvements were made to the shortcomings of the traditional A* algorithm, and more safety factors were considered, so that the algorithm can be "adaptively" adjusted during operation. When the UAV approaches the mountain, it will automatically detour, making the trajectory safer and smoother. At the same time, the search speed is accelerated by constructing the HFM matrix to avoid repeated calculations, jump point search and other methods, making the algorithm more efficient. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] Figure 1A schematic diagram of a no-fly threat according to the present invention;

[0018] Figure 2 It is a schematic diagram of the atmospheric threat of the present invention;

[0019] Figure 3 It is a schematic diagram of the minimum track segment constraint of the present invention;

[0020] Figure 4 It is a schematic diagram of the grid division method of the DEM of the present invention;

[0021] Figure 5 (a) is a track deviation diagram when the speed direction of the UAV of the present invention remains unchanged;

[0022] Figure 5 (b) a trajectory diagram for correcting the speed direction of the UAV of the present invention;

[0023] Figure 6 This is the improved effect of the A* algorithm of the present invention;

[0024] Figure 7 (a) is a node graph in the OPEN set of the present invention;

[0025] Figure 7 (b) is a schematic diagram of a grid expansion mode of the present invention;

[0026] Figure 8 A schematic diagram of a jump point search according to the present invention;

[0027] Fig. 9 The improved A* algorithm of the present invention is a flow chart for solving the UAV trajectory;

[0028] Fig.10 The calculation residual graph of the present invention;

[0029] Fig.11 (a) is the distribution at a wind speed of 3500 m of the present invention;

[0030] Fig.11 (b) is the distribution at a wind speed of 4500 m of the present invention;

[0031] Fig.12 (a) is the track diagram under viewing angle 1 obtained by the original A* algorithm of the present invention;

[0032] Fig.12 (b) is the track diagram under view angle 2 obtained by the original A* algorithm of the present invention;

[0033] Fig.12 (c) is the track diagram under viewing angle 3 obtained by the original A* algorithm of the present invention;

[0034] Fig.13(a) is the track diagram under viewing angle 1 obtained by the IA* algorithm of the present invention;

[0035] Fig.13 (b) is the track diagram under viewing angle 2 obtained by the IA* algorithm of the present invention;

[0036] Fig.13 (c) is the track diagram under viewing angle 3 obtained by the IA* algorithm of the present invention;

[0037] Fig.14 (a) is the track diagram under viewing angle 1 obtained by the IIA* algorithm of the present invention;

[0038] Fig.14 (b) is the track diagram under viewing angle 2 obtained by the IIA* algorithm of the present invention;

[0039] Fig.14 (c) is the track diagram under viewing angle 3 obtained by the IIA* algorithm of the present invention;

[0040] Fig.15 (a) is a comparison of the track diagrams obtained under the viewing angle 1 by different algorithms of the present invention;

[0041] Fig.15 (b) is a comparison of the track diagrams obtained under viewing angle 2 by different algorithms of the present invention;

[0042] Fig.15 (c) is a comparison of the track diagrams obtained under viewing angle 3 by different algorithms of the present invention;

[0043] Fig.16 (a) is the track diagram under view angle 1 obtained by using the IIA* algorithm in the present invention with consideration of energy consumption and wind;

[0044] Fig.16 (b) is the track diagram under view angle 2 obtained by using the IIA* algorithm in the present invention with consideration of energy consumption and wind;

[0045] Fig.16 (c) is the track diagram under viewing angle 3 obtained by using the IIA* algorithm in the present invention with consideration of energy consumption and wind;

[0046] Fig.17 (a) is a comparison between the initial track of the present invention and the track of view 1 considering energy consumption and wind conditions;

[0047] Fig.17 (b) is a comparison between the initial track of the present invention and the track of view 2 considering energy consumption and wind conditions;

[0048] Fig.17 (c) is a comparison between the initial track of the present invention and the track of view 3 considering energy consumption and wind conditions;

[0049] Fig.18 It is a schematic diagram of the process of the present invention. DETAILED DESCRIPTION

[0050] The path length, threats, energy consumption, etc. in UAV operation are comprehensively considered as optimization objectives, and the trajectory planning model is established considering the UAV's own performance constraints, flight environment constraints and other constraints. The wind data is simulated by the CFD method, and the impact of low-altitude wind on trajectory planning is considered. The A* algorithm is improved to improve the safety of trajectory planning and the speed of model calculation, so as to obtain a safe and reliable UAV trajectory and complete the UAV trajectory planning work.

[0051] The following assumptions are made for the UAV trajectory planning model:

[0052] (1) The release location of the UAV and the search and rescue points assigned to each UAV are known and determined, that is, the trajectory planning is carried out based on the UAV task allocation results in Chapter 3;

[0053] (2) The energy consumption of UAV-carried equipment such as camera shooting activities is not taken into account, and only the energy consumption generated by flight activities is taken into account;

[0054] (3) The UAV maintains a consistent speed in each flight path segment;

[0055] (4) Ignore the impact of dynamic obstacles on UAV trajectory planning;

[0056] (5) Each target search point is searched by only one UAV and only once;

[0057] (6) After the UAV completes all its assigned search missions, it returns to the location where the general aviation helicopter released the UAV;

[0058] (7) The drone model is consistent;

[0059] The quality of the UAV's trajectory needs to be evaluated. The key to the trajectory evaluation standard is to construct an objective function. The trajectory evaluation function is used to calculate the suitability of the trajectory, which is an important criterion for judging the quality of the path and the key to guiding the search algorithm to achieve the optimal solution. The evaluation of the trajectory cost also needs to consider various constraints.

[0060] The explanation of the variables in the UAV trajectory planning model is shown in Table 1:

[0061] Table 1 Variables of UAV trajectory planning model

[0062]

[0063]

[0064] Objective Function

[0065] Finding a feasible UAV flight path is a complex optimization problem with multiple objectives. Considering the path length, environmental constraints and collision avoidance, the evaluation function for UAV path planning can be written as follows:

[0066] f=f L +f T +f w +f e (1)

[0067] where f L is the flight path length cost, f T is the terrain cost, i.e. the cost of safe flight over the mountain, f w is the threat cost, f e is the energy consumption cost of the drone.

[0068] The objective function can be expressed as formula (2):

[0069] minf (2)

[0070] (1) Minimum flight path length cost

[0071] The minimum flight path cost is used to measure the impact of the UAV's range on the time required to perform a search and rescue mission. For flight missions, a shorter path length is always preferred over a longer path length, because a shorter path usually reduces energy consumption and can reduce the possibility of encountering unexpected threats. In order to characterize this factor, the path length ratio is used to evaluate the track length, as shown in formula (3):

[0072]

[0073] Where N w is the number of track segments formed by the total track waypoints, including the starting waypoint P0 = (x0, y0, z0), the destination waypoint and intermediate waypoints. Then there is always f L ≥1, smaller values ​​mean shorter tracks.

[0074] (2) Mountain terrain costs

[0075] In order for the drone to fly safely, the planned path should not cross the terrain and must avoid all mountains in the 3D environment. If at least one point in the path is inside a mountain, the solution should be penalized. The cost function of the mountain constraint can be expressed as Equations (4) and (5):

[0076]

[0077]

[0078] Among them Ai,j is a binary value. For the segment between the (i-1)th and i-th waypoints, it is divided into ζ segments, where the jth grid coordinate is represented by (x i,j ,y i,j ,z i,j ). The horizontal position of the kth mountain grid (x i,j ,y i,j ) The number of segments ζ is set to the number of grids passed between the (i-1)th and i-th waypoints.

[0079] (3) Threat Cost

[0080] 1) No-fly threat (cylindrical)

[0081] Cylindrical drone no-fly zones (areas with weak signals when flying drones in mountainous areas) are listed here in a matrix, where the kth row of the matrix has coordinates (x k ,y k ,r k ) represents the central abscissa and ordinate of the kth no-fly zone and the prohibited range, then the no-fly zone matrix can be constructed as (6):

[0082]

[0083] The no-fly threat level is set to T m , let the distance between the drone and the i-th no-fly threat center be dis i , the no-fly threat radius is r i Then we have equation (7) and equation (8):

[0084]

[0085] 2) Atmospheric threats (spherical)

[0086] The harsh atmospheric environment can be abstracted as a sphere, assuming that the range of atmospheric threat is [R Tmin ,R Tmax ], when the distance between the drone and the threat center is greater than R Tmax When the UAV is not affected by the atmospheric threat, the atmospheric threat level is 0. The closer the UAV is to the threat center, the greater the threat level. The atmospheric threat level is set to T w Then we have formula (8):

[0087]

[0088] Generally speaking, at a certain point, the threat may be either a no-fly threat or an atmospheric threat. At this time, the threat cost is the corresponding threat. If both threats appear at the same time, the one with the larger value will be used for calculation.

[0089] Assume B i,j =max(T W ,T m ) then the threat cost f w It can be expressed as formula (9):

[0090]

[0091] (4) Energy consumption cost

[0092] References

[70] The power P from the UAV's i-1th track point to the i-th track point i-1,i The expression of is as follows:

[0093]

[0094] Among them C D is the air resistance coefficient, a1 is the front area of ​​the drone, ep is the empty weight of the drone, and f y is the payload weight carried by the drone, ρ is the density of air, and b1 is the width of the drone. The airspeed v of the drone is defined by the following formula (11):

[0095] v 2 =(v g,i-1,i cosθ i-1,i -WS i-1,i cosκ i-1,i ) 2 +(v g,i-1,i sinθ i-1,i -WS i-1,i sinκ i-1,i ) 2 (11)

[0096] Among them, v g,i-1,i Indicates the ground speed from the i-1th track point to the i-th track point. In windless conditions, the ground speed is the speed of the drone; WS i-1,i represents the wind speed from the i-1th track point to the i-th track point; θ i-1,i represents the ground speed v g,i-1,i The vector angle of i-1,i Indicates the wind direction from the i-1 track point to the i-th track point. i-1 ,y i-1 ,z i-1 )~(x i ,y i ,z i ) Energy consumption E in the flight path i-1,i It is expressed as formula (12):

[0097]

[0098] Assume the total battery capacity of the drone is E max , then the energy consumption ratio is used to evaluate the energy consumption cost as formula (13):

[0099]

[0100] 4.3.3 UAV performance constraints

[0101] Since this article uses a rotorcraft UAV, there are no maximum turning angle constraints and maximum pitch angle constraints. The performance constraints of the UAV in this article mainly include flight altitude constraints, minimum path length constraints, etc.

[0102] (1) Flight altitude constraints

[0103] When the drone flies too low, it is easy to cause damage to the drone itself, which is not conducive to flight. During the flight of the drone, although there is no requirement for hovering to shoot, you should try to capture the scene. When the flight altitude is too high, the shooting equipment cannot meet the shooting accuracy requirements. Therefore, the flight altitude should be kept within a certain range. The minimum flight altitude is H min (During normal flight, the drone should maintain a safe distance from the ground), the maximum flight altitude (i.e. the drone ceiling) is H max , let the flying height of the UAV at the i-th grid point be H i , the corresponding terrain height is h i , the minimum relative altitude of flight is h safe , then the relative flight height of the UAV needs to satisfy equations (14) and (15):

[0104] H min ≤H i ≤H max (14)

[0105] H min =h i +h safe (15)

[0106] (2) Minimum path length

[0107] Before the rescue UAV changes its flight direction or flight posture, there is a minimum inertial distance constraint, that is, the performance of the UAV itself determines that it must fly a certain distance. This shortest distance is defined as the minimum path length, set to l min , the length of the i-th path is l i , then the constraint (16) needs to be satisfied:

[0108] l i ≥l min (16)

[0109] 4.3.4 UAV flight environment constraints

[0110] When planning a path for each drone, some important factors need to be considered, such as the environment of the mission area, the safety of the drone, and the cost of each path. Generally speaking, the mission area environment mainly includes the flight range, terrain conditions, threat level, etc. The target function should consider the necessary environmental factors. The environmental constraints are introduced below.

[0111] (1) Terrain modeling and computational domain grid division

[0112] 1) In this paper, we assume that mountains are the main terrain features and each drone cannot cross them. However, the known data of these mountains, that is, the elevation data, have been obtained. The grid method can be used to divide the studied space and assign corresponding geographic information to each grid.

[0113] In order to prevent data distortion, the data in the X and Y directions are stored in the form of coordinates, where the coordinates are the center coordinates of each grid, and the information in the grid is the height information, that is, the coordinates of the mountain in the Z axis direction.

[0114] 2) For the calculation domain, compare the accuracy of the elevation data and the size of the drone body to determine the grid size. The X- and Y-coordinates of the elevation data point correspond to the X- and Y-coordinates of the center of the grid. The Z-axis coordinate of the center of the grid is determined by its height. The altitude of the drone will not be lower than the sum of the height of the mountain corresponding to its location and the minimum relative height. Suppose the length, width, and height of the drone body are xx1, yy1, and zz1, respectively. Suppose size2 = max(xx1, yy1, zz1), and expand the size of the drone body to a cube grid of size2×size2×size2. The accuracy of the elevation data is set to size1, and the unit grid size to be determined is at least size×size×size, then size is defined as formula (17):

[0115]

[0116] If the size of the unit grid is larger than size×size×size, such as 2size×2size×2size, considering the safety of UAV flight, the elevation data takes the maximum value of the original four adjacent elevation data. At this time, the accuracy of track planning is also reduced. Conversely, when the size of the unit grid is smaller, although the accuracy of track planning is higher, the planning time is longer.

[0117] (2) Flight space constraints

[0118] The goal of UAV path planning is to find a feasible and optimal path to the target area under complex environment and state constraints. i ,y i ,z i )(unit: km) is the coordinate of the path waypoint in the three-dimensional environment. The flight space of the UAV can be expressed as (18):

[0119] {x i ,y i , z i |x min <x i <x max ,y min <y i <y max ,z min <z i <z max} (18)

[0120] where x min ,x max ,y min ,y max ,z min ,z max They respectively represent the flight boundaries of the drone.

[0121] For the convenience of calculation, the longitude and latitude coordinates of the terrain elevation data are converted into rectangular coordinates. Suppose the longitude and latitude coordinates of a point are (xx i ,yy i ,zz i )(zz i The unit is m), and the coordinates after conversion to rectangular coordinates are (x i ,y i , z i ), then we have formula (19):

[0122]

[0123] The horizontal and vertical coordinates of the rectangular coordinate system can be obtained by the above conversion formula, and the radius of the earth is 6371 km. 4.4 Track correction model under the influence of low-altitude wind

[0124] Based on the wind field simulation data of the study area and the initial track of the UAV, the influence of low-altitude wind on the flight of the UAV is introduced, and a track correction model is established according to the velocity vector relationship to correct the initial track. In three-dimensional space, WS represents the wind speed direction in the three-dimensional scene, and v represents the flight speed of the UAV under windless conditions.

[0125] The following discusses the most general case. In the space plane where WS and v are located, under windless conditions, the UAV flies along the direction A1→B1 to obtain the initial flight track. Under the action of WS, the UAV will fly along the direction A1→B2, thus deviating from the original track. In order to make the UAV fly according to the original track, the direction of v is changed, and the actual ground speed v is synthesized according to the velocity vector g ,θ cr For v and v g The angle between the two, ΔL is the track correction, t x is the time taken by the UAV to pass A1→B1 after the track is corrected, where v g is the ground speed v at point A1 gA1 Half of the sum of the vectors of the ground speed at point B1.

[0126] Then the modified model is formula (20):

[0127]

[0128] In order to facilitate the calculation in three-dimensional coordinates, v and WS are decomposed into components in the x direction, y direction, and z direction for calculation, where the components of v in the x direction, y direction, and z direction are v x 、v y and v z , the components of WS in the x-direction, y-direction, and z-direction are WS x , WS y and WS z , then v g The components in the x-direction, y-direction and z-direction are synthesized by the components of v and WS in each direction, so the correction model of the velocity component is obtained as equations (21) and (22):

[0129]

[0130]

[0131] Among them, ΔL x , ΔL y , ΔL z The directions of v gx 、v gy 、v gz The direction is consistent, and the direction of ΔL is the same as v g The key to the traditional A* algorithm is to construct the heuristic cost F(n), which is generally expressed as formula (23):

[0132] F(n)=G(n)+H(n) (23)

[0133] Among them, G(n) is generally set as the objective function value, and H(n) is the heuristic function, which can be set as the Euclidean distance from the current node to the target point. However, the traditional setting method still has shortcomings. It only optimizes the length of the path and ignores the safety of the path. When the object approaches an obstacle, it does not seek the shortest distance but is very close to the obstacle, which increases the operational risk for the aircraft. At the same time, near the more "sharp" track section, the aircraft needs to make "sharp turns" to fly according to the planned track, which also increases the safety risk, and the efficiency of the traditional A* algorithm needs to be improved.

[0134] In view of the shortcomings of the classic A* algorithm, some improvements are proposed. Among them, when only the heuristic cost and algorithm speed are improved, the algorithm is called IA*, and the improved algorithm that adds a track smoothing step on the basis of IA* is called IIA* algorithm.

[0135] For each obstacle in the terrain map, we assign it a weight w, which is the distance from the current node to the closest inaccessible obstacle. The actual cost G(n) is replaced by formula (24):

[0136] G'(n)=G(n) / w (24)

[0137] Then the heuristic cost of the improved method is formula (25):

[0138] F'(n)=G'(n)+H(n) (25)

[0139] Since the terrain map contains altitude information, the allowed altitude should be added to the search strategy to obtain a logically optimal approach.

[0140] The safety factor is introduced to improve the heuristic function, and the angle evaluation function is introduced to improve the heuristic function of the A* algorithm. To meet the needs of this article, it is extended to three-dimensional planning.

[0141] (1) Introducing path safety factor

[0142] In the heuristic function of the traditional A* algorithm, only the Euclidean distance from the extended node to the target node is included. This path has a shorter distance, but is usually close to obstacles, so it is more dangerous and not suitable for drone flight. The improved heuristic function includes the target guidance factor and the path safety factor, as shown in equations (26), (27) and (28):

[0143] H(n)=h goal (n)+h obs (n) (26)

[0144]

[0145] hobs (n) = max(p obst1 ,p obst2 ,...) (28)

[0146] The improved algorithm should not only keep the drone at an appropriate safe distance from obstacles, but also plan feasible paths in narrow terrains. The improved method is to increase the distance of the node to the nearest obstacle in the heuristic function of each feasible node.

[0147] This paper designs the relationship between safety distance and danger level. A nonlinear function is used to convert the safety distance into danger level p obsti The range of the danger level is from 0% to 100% (0% means no collision will occur, and 100% means a collision will occur). The function equation is (29):

[0148]

[0149] The parameter τ is the steepness coefficient, d obsti is the distance to obstacle i. τ varies in the positive range. When the distance is zero, the risk p obsti The steepness factor τ is set to 100%. The larger the steepness factor τ is set, the faster the danger level decreases with distance. Through many simulation experiments, it is appropriate to set the steepness factor τ to 5.

[0150] (2) Improvement of the heuristic cost function by the angle evaluation function

[0151] The selected heuristic function is Euclidean distance, then h1(n) (original H(n)) is formula (30):

[0152]

[0153] The biggest problem with the cost function search of the A* algorithm is that it cannot select the optimal path with the smallest inflection point among many symmetric paths with the same path length. This paper improves the A* algorithm by combining the angle cost estimation method. Assume that the coordinates of node n are (x n ,y n ,z n ), then the coordinates of the previous node are (x n-1 ,y n-1 ,z n-1 ), the coordinates of the next node are (x n+1 ,y n+1 ,z n+1 ). Then the angle evaluation function of node n is h2(n) as formula (31):

[0154]

[0155] Combining the heuristic function and the angle evaluation function, the search heuristic cost of the improved A* algorithm is as follows:

[0156] F(n)=G(n)+k1h1(n)+k2h2(n) (32)

[0157] References [72 ,73], introduced the changing greedy weight coefficient, heuristic memory, and jump point search method to improve the real-time performance of the algorithm and speed up the search speed.

[0158] (1) Introducing the greedy weight coefficient

[0159] In the weighted A* algorithm, the heuristic function H(n) is multiplied by the greedy weight coefficient ε to quickly find a feasible path, and then we get formula (33):

[0160] H'(n)=ε×H(n) (33)

[0161] However, the approach speed of the OPEN set to the target node is not constant, and the search process usually falls into a local optimal situation. The algorithm usually expands many unnecessary nodes, which wastes a lot of time. Therefore, the greedy weight coefficient designed in this paper is not a predetermined fixed value, but is related to the rate of change Δd(n) of the distance between the current OPEN set of the algorithm and the target node. When the algorithm falls into a local optimal situation, the algorithm does not expand to the target node for a long time, so the greedy coefficient can be set to formula (34):

[0162]

[0163] (2) Heuristic function memory speeds up search

[0164] In this paper, a heuristic function memory (HFM) is set to improve the real-time performance of the algorithm. Since the heuristic function of each node must be repeatedly calculated when selecting an expansion node, the computational cost of the A* algorithm is relatively high. This study uses the HFM matrix to store the heuristic function information of each node. When the node is expanded again, the function value can be directly read from the HFM matrix, avoiding repeated calculations and greatly improving the algorithm speed.

[0165] (3) Jump point search

[0166] Network expansion mode: The square in the middle represents the node in the OPEN set that will be expanded into the CLOSE set, and the squares around it represent the nodes that will be added to the OPEN set.

[0167] Jump point search is to find points that can jump between the starting point and the target point, and find the path by searching between the jump points. This is the search strategy of the jump point search algorithm. As shown in the figure, the G1 point can be searched directly from the S1 point without passing through the nodes n2 or n4. The intermediate nodes on the path formed by the jump points (such as S1 and G1) will not be expanded, which also greatly reduces the amount of calculation and storage in the expansion process and speeds up the search.

[0168] After integration, the heuristic cost of the improved algorithm is obtained as formula (35):

[0169]

[0170] The letters in the formula have exactly the same meanings as described above.

[0171] ,Smoothing the path,The smoothing in the literature only considers deleting nodes from a geometric perspective, while here it is applicable to the scenario of this article, and the strategy of deleting nodes is changed.

[0172] Let P = {P0, P1, P2, ..., P Nw-1 ,P Nw} is the set of unsmoothed track points calculated, and the path smoothing algorithm is as follows:

[0173] Input: P = {P0, P1, P2, ..., P Nw-1 ,P Nw};

[0174] Output: the newly set P;

[0175] Step 1: i=0;

[0176] Step 2: Calculate InterSet1, ie P i P i+2 and the intersection of surrounding obstacles and threats; when energy consumption and wind are considered, InterSet2 is calculated, that is, P i P i+2 Is the calculated energy consumption after introducing wind information greater than P i P i+1 and P i+1 P i+2 If yes, then InterSet2 = 1, otherwise InterSet2 = 0;

[0177] Step 3: When energy consumption and wind are not considered, if InterSet1 is not empty, then i=i+2, go to step 2, otherwise, execute step 4; When energy consumption and wind are considered, if InterSet1 is not empty or InterSet2=1, then i=i+2, go to step 2, otherwise, execute step 4;

[0178] Step 4: Delete P from PointSet i+1 To obtain a new set, P = P0, P1, P2, P3, ..., P i ,P i+2 ,...,P Nw}, then go to step 1;

[0179] Step 5: If i=n-2, output the final set P.

[0180] The final point set only includes key points belonging to the route. In this way, the zigzags in the path can be effectively reduced.

[0181] The steps of improving the A* algorithm to solve the track considering energy consumption and wind influence are as follows:

[0182] Step 1: Consider the energy consumption setting G(n), and also set H(n). According to the algorithm improvement steps, G'(n) and H'(n) are obtained, thereby obtaining the expression of F(n), calculating HFM, and using the CFD method to calculate the wind data;

[0183] Step 2: Put the starting point into openlist;

[0184] Step 3: Use the jump point search method to find the points around the node that can be reached and skip the points in the closelist. This point is used as the parent point of these points, and the point is deleted from the openlist and added to the closelist.

[0185] Step 4: Calculate the ground speed according to the track correction model, correct the time of passing the point, and then calculate the energy consumption of the track segment according to the energy consumption calculation model, and then calculate the F(n) of the point.

[0186] Step 5: Check whether F(n) is the minimum value in the openlist and satisfies the constraints. If not, add the point to the closelist, otherwise proceed to the next step.

[0187] Step 6: If the node is found successfully, delete it from the openlist table and add it to the closelist;

[0188] Step 7: Determine whether the point is the target node. If not, return to step 3. If yes, proceed to the next step.

[0189] Step 8: Smooth the track according to the track smoothing method in the algorithm improvement to obtain the final track.

[0190] Two points located near the observation point in Xiaojin County, Sichuan Province are selected for trajectory planning example analysis. The two points are (102.45013, 31.20414, 2733) and (102.50767, 31.43384, 3015.00000), where the first item of the coordinates is longitude, the second item is latitude, and the third item is altitude (m).

[0191] The raw data obtained include drone performance data, wind data of the studied area, and elevation data (with an accuracy of 20m). The drone performance data used mainly include energy consumption data under various flight conditions and the total battery capacity of the drone.

[0192] In addition, in order to verify the correctness of the wind field simulation, it is also necessary to obtain the wind observation data of the Xiaojin County Observation Station as shown in Table 2:

[0193] Table 2 Hourly wind observation data of Xiaojin County meteorological station

[0194]

[0195] 2 Wind field simulation

[0196] 1) Build the terrain

[0197] Before simulating the wind farm, the terrain needs to be physically constructed. The GH plug-in in Rhino is used to import elevation data, build a terrain solid model, and export the model to stl format. The ICEM plug-in in ANSYS is used to process the terrain.

[0198] 2) Terrain topology and computational domain

[0199] The terrain used in the current numerical simulation is the elevation data product from Google Earth. The accuracy is 20m. The longitude range of the simulated terrain is 102.1505°~102.5494°, the latitude range is 30.9001°~31.4995°, the maximum altitude is 5.020km, and the minimum altitude is 2.175km. The grid size is 20m×20m×20m, and the calculation domain is based on the horizontal plane where the highest point is located and increases by 5km upward.

[0200] 3) Observational data

[0201] The longitude of the Xiaojin County Observation Station is 102°21′, the latitude is 31°00′, and the altitude is 2369.2m. The wind data are obtained from the China Meteorological Data Sharing Center.

[0202] By statistically analyzing the hourly data recorded by the Xiaojin County Meteorological Station, it was found that the average maximum wind speed observed at this station was 4.25m / s, and the average wind force was 1.46, no higher than level three, so the scene is suitable for drone flight.

[0203] 4) Wind field simulation

[0204] Using the method described above, the turbulence model is selected in the Fluent plug-in, and the boundary conditions are set. The inlet velocity is adjusted continuously to make the wind speed near the Xiaojin County Meteorological Observatory close to the average maximum wind speed. The wind field data at the height of the meteorological observatory is extracted to obtain Table 3:

[0205] Table 3 Wind speed information table of some grid nodes at Xiaojin County meteorological station height (section data)

[0206]

[0207] It can be seen from the above table that when the longitude and latitude are 102.3502° and 31.0017° respectively, that is, at the location of the Xiaojin County Observation Station, the wind speed obtained is 4.2521m / s, which is closest to the maximum wind speed of 4.25m / s. Therefore, it is more appropriate to adjust the inlet speed to 5.3m / s, and the residual graph is obtained.

[0208] According to the residual graph, when Fluent calculates the wind field, the residual graph is obtained by iterating 300 steps. After 300 iterations, the residuals of all parameters are less than 10e-3, which means that convergence has been achieved. The section velocity graphs at 3500m and 4500m are given as shown in the figure. Due to the large amount of data, not all data are listed here.

[0209] 3 Simulation parameter settings

[0210] According to the data acquisition, the relevant parameters of the simulation are set as shown in Table 4:

[0211] Table 4 Simulation parameters

[0212]

[0213]

[0214] 4 Solution results and analysis

[0215] (1) Comparison of algorithm advantages and disadvantages

[0216] Since we want to compare the advantages and disadvantages of the algorithms, the energy consumption and wind factors of the drone are not considered, and the objective function does not include energy consumption costs.

[0217] First, the traditional A* algorithm is used to calculate the UAV track map. It is not difficult to see that the track obtained before the improved A* algorithm is more "sharp" in shape, and most of the places are very close to mountains and threats, which obviously increases the risk of UAV operation.

[0218] The IA* algorithm is used to calculate and obtain the UAV track. After improving the heuristic cost, the track obtained will try to stay away from the "mountain", no-fly threats and atmospheric threats, avoid collisions with mountains and other threats, and greatly improve the safety of the track. Although the track is safe, it is obviously longer than the track obtained by the unimproved algorithm. At the same time, there is still room for optimization of the track, that is, smoothing the track to make it shorter and save the UAV's power.

[0219] The track is calculated using the IIA* algorithm. As can be seen from the figure, the smoothed graph looks smoother subjectively, and the path is shorter and has fewer "bumps".

[0220] The tracks obtained by the IIA* and IA* algorithms are compared with those obtained by the unimproved A* algorithm and the RRT algorithm.

[0221] Intuitively, it can be seen that the track obtained by the IIA* algorithm is the safest and can avoid mountain obstacles as much as possible. Objectively, in order to quantitatively analyze the advantages and disadvantages of the algorithm optimization degree, the table of the path length, path points and running time calculated by the IIA* algorithm, IA* algorithm, A* algorithm and RRT algorithm is shown in Table 5, where the path change rate of IIA* relative to other algorithms is Path point change rate and the running time rate of change Calculated by formula (36), the positive and negative signs represent rising and falling respectively.

[0222]

[0223] Table 5 Comparison of calculation paths and running time results

[0224]

[0225] As can be seen from Table 4.5, the average path length of the IIA* algorithm increased by 6.04%, decreased by 2.70%, and decreased by 0.88% compared with the IA* algorithm, A* algorithm, and RRT algorithm, respectively. This is because the IA* algorithm takes more safety factors into consideration. In order to make the flight of the UAV safer, the UAV makes adaptive adjustments when it wants to approach the mountain and tries to "stay away" from the mountain, adding some paths. After path smoothing (IIA*), some unnecessary paths are reduced; the change in the number of path points of the IIA* algorithm is reduced by 74.58%, 72.22%, and 71.15% respectively compared with the other three algorithms, indicating that the IIA* algorithm greatly reduces the complexity of the UAV path; the running time of the IIA* algorithm is increased by 10.31%, decreased by 20.20%, and decreased by 15.46% respectively compared with the other three algorithms. This is because the IIA* algorithm adds the calculation of path smoothing, which increases the time compared with IA*, but is better than other unimproved algorithms.

[0226] (2) IIA* algorithm solves the UAV path considering wind and energy consumption

[0227] In order to illustrate the difference between the models, the IIA* algorithm is used to compare the initial track (the track needs to be calculated first and then the wind information is added to calculate the energy consumption) with the track solution considering wind and energy consumption. The wind field information at some nodes is extracted using the CFD method mentioned above as shown in Table 6:

[0228] Table 6 Wind field information at some nodes

[0229]

[0230] When considering energy consumption, the cost function contains energy consumption cost when planning the UAV trajectory, thus obtaining the solved UAV trajectory. After smoothing, the trajectory under the condition of considering energy consumption and wind still appears "tortuous" in some track sections compared with the trajectory obtained without considering energy consumption and wind. This is because when smoothing, not only the shorter trajectory is sought in geometric space, but also some points are deleted to consider whether the energy consumption will be reduced. The UAV trajectory under energy consumption and wind is compared with the trajectory without considering energy consumption and wind.

[0231] In order to quantitatively analyze the advantages and disadvantages of the model before and after improvement, the path length, energy consumption value and its corresponding change rate are compared, as shown in Table 7. The energy consumption of the initial track is calculated based on the track and wind data obtained, according to equations (12)-(14) and the track correction model, and its path change rate is and energy consumption change rate The calculation method of is as follows:

[0232]

[0233] Table 7 Comparison of energy consumption and wind conditions under the solution results of IIA* algorithm

[0234]

[0235] It can be seen from Table 7 that when the energy consumption is taken into account in solving the model, the path obtained will increase, but the energy consumption value will decrease, energy will be saved, the path length will increase by 2.94%, and the energy consumption of the path will decrease by 5.76%. When the wind factor is considered, the effectiveness of the model is verified from the perspective of energy consumption.

Claims

1. A trajectory planning method considering low-altitude wind and UAV energy consumption, characterized in that: The method comprises the following steps: 1) Obtain the performance data of the UAV, including the maximum range, minimum flight path, ceiling, minimum relative altitude and power when flying at a certain speed; the terrain data of the studied area and divide it into grids, giving each grid corresponding geographical information; wind data of the observation point; 2) Establish an objective function for initial trajectory planning taking into account the energy consumption cost of the UAV; 3) Establish constraints to constrain the objective function; 4) Introducing the influence of low-altitude wind on UAV flight, establishing a track correction model, and correcting the objective function obtained in step 3); wherein the correction model is expressed as: Among them, t x v is the time taken by the UAV to pass A1→B1 after the track is corrected; g is the ground speed of the drone; is the path length from point A1 to point B1; ΔL is the track correction; θ cr is the airspeed v of the drone and the ground speed v of the drone g WS is the wind speed direction; is the ground speed vector of the UAV; is the ground speed vector of the drone at point A1; is the ground speed vector of the UAV at point B1; 5) The improved A-Star algorithm is used to solve the corrected objective function in step 4) to obtain the shortest path of the UAV.

2. A method for trajectory planning considering low-altitude wind and energy consumption of unmanned aerial vehicles as claimed in claim 1, characterized in that: In step 2), the evaluation function of the UAV trajectory planning is first established as: f=f L +f T +f w +f e ; Then the objective function is: minf; f L is the flight path length cost, specifically: Where N w is the number of track segments formed by the track points of the total track, the starting track point P0 = (x0, y0, z0), the destination track point x i ,y i , z i is the track point coordinate, x i-1 ,y i-1 , z i-1 is the coordinates of the previous track point; f T is the terrain cost, i.e. the cost of safe flight over the mountain, specifically: Among them A i,j is a binary value. For each segment between the (i-1)th and i-th track points, it is divided into ζ segments, represented as (x i,j ,y i,j ,z i,j ), z i,j is the horizontal position of the drone (x i,j ,y i,j )'s return altitude, is the horizontal position of the kth grid mountain (x i,j ,y i,j ), the segment number ζ is set to the number of grids passed between the (i-1)th and i-th track points; f w is the threat cost. Assume that the threat level of the drone passing through the jth grid from track point i-1 to i is B i,j but: B i,j =max(T W ,T m ); Where T W is the atmospheric threat level, T m is the no-fly threat level, max means the maximum value, T W , T m The expression is as follows: where dis i is the distance between the UAV and the center of the i-th no-fly threat; r i is the threat radius of the i-th no-fly threat center; L UAV is the distance from the drone to the threat center; R Tmax is the maximum threat radius; From the above, we know that the threat cost f w for: f e is the energy consumption cost of the drone, specifically: Where E max is the total battery capacity of the drone; E i-1,i For drones in (x i-1 ,y i-1 ,z i-1 )~(x i ,y i ,z i ) Energy consumption within the track segment; E i-1,i The specific expression is: where v g,i-1,i represents the ground speed of the UAV from the i-1th track point to the i-th track point. In the absence of wind, the ground speed is the airspeed of the UAV; t i-1,i represents the flight time from the i-1th track point to the i-th track point; P i-1,i represents the power from the i-1th track point to the i-th track point; l i represents the length of the drone’s path in segment i; P i-1,i The expression is: Among them C D is the air resistance coefficient, a1 is the front area of ​​the drone, ep is the empty weight of the drone, and f y The weight of the payload carried by the drone, ρ is the density of the air, b1 is the width of the drone, and v is the airspeed of the drone; The airspeed v of the drone is defined as: v 2 =(v g,i-1,i cosθ i-1,i -WS i-1,i cosκ i-1,i ) 2 +(v g,i-1,i sinθ i-1,i -WS i-1,i sinκ i-1,i ) 2 ; Among them, WS i-1,i represents the wind speed from the i-1th track point to the i-th track point; θ i-1,i Indicates ground speed v g,i-1,i The vector angle of i-1,i Indicates the wind direction from the i-1 track point to the i-th track point.

3. A method for trajectory planning considering low-altitude wind and energy consumption of unmanned aerial vehicles as claimed in claim 2, characterized in that: The constraints in step 3) include flight height constraint and minimum path length constraint, where the flight height H i The constraints are as follows: H min ≤H i ≤H max , H min =h i +h safe ; H min is the minimum flight altitude, H max is the maximum flight altitude, and the flight altitude of the UAV at the i-th grid point is H i , the corresponding terrain height is h i , the minimum relative altitude of flight is h safe ; Before the drone changes its flight direction or flight posture, there is a minimum inertial distance constraint, that is, the drone's own performance determines that it must fly a certain distance. This shortest distance is defined as the minimum path length, set to l min ; l i is the length of the i-th path; The minimum path length constraint is as follows: l i ≥l min ； Compare the elevation data accuracy size1 and the size of the UAV body size2×size2×size2 to determine the grid size and improve the accuracy of track planning; Assume that the length, width and height of the drone body are xx1, yy1 and zz1 respectively, and size2 = max(xx1, yy1, zz1). The size of the drone body is expanded to a cube grid of size2×size2×size2. Assume that the unit grid size to be determined is at least size×size×size, then define the size function as: size2=max(xx1,yy1,zz1) The flight space constraints of the UAV (x i ,y i , z i ), expressed as: {x i ,y i ,z i |x min <x i <x max ,y min <y i <y max ,z min <z i <z max }; where x min ,x max ,y min ,y max ,z min ,z max They respectively represent the limits of drone flight; To facilitate calculation, the longitude and latitude coordinates of the terrain elevation data are converted into rectangular coordinates. Suppose the longitude and latitude coordinates of a certain point are (xx i ,yy i ,zz i ), the coordinates after conversion to rectangular coordinates are (x i ,y i , z i ), 4. A method for trajectory planning considering low-altitude wind and energy consumption of unmanned aerial vehicles as claimed in claim 1, characterized in that: In order to facilitate the calculation in three-dimensional coordinates, v and WS are decomposed into components in the x direction, y direction, and z direction for calculation, where the components of v in the x direction, y direction, and z direction are v x 、v y and v z , the components of WS in the x-direction, y-direction, and z-direction are WS x , WS y and WS z , then v g The components in the x-direction, y-direction, and z-direction respectively The modified model of velocity component is obtained by synthesizing the components of v and WS in each direction: Where, ΔL x , ΔL y , ΔL z The actual track correction in the x, y, and z directions, respectively, and ΔL x , ΔL y , ΔL z The directions of v gx 、v gy 、v gz The direction is consistent, and the direction of ΔL is the same as v g Consistent.

5. A method for trajectory planning considering low-altitude wind and energy consumption of unmanned aerial vehicles as claimed in claim 4, characterized in that: Step 5) The heuristic cost of the improved A-Star algorithm is expressed as: Where G(n) is the actual cost from the initial node to the target node n; w is the weight of each obstacle; Δd(n) is the rate of change of the distance between the OPEN set and the target node n; h goal (n) is the target guiding factor; h obs (n) is the path safety factor; k1 and k2 are weight factors; h2(n) is the angle evaluation function of the target node n; h obs (n)=max(p obst1 ,p obst2 ,...); where h goal The coordinates of the nodes in (n) are (x, y, z); x goal ,y goal 、z goal is the target guidance value in the x, y, and z directions, h obs (n) p obsti is the degree of danger, and the coordinates of the nodes in h2(n) are (x n ,y n ,z n ), the coordinates of the previous node are (x n-1 ,y n-1 ,z n-1 ), the coordinates of the next node are (x n+1 ,y n+1 ,z n+1 ), τ is the steepness coefficient, d obsti is the distance to obstacle i.

6. A method for trajectory planning considering low-altitude wind and energy consumption of unmanned aerial vehicles as claimed in claim 5, characterized in that: Solving the objective function using the improved A-Star algorithm includes the following steps: Step 1: Let minf = F(n), calculate HFM, HFM is the set heuristic function memory, the heuristic function information of each node is stored through the HFM matrix, when the node is expanded, the function value is directly read from the HFM matrix, avoiding repeated calculations, and the wind data is calculated using the CFD fluid mechanics method; Step 2: Add the initial node to the openlist status table; Step 3: Use the jump point search method to skip the points in the closelist state table to find the points that can be reached around the node. This point is used as the parent point of these points. The point is deleted from the openlist state table and added to the closelist state table. Step 4: Calculate the ground speed according to the track correction model, correct the passing time, calculate the energy consumption of the track segment according to the energy consumption cost, and then calculate F(n) at the point; Step 5: Check whether F(n) is the minimum value in the openlist state table and satisfies the constraints. If not, add the point to the closelist state table, otherwise proceed to the next step. Step 6: If the node is found successfully, delete the node from the openlist state table and add it to the closelist state table; Step 7: Determine whether the node is the target node. If not, return to step 3. If yes, proceed to the next step. Step 8: Smooth the calculated track according to the track smoothing method, delete the zigzag track, and obtain the final track.

7. A method for trajectory planning considering low-altitude wind and energy consumption of unmanned aerial vehicles as claimed in claim 6, characterized in that: The specific steps of the track smoothing method in step 8) are: Let P = {P0, P1, P2, ..., P Nw-1 ,P Nw } is the set of unsmoothed track points, Step 1: i = 0; P i represents the i-th unsmoothed track point; Step 2: Calculate InterSet1, InterSet1 represents P i P i+2 The set of intersections with surrounding obstacles and threats; When energy consumption and wind are considered, InterSet2 is calculated. InterSet2 represents P i P i+2 Is the calculated energy consumption after introducing wind information greater than P i P i+1 and P i+1 P i+2 The total energy consumption of the set, if yes, then InterSet2 = 1, otherwise InterSet2 = 0; Step 3: When energy consumption and wind are not considered, if InterSet1 is not empty, then i=i+2, go to step 2, otherwise, execute step 4; when energy consumption and wind are considered, if InterSet1 is not empty or InterSet2=1, then i=i+2, go to step 2, otherwise, execute step 4; Step 4: From P = {P0, P1, P2, ..., P Nw-1 ,P Nw }Delete P i+1 To obtain a new set, P = {P0, P1, P2, P3, ..., P i ,P i+2 ,...,P Nw }, then go to step 1; Step 5: If i=n-2, output the final set P.

Images