Multi-unmanned aerial vehicle task allocation method based on dynamic route set description algorithm

Abstract

The invention relates to a multi-unmanned aerial vehicle task allocation method based on a dynamic route set description algorithm, and the method comprises the steps: constructing an initial task arrangement model, setting a space partitioning model of the dynamic route set description algorithm, completing the task allocation integrated with a route optimization method, screening a most preferable task allocation set through a decision model, and completing the task allocation. And selecting the optimal number of unmanned aerial vehicles and calculating the voyage. According to the method, the flight task allocation problems of optimal unmanned aerial vehicle number judgment, flight task segmentation, flight safety and the like are considered. Example results of the method show that flight safety and task division fair multi-unmanned-aerial-vehicle task allocation can be guaranteed under the multi-unmanned-aerial-vehicle and multi-task conditions through the dynamic route set description algorithm, and the method can provide an effective cooperation mode for civil unmanned aerial vehicles and improve task execution efficiency and economic benefits at the same time.

Description

Technical Field

[0001] The present invention belongs to the field of UAV mission allocation, and involves designing an optimal allocation plan using constraints such as the economic performance and rationality of mission allocation. Specifically, it is a multi-UAV mission allocation method with a dynamic route set characterization algorithm. Background Art

[0002] So far, in the research on mission allocation and route planning of multi-UAV emergency monitoring in complex environments at home and abroad, a unified design model has not been formed yet, and there is no solution to optimize the rapid integration of emergency tasks. Only the emergency mission allocation is most significantly affected by the allocation timeliness. The low mission allocation efficiency will limit the UAV monitoring time, thus affecting the accuracy of rescue information, making the rapid emergency rescue more difficult. How to efficiently allocate mission points and form flight routes to achieve rapid and safe rapid planning of emergency mission flight routes, thereby improving the UAV monitoring efficiency, still requires in-depth research and exploration. The solution to this problem has extremely important practical significance for broadening the application of civilian UAVs in both military and military application effectiveness. Summary of the Invention

[0003] (1) Construct an initial task arrangement model, and the mathematical description is as follows:

[0004] z = xrd + 1i * yrd (1)

[0005] xrd = [xrdi], i = 1, 2,..n

[0006] yrd = [yrdi], i = 1, 2,..n

[0007] In formula (1), z, i, xrd, and yrd represent the model values of the initial task points to be allocated, complex identifiers, abscissa vectors, and ordinate vectors, with the unit of kilometer. xrdi, yrdi, and n represent the abscissa and ordinate of the i-th task point in the planned airspace and the number of task points, with the unit of kilometer and piece respectively; "[]" represents the matrix symbol;

[0008] (2) Establish a spatial partition model of the dynamic route set characterization algorithm:

[0009] First, use formulas (2) to (3) to obtain the abscissa and ordinate arrays mcz of the clustering center points of the task points:

[0010] [~, cls] = kmeans_func([xrd; yrd]' ncl] (2)

[0011] mcz = mean(cls) (3)

[0012] In Formulas (2) to (3), kmeans_func, ncl, cls, "~", and "′" are the function for solving the cluster centers, the number of its classifications, the planar coordinates of the obtained cluster center points, the symbol for not outputting the number of categories, and the transpose symbol respectively; mcz and mean are the horizontal and vertical coordinate arrays of the cluster center points and the function for calculating the average value respectively;

[0013] Secondly, use Formulas (4) to (7) to obtain the spatial partitioning model cm corresponding to the dynamic route set characterization algorithm:

[0014] r2 = [rand(c4, 1)*c1, rand(c4, 1)*c2, rand(c4, 1)*c3] (4)

[0015]

[0016] p0 = [x0, y0]

[0017] pe = [xe, ye]

[0018] In Formulas (4) to (5), r2 and rand are the radius set corresponding to the dynamic route chord set and the function for generating random numbers respectively; c1, c2, c3, and c4 are the radius dynamic adjustment coefficients corresponding to the dynamic route chords 1, 2, 3, and 4 respectively; x, y, cm(x, y), cu(x, y), and cd(x, y) are the planar coordinates, the value of the spatial partitioning model of the dynamic route set characterization algorithm at the coordinate (x, y), and the sub-model values of its upper left triangle and lower right triangle regions respectively; in Formulas (6) to (7), p0, x0, y0, pe, xe, and ye are the lower left corner point of the planning area and its abscissa and ordinate, the upper right corner point and its abscissa and ordinate respectively; uxi and uyi are the abscissa and ordinate of the curve part of the dividing line in the upper left triangle region, with the unit of km; fliplr and end are the function for inverting the array elements from left to right and the symbol for taking the last element in the array respectively; nu is the number of dividing regions in the upper left triangle region; dxi and dyi are the abscissa and ordinate of the curve part of the dividing line in the lower right triangle region, with the unit of km; nd is the number of dividing regions in the lower right triangle region;

[0019] Finally, use Formulas (8) to (9) to divide the task sets in each characterization region;

[0020] [in, on] = inpolygon(rw, cm(x, y)) (8)

[0021] frw = [frwu, frwd] = rw(in|on) (9)

[0022] In Formulas (8) to (9), rw, inpolygon, in, and on are respectively the task point coordinate set, the function of the point set included in the selected polygon, and the point sets depicting the interior and boundary of the self-region selected by it. "|" represents the "or" operator in relational operations. frw, frwu, and frwd respectively represent the divided task set, its upper left triangular subset, and its lower right triangular subset;

[0023] (3) Complete the task assignment incorporating the route optimization method:

[0024] Use Formulas (10) to (12) to obtain the adaptive generation route set frwp of the task point set integrating the nearest neighbor sorting algorithm:

[0025] frwdp(i) = fitsortfun([p0, frwu(i), pe]) (10)

[0026] frwdp(i) = fitsortfun([p0, frwu(end - i + 1), pe]) (11)

[0027] frwp = [frwup(i), frwdp(i)], i = 1, 2,... nu + nd (12)

[0028] In Formulas (10) to (12), frwp, frwup, and frwdp are respectively the adaptive generation route set of the task point set integrating the nearest neighbor sorting algorithm, its upper triangular subset, and its lower triangular subset. fitsortfun is the function for finding the integrated nearest neighbor sorting;

[0029] (4) Use the decision model to screen the optimal task assignment set:

[0030] Use Formulas (13) to (16) to obtain the optimal task assignment and route set zhp:

[0031] nb = numel(frwp) (13)

[0032] nba = std(nb) (14)

[0033] [mn, mnc, mnl] = sort(nbz) (15)

[0034] zhp = frwp(mnc(1), mnl(1)) (16)

[0035] In formulas (13) to (16), numel and nb are the function for counting the number of array elements and the number of task points characterizing the partition obtained by it respectively, std and nbz are the function for calculating the standard deviation of the array and the standard deviation of the number of task points within the partition respectively; sort, mn, mnc, and mnl are the function for sorting array elements, the task point route sequence after sorting, the row number and column number corresponding in the original sequence respectively, and zhp is the optimal task allocation and route set.

[0036] (5) Select the optimal number of UAVs and calculate the flight range:

[0037] First, use formula (17) to obtain the number of necessary UAVs selected:

[0038] flyn = flyn + 1, length(zhp) > 0 (17)

[0039] In formula (17), flyn and length are the number of necessary UAVs and the function for calculating the maximum dimension number of the array respectively;

[0040] Secondly, use formulas (18) to (20) to obtain the flight ranges of the sub-task points included in the UAVs corresponding to each frame:

[0041] tp = dist(zhp(j), zhp(j + 1)]), j = 1, 2,... np (18)

[0042] Dt = [Dt tp] (19)

[0043] Dl(i) = sum(Dt), i = 1, 2,... nu + nd (20)

[0044] In formulas (18) to (20), np, dist, tp, and Dt are the number of track points within one route, the function for calculating the distance between two points, the distance between track points j and j + 1, and the flight range of each UAV respectively, sum and Dl are the function for calculating the sum of elements within the array and the sum of the flight ranges of all UAVs respectively; This method realizes the fine optimization of flight task allocation and provides an effective method for improving the task completion efficiency of UAVs.

[0045] Compared with the prior art, the advantages of the present invention are as follows:

[0046] a. It realizes the multi-dimensional dynamic allocation that meets the UAV monitoring tasks. By combining the dynamic route set characterization algorithm and the dynamic adjustment coefficient, taking into account the task space distribution characteristics, it provides an effective method for improving the pertinence, allocation efficiency, and monitoring timeliness of multi-UAV and multi-task.

[0047] b. Fine optimization of the self - decision monitoring routes of multiple UAVs using the integrated route optimization method. By using the integrated design scheme of task allocation, route generation, and optimization, and combining the fine - optimization algorithm for the balance of UAV tasks, the flight routes are finely optimized, realizing the characteristics of self - decision multi - UAV collaborative emergency monitoring routes and comprehensive route optimization. The integration of task allocation and route optimization is achieved in a complex dynamic emergency environment, which plays an important role in quickly and efficiently completing the emergency monitoring task. Detailed implementation mode

[0048] The technical solution of the present invention will be further described below in conjunction with embodiments. In the embodiment, 5 UAVs fly. The coordinates of the lower - left corner of the planned airspace are [0 km, 0 km, 20 m], and the coordinates of the upper - right corner are [20 km, 20 km, 50 m]. There are 30 task points at (16 km, 2 km), (6 km, 2 km), (11 km, 11 km), (3 km, 1 km), (6 km, 7 km), (17 km, 4 km), (6 km, 16 km), (17 km, 7 km), (5 km, 11 km), (19 km, 4 km), (7 km, 13 km), (4 km, 6 km), (6 km, 14 km), (13 km, 14 km), (10 km, 15 km), (8 km, 10 km), (17 km, 2 km), (12 km, 5 km), (11 km, 19 km), (12 km, 9 km); The performance parameters of the UAVs are: the maximum flight time is 35 hours, the minimum cruising speed is 300 km / h, and the range is greater than 3000 km; The coordinates of the starting point and the target point are (1 km, 1 km, 50 m) and (20 km, 20 km, 50 m) respectively.

[0049] (1) Construct an initial task arrangement model, and the mathematical description is as follows:

[0050] z = xrd + 1i * yrd (1)

[0051] xrd = [xrdi], i = 1, 2,..30

[0052] yrd = [yrdi], i = 1, 2,..30

[0053] In formula (1), z, i, xrd, and yrd represent the model values of the initial task points to be allocated, complex identifiers, abscissa vectors, and ordinate vectors, with the unit of kilometers. xrdi, yrdi, and n represent the abscissa and ordinate of the i - th task point in the planned airspace and the number of task points respectively, with the unit of kilometers and pieces; "[]" represents the matrix symbol.

[0054] (2) Establish a spatial partition model for the dynamic route set characterization algorithm:

[0055] First, use formulas (2) to (3) to obtain the horizontal and vertical coordinate arrays mcz of the clustering center points of the task points:

[0056] [~, cls] = kmeans_func([xrd; yrd]', 2) (2)

[0057] mcz = mean(cls) (3)

[0058] In formulas (2) to (3), kmeans_func, ncl, cls, "~", and "'" are respectively the function for solving the clustering center, its number of classifications, the planar coordinates of the obtained clustering center points, the symbol for not outputting the number of categories, and the transpose symbol; mcz and mean are respectively the horizontal and vertical coordinate arrays of the clustering center points and the function for calculating the average value;

[0059] Secondly, use formulas (4) to (7) to obtain the spatial partition model cm corresponding to the dynamic route set characterization algorithm:

[0060] r2 = [rand(50, 1)*5, rand(50, 1)*7.8, rand(50, 1)*8.1] (4)

[0061]

[0062] p0 = [0, 0]

[0063] pe = [20, 20]

[0064] In formulas (4) to (5), r2 and rand are respectively the radius set corresponding to the dynamic route chord set and the function for generating random numbers, and c1, c2, c3, and c4 are respectively the radius dynamic adjustment coefficients corresponding to the dynamic route chords 1, 2, 3, and 4; x, y, cm(x, y), cu(x, y), and cd(x, y) are respectively the planar coordinates, the value of the spatial partition model of the dynamic route set characterization algorithm at the coordinate (x, y), and the sub-model values of its upper left triangle and lower right triangle regions; in formulas (6) to (7), p0, x0, y0, pe, xe, and ye are respectively the lower left corner point of the planning area and its abscissa and ordinate, the upper right corner point and its abscissa and ordinate; uxi and uyi are respectively the abscissa and ordinate of the curve part of the dividing line in the upper left triangle region, with the unit of km; fliplr and end are respectively the function for reversing the elements of the array from left to right and the symbol for taking the last element in the array, nu is the number of dividing regions in the upper left triangle region; dxi and dyi are respectively the abscissa and ordinate of the curve part of the dividing line in the lower right triangle region, with the unit of km; nd is the number of dividing regions in the lower right triangle region;

[0065] Finally, use formulas (8) to (9) to divide the task sets in each characterization region;

[0066] [in, on] = inpolygon(rw, cm(x, y)) (8)

[0067] frw = [frwu, frwd] = rw(in|on) (9)

[0068] In formulas (8) to (9), rw, inpolygon, in, and on are respectively the task point coordinate set, the function for selecting the point set included in the selected polygon, and the point sets depicting the interior and boundary of the self-region selected by it. "|" represents the "or" operator in relational operations. frw, frwu, and frwd respectively represent the task set after segmentation and its upper left triangular subset and lower right triangular subset;

[0069] (3) Complete the task assignment incorporating the route optimization method:

[0070] Use formulas (10) to (12) to obtain the adaptive generation route set frwp of the task point set integrating the nearest neighbor sorting algorithm:

[0071] frwup(i) = fitsortfun([0, frwu(i), 20]) (10)

[0072] frwdp(i) = fitsortfun([0, frwu(end - i + 1), 20]) (11)

[0073] frwp = [frwup(i), frwdp(i)]frwp = [frwup(i), frwdp(i)], i = 1, 2,... 6 (12)

[0074] In formulas (10) to (12), frwp, frwup, and frwdp are respectively the adaptive generation route set of the task point set integrating the nearest neighbor sorting algorithm and its upper triangular and lower triangular subsets. fitsortfun is the function for finding the integrated nearest neighbor sorting;

[0075] (4) Use the decision-making model to screen the optimal task assignment set:

[0076] Use formulas (13) to (16) to obtain the optimal task assignment and route set zhp:

[0077] nb = numel(frwp) (13)

[0078] nbz = std(nb) (14)

[0079] [mn, mnc, mnl] = sort(nbz) (15)

[0080] zhp = frwp(mnc(1), mnl(1)) (16)

[0081] In formulas (13) to (16), numel and nb are the function for counting the number of array elements and the number of partition task points counted by it respectively, std and nbz are the function for calculating the standard deviation of the array and the standard deviation of the number of task points within the partition respectively; sort, mn, mnc, and mnl are the function for sorting array elements, the task point route sequence after sorting, the row number and column number corresponding in the original sequence respectively, and zhp is the optimal task allocation and route set.

[0082] (5) Select the optimal number of UAVs and calculate the flight range:

[0083] First, use formula (17) to obtain the number of necessary UAVs selected:

[0084] flyn = flyn + 1, length(zhp) > 0 (17) In formula (17), flyn and length are the number of necessary UAVs and the function for calculating the maximum dimension number of the array respectively;

[0085] Secondly, use formulas (18) to (20) to obtain the flight ranges of the UAVs corresponding to each UAV including sub-task points:

[0086] tp = dist(zhp(j), tzhp(j + 1)), j = 1, 2,... 11 (18)

[0087] Dt = [Dt tp] (19)

[0088] Dl(i) = sum(Dt), i = 1, 2,... 6 (20)

[0089] In formulas (18) to (20), np, dist, tp, and Dt are the number of track points within a route, the function for calculating the distance between two points, the track point j, the distance between track points j + 1, and the flight range of each UAV respectively, sum and Dl are the function for calculating the sum of elements within the array and the sum of the flight ranges of all UAVs respectively.

[0090] The experimental results are as follows: The maximum flight range of the allocated UAV is 510.85 kilometers, the planning time is 4.09 seconds, and the mean square error of the tasks allocated to the UAVs is only 0.82. Using the dynamic route set characterization algorithm in this method can also quickly obtain the UAV task allocation scheme, and can also plan the flight route with the shortest total distance, and can achieve the rapid and efficient completion of the UAV cooperative flight task, providing technical support for emergency event response.

Claims

1. A multi-UAV task allocation method for a dynamic route set characterization algorithm, characterized by using The following implementation steps: (1) Construct an initial task arrangement model, the mathematical description is as follows: z = xrd + 1i * yrd (1) xrd = [xrdi], i = 1, 2,..n yrd = [yrdi], i = 1, 2,..n In formula (1), z, i, xrd, and yrd represent the model values of the initial task points to be allocated, complex identifiers, abscissa vectors, and ordinate vectors, with the unit of kilometers. xrdi, yrdi, and n represent the abscissa and ordinate of the i-th task point in the planned airspace and the number of task points, with the unit of kilometers and pieces respectively; "[]" represents the matrix symbol; (2) Establish a spatial partitioning model for the dynamic route set characterization algorithm: First, use formulas (2) - (3) to obtain the abscissa and ordinate arrays mcz of the clustering center points of the task points: [~, cls] = kmeans_func([xrd; yrd]', ncl) (2) mcz = mean(cls) (3) In formulas (2) - (3), kmeans_func, ncl, cls, "~", and "'" are the function for solving the clustering center, its classification number, the planar coordinates of the obtained clustering center points, the symbol for not outputting the category number, and the transpose symbol respectively. mcz and mean are the abscissa and ordinate arrays of the clustering center points and the function for calculating the average value respectively; Secondly, use formulas (4) - (7) to obtain the spatial partitioning model cm corresponding to the dynamic route set characterization algorithm: r2 = [rand(c4, 1) * c1, rand(c4, 1) * c2, rand(c4, 1) * c3] (4) p0 = [x0, y0] pe = [xe, ye] In formulas (4) - (5), r2 and rand are the radius set corresponding to the dynamic route chord set and the function for generating random numbers respectively. c1, c2, c3, and c4 are the radius dynamic adjustment coefficients corresponding to dynamic route chords 1, 2, 3, and 4 respectively; x, y, cm(x, y), cu(x, y), and cd(x, y) are the planar coordinates, the value of the spatial partitioning model of the dynamic route set characterization algorithm at the coordinate (x, y), and the sub-model values of its upper left triangle and lower right triangle regions respectively; in formulas (6) - (7), p0, x0, y0, pe, xe, and ye are the lower left corner point of the planned area and its abscissa and ordinate, the upper right corner point and its abscissa and ordinate respectively; uxi and uyi are the abscissa and ordinate of the curve part of the dividing line in the upper left triangle region, with the unit of km; fliplr and end are the functions for reversing the elements of the array from left to right and the symbol for taking the last element in the array respectively, nu is the number of dividing regions in the upper left triangle region; dxi and dyi are the abscissa and ordinate of the curve part of the dividing line in the lower right triangle region, with the unit of km; nd is the number of dividing regions in the lower right triangle region; Finally, use formulas (8) - (9) to divide the task sets in each characterization region; [in, on] = inpolygon(rw, cm(x, y)) (8) frw = [frwu, frwd] = rw(in|on) (9) In formulas (8) to (9), rw, inpolygon, in, and on are respectively the set of task point coordinates, the function of the set of points included in the selected polygon, and the sets of points depicting the interior and boundary of the self-region selected by it. "|" represents the "or" operator in relational operations. frw, frwu, and frwd respectively represent the task set after segmentation, its upper left triangular subset, and its lower right triangular subset; (3) Complete the task allocation of the integrated route optimization method: Use formulas (10) to (12) to obtain the adaptive generation route set frwp of the task point set integrating the nearest neighbor sorting algorithm: frwup(i) = fitsortfun([p0, frwu(i), pe]) (10) frwdp(i) = fitsortfun([p0, frwu(end - i + 1), pe]) (11) frwp = [frwup(i), frwdp(i)], i = 1, 2,... nu + nd (12) In formulas (10) to (12), frwp, frwup, and frwdp are respectively the adaptive generation route set of the task point set integrating the nearest neighbor sorting algorithm, its upper triangular and lower triangular subsets, and fitsortfun is the function for finding the integrated nearest neighbor sorting; (4) Use the decision model to screen the optimal task allocation set: Use formulas (13) to (16) to obtain the optimal task allocation and route set zhp: nb = numel(frwp) (13) nbz = std(nb) (14) [mn, mnc, mnl] = sort(nbz) (15) zhp = frwp(mnc(1), mnl(1)) (16) In formulas (13) to (16), numel and nb are respectively the function for counting the number of array elements and the number of task points depicting the partition counted by it. std and nbz are respectively the function for calculating the standard deviation of the array and the standard deviation depicting the number of task points in the partition; sort, mn, mnc, and mnl are respectively the array element sorting function and the task point route sequence, the row number and column number corresponding to the original sequence after sorting, and zhp is the optimal task allocation and route set; (5) Select the optimal number of UAVs and calculate the flight range: First, use formula (17) to obtain the selected number of necessary UAVs: flyn = flyn + 1, length(zhp) > 0 (17) In formula (17), flyn and length are respectively the number of necessary UAVs and the function for calculating the maximum dimension number of the array; Secondly, use formulas (18) to (20) to obtain the flight range of the UAVs corresponding to each frame including sub-task points: tp = dist(zhp(j), zhp(j + 1)), j = 1, 2,... np (18) Dt = [Dt tp] (19) Dl(i) = sum(Dt), i = 1, 2,... nu + nd (20) In formulas (18) to (20), np, dist, tp, and Dt are the number of waypoints in a flight path, the function for calculating the distance between two points, the distance between waypoints j and j + 1, and the flight range of each UAV, respectively. sum and Dl are the function for summing the elements in an array and the sum of the flight ranges of all UAVs, respectively.

Images