Constraint optimization method for solving heterogeneous task allocation of multiple unmanned aerial vehicles

Abstract

The invention belongs to the technical field of intelligent optimization of intelligent robots, and particularly relates to a constraint optimization method for solving heterogeneous task allocation of multiple unmanned aerial vehicles. According to the method, a composite real number coding mode is adopted, and more information can be recorded in short individual codes; a self-adaptive alpha constraint level comparison strategy is adopted in the algorithm, relative feasible solution information can be effectively utilized in the iteration process, the convergence speed of the alpha level is adaptively influenced through the proportion of a relative feasible solution, the safety level is gradually improved, the search range of an optimal allocation scheme is limited within a satisfactory scheme, and the optimal allocation efficiency is improved. Therefore, an effective distribution scheme can be planned.

Description

Technical Field

[0001] The present invention belongs to the technical field of intelligent robot intelligent optimization, and specifically relates to a constrained optimization method for solving heterogeneous task allocation of multiple unmanned aerial vehicles. Background Art

[0002] In the increasingly complex and severe battlefield environment, the reconnaissance and monitoring of hostile targets and battlefield environment is an important basis for further decision-making. In recent years, the superior performance of drones in the battlefield has increasingly reflected their importance in the modern battlefield. However, the ability of a single drone to carry resources and perform tasks is limited, and it is difficult to meet the complex and changing task requirements. Therefore, the intelligent decision-making and planning of drone clusters composed of multiple drones and collaborative execution of tasks has become a research hotspot in the current drone field. And the task allocation of multiple drones has become the key to improving the autonomy and combat efficiency of the cluster. UAV task allocation refers to the allocation of one or a group of ordered tasks to individuals in a drone cluster based on certain environmental knowledge and task requirements, so as to achieve the best overall efficiency of the cluster while completing the maximum possible number of tasks. With the rapid development of drone technology and air defense countermeasure technology, relying solely on manual task allocation and flight planning can no longer meet the actual needs of modern warfare. Therefore, a large number of intelligent optimization algorithms that model this problem as an optimization problem have emerged, which have become an effective method to solve the problem of drone task allocation.

[0003] The differential evolution algorithm was first proposed by Storn and Price in 1995. It is a heuristic algorithm based on swarm intelligence with fast convergence speed and good robustness. It has maintained excellent performance in various constrained / unconstrained optimization problems. The basic steps of the algorithm are as follows: randomly initialize the population, generate a new population through mutation, crossover and selection operations, and is a random parallel global search algorithm with the characteristics of memorizing individual optimal solutions and sharing information within the population, that is, solving the optimization problem through cooperation and competition among individuals in the population. Its basic evolutionary idea is: generate mutation vectors according to the prescribed mutation rules, retain excellent solutions through greedy selection to avoid population degradation, generate new populations through selection operations, and iterate continuously until the optimal solution is approached.

[0004] Although this method has its advantages, it also has inevitable defects, which are mainly reflected in: the selection operation adopts a greedy strategy, discarding some feasible non-optimal solutions; the algorithm is prone to fall into local optimality and lacks a mechanism to jump out of the local optimality. Summary of the invention

[0005] The purpose of the present invention is to provide a constrained optimization method for solving heterogeneous task allocation of multiple unmanned aerial vehicles in view of the defects of the prior art.

[0006] The specific technical solution adopted by the present invention is as follows: A constrained optimization method for solving heterogeneous task allocation of multiple UAVs comprises the following steps:

[0007] Preparation: Assume that the flight speed of each drone in the drone cluster is v = 1m / s, and the parameter values ​​are as follows: Maximum number of iterations NC max =500; time window T for point target win =150; total time threshold T sum =450; Number of differential evolution populations N pop =10; scaling factor F(t) is a random number between [0.2,0.9], crossover factor C r , cross factor upper and lower limits C rmax =0.9, C rmin = 0.6, the coordinates of the starting point and the flight direction of the drone are as follows Figure 7 As shown, the number of drones M = 3, the number of tasks N = 8, and the starting point of the drone cluster {S 1 ,S 2 ,…S M}、Target point{T 1 ,T 2 ,…T N}, the type of each target and the turning radius of the drone are known, and the maximum number of iterations of the algorithm is NC max ,

[0008] Step 1: Calculate the dimension D = 2*N of the population individuals according to the target number N;

[0009] The dimension D = 2*N of the population individuals is calculated based on the number of targets N. Each target occupies two dimensions. The scope and meaning of each dimension are determined by the type of heterogeneous targets: an individual in the population is X(t) = {(x 1 ,x 2 ),(x 3 ,x 4 ),…(x D-1 ,x D )} Starting from the first dimension, every two dimensions represent the distribution of a target, and have different meanings depending on the target type. For line targets: (x 2i-1 ,x 2i ) After scaling, the integer part of the first dimension data represents the corresponding drone ID, and the size of the decimal part describes the sequence of the task performed by the drone. The second dimension data is scaled to determine the direction of its entry into the line target; for the surface target: (x 2i-1 ,x 2i) After scaling, the first dimension data is the same as the line target, and the second dimension data is scaled to determine the direction in which it enters the surface target. The difference between the line target and the surface target is that the entrance direction of the surface target is four, while the line target has two.

[0010] Step 2: Initialize the number of iterations t = 1, initialize the differential evolution population, and the number of populations is N pop , the individual dimension of the differential evolution population is D, and N is randomly generated pop The multi-UAV task allocation scheme is represented by the differential evolution population as a vector group X(t) = {X 1 (t),…,X i (t),...,X Npop (t)}, where X i (t) = {(x 1 ,x 2 ),(x 3 ,x 4 ),…(x D-1 ,x D )} represents the i-th individual, which corresponds to the i-th complete task allocation plan, i = 1, 2…, N pop , x j (t) The meaning is determined by the type and specific dimension of the corresponding target, j = 1, 2, ..., D;

[0011] Step 3: Obtain the mission planning sequence of the individual drone cluster through decoding operation as P = {[P 1,1 ,P 1,2 ,…,P 1,n1 ],[P 2,1 ,P 2,2 ,…,P 3,n2 ],…[P M,1 ,P M,2 ,…,P M,nM ]}, where [P 1,1 ,P 1,2 ,…,P 1,n1 ] represents the task allocation result of UAV 1; if P 1,1 =4, indicating that the first mission executed by UAV 1 is numbered 4; the subscript n1 indicates the total number of targets assigned to UAV 1 in this mission;

[0012] Step 4: Calculate the task execution time f of the task allocation scheme for the Durbin four-curve corresponding to the differentially evolved individuals T (X), satisfaction level μ(X);

[0013] Step 5: According to the comparison of the satisfaction level and safety level α of the individuals in the population, the proportion of the relatively feasible solutions of the population is obtained, and the adaptive safety level α is updated;

[0014] Step 6: Calculate the population individual X by the mutation operation of the differential evolution algorithm i The variation vector V of (t) i (t);

[0015] Step 7: According to the mapping rule, the mutation vector dimension that exceeds the boundary is mapped to the value range through random transformation;

[0016] Step 8: Mix the target vector and the mutation vector according to the crossover operation of the differential evolution algorithm to generate each individual X i The corresponding experimental vector U i (t);

[0017] Step 9: Calculate the fitness value f of the experimental vector U(t) T (U(t)) and satisfaction level μ(U(t));

[0018] Step 10: Perform the selection operation of the differential evolution algorithm through the α-level comparison method to form a new differential evolution population, and select the optimal individual X through the α-level comparison method best (t);

[0019] Step 11: Obtain the relative feasible solution ratio based on the satisfaction level of individuals in the population and the current adaptive safety level α, and update the cross factor C based on this value r (t) and control parameter β t ;

[0020] Step 12: Let t = t + 1, and determine whether t> NC max If it is not satisfied, return to step 4 for iteration; if it is satisfied, the operation is considered to be completed and the global optimal solution X is output. best (t);

[0021] Step 13: The global optimal solution X best (t) Decoding; outputting the obtained task allocation plan.

[0022] A constrained optimization method for solving heterogeneous task allocation of multiple UAVs as described above, wherein:

[0023] The encoding method described in step 1 is introduced in detail as follows: In order to solve the problem of allocating three types of heterogeneous reconnaissance targets, the present invention adopts a composite encoding method. Each target will be represented by two-dimensional data, and different information will be represented according to different target types. For point targets, the two-dimensional data respectively represent the labels of two drones targeting the point target; for line targets, the first dimension data represents the label of the drone, and the second dimension determines the direction of the drone entering the line target; similarly, the two-dimensional data representing the surface target respectively represent the label of the drone and the direction of entering the surface target. The entry direction for the line target includes two, and the entry direction for the surface target includes four different entry directions. This selection will further affect the exit position and exit direction for the target. The specific allocation is shown in the following table:

[0024]

[0025] The decoding calculation formula for the task execution order described in step 3 is as follows:

[0026]

[0027] Each individual obtains the execution task list of each drone according to the decoding formula (1), and then sorts the data representing the task in the list from large to small according to the decimals, and finally obtains the order in which the drone executes the task.

[0028] A constrained optimization method for solving heterogeneous task allocation of multiple UAVs as described above, wherein:

[0029] The UAV trajectory model based on the Dubins curve described in step 4 is introduced in detail as follows:

[0030] In the heterogeneous task allocation problem considered, multiple UAVs need to conduct reconnaissance on heterogeneous targets under kinematic constraints. To simplify the calculation, we make the following assumptions:

[0031] (1) The speed of each drone is constant;

[0032] (2) The UAV performs a reconnaissance mission at a given altitude;

[0033] (3) Drones fly at different altitudes to avoid collisions between drones;

[0034] (4) The fuel carried by each UAV is limited by the final mission time. In addition, the Dubins aircraft model is used to describe the motion constraints of each UAV as follows:

[0035]

[0036] Where (x, y) is the horizontal position in the Cartesian inertial reference frame; θ is the heading angle; v is the flight speed; rk is the minimum turning radius; c is the control input, c>0 means the UAV turns left, c<0 means the UAV turns right, c=0 means the UAV goes straight, and c=±1 means the UAV turns with the minimum turning radius. Therefore, the vector (x, y, θ) can be used to represent the state of the UAV in the battlefield environment.

[0037] The Dubins path is used to generate the connection path between targets. Considering the requirements of each target for the terminal heading in the modeling scenario of this paper, the Dubins curve in this paper will take into account the heading constraint of the terminal.

[0038] The Dubin curve with terminal heading constraint refers to the shortest path for a machine with the minimum turning radius to move from the initial state to the next state. Its shortest path contains two segments with radius r k An arc and a straight line, and it must be one of the following path sets {RSL,RSR,LSR,LSL,RLR,LRL}, where R represents a clockwise arc, L represents a counterclockwise arc, and S represents a straight line segment.

[0039] As described above, a constrained optimization method for solving heterogeneous task allocation of multiple UAVs is described, wherein the task execution time f of the task allocation scheme corresponding to the differentially evolved individuals is calculated. T (X), satisfaction level μ(X), which is calculated as follows:

[0040]

[0041] The above formula L ij It represents the path length that drone i flies when executing the jth target, which includes two parts: the path length from the previous target or the starting point to the target and the flight time required to execute the task. The flight time required to execute three heterogeneous targets is different. In order to improve the calculation reliability of the flight path length of drones in the task allocation scenario, the present invention abandons the Euclidean distance used in task allocation under normal circumstances, and calculates the path length value strictly according to the Dubins curve;

[0042]

[0043] in represents the flight distance of UAV i from the last target point or starting point to perform task j, It represents the flight distance required to perform mission j. Different target types have different flight distances. To perform a point target mission, the drone needs to fly around the target point. The flight distance of a line target is determined by the length of the line target. The flight distance of a surface target is determined by the length and width of the surface target. The basic reconnaissance model is as follows: Figure 2 , 3 , 4,

[0044] The equality constraints and inequality constraints described in step 4 are described in detail as follows:

[0045] Time window constraint: The exploration of point targets requires two drones to conduct reconnaissance. The arrival time of the two drones is t i , t j A certain time window must be satisfied. The inequality constraint is described as g 1 ,

[0046] g 1 =|t i -t j |-T win (8)

[0047] Total time threshold constraint: There is also a limit on the completion time of the entire task. This constraint can be understood as the range constraint of the drone due to the limited fuel it carries. This inequality constraint is described as g 2 , the inequality constraint calculation is given by formula (9),

[0048] g 2 =|F|-T sum (9)

[0049] Redundancy constraint: From the task allocation scenario setting, it can be seen that the allocation between drones and targets must meet certain redundancy requirements. Each drone must be assigned to one target. At the same time, the number of drones assigned to each target must meet the requirements of its task type. For example, a line target requires two drones while a surface target only requires one drone. It is necessary to determine whether there is redundancy or lack in the number of drones required for each task. The equality constraint is described as h:

[0050]

[0051] where c ij is the decision variable representing the allocation relationship between UAV i and task j, c ij If it is 1, it means that drone i performs task j, N limj represents the number of drones required for task j, where point target N lim =2, Line target and area target N lim =1; M represents the total number of drones performing reconnaissance missions, N represents the total number of reconnaissance missions, 5. A constrained optimization method for solving heterogeneous task allocation of multiple drones according to claim 1, characterized in that: safety level α t Adaptive convergence rule for :

[0052]

[0053] To make the safety level α tIt can slow down or accelerate convergence according to the evolution of the population, improve global and local search capabilities, and adjust the parameter β t Parameters and the relative feasible solution ratio k of the population flag Establish a connection, the relatively feasible solution refers to individual X i The satisfaction level μ(X i )≥α, the infeasible solution refers to the individual X i The satisfaction level μ(X i )<α, β t To control the safety level α t The key coefficient of the convergence rate, its value will directly determine α t Compared with the increment of the parent generation, the individual satisfaction level μ in the current population is greater than the safety level α t The proportion of the number of species in a population.

[0054] A constrained optimization method for solving heterogeneous task allocation of multiple UAVs as described above, wherein the adaptive crossover factor and selection operation in the differential evolution algorithm are improved:

[0055] The crossover operation of the differential evolution algorithm described in step 8 is described in detail as follows:

[0056]

[0057] Among them C r (t) represents the crossover factor, r n The random number index that changes in each round ensures that there is at least one crossover operation in each round, otherwise the population will not change each time. rmax With C rmin represents the upper and lower bounds of the crossover factor, and C r The size of (t) determines the probability of crossover, and the selection of its size has a huge impact on the evolution rate of the population. In this invention, the feasibility ratio is continued to be used to determine the crossover constant k. flag Make adaptive changes. When the proportion of relatively feasible solutions in the population is large, it means that there are fewer violations of constraints in the current population, and the population may converge to the optimal solution from a larger possibility. Then, when performing crossover operations, performing smaller mutations and performing detailed dimensional searches are more conducive to fine convergence to the optimal solution. On the contrary, if the proportion of relatively feasible solutions in the current population is small, it means that the current population has serious violations of constraints. Then the population's evolutionary step size and the crossover probability of the mutation vector should tend to be larger, which is more conducive to jumping out of infeasible solutions.

[0058] The selection operation of the differential evolution algorithm under the α level comparison described in step 10 is described in detail as follows:

[0059]

[0060] Assume f(X 1 )、f(X 2 ) and μ(X 1 ), μ(X 1 ) are X 1 and X 1 The fitness value and satisfaction of individual X 1 ,X 2 The alpha level comparison between is defined as follows:

[0061]

[0062] According to the α-level comparison, only when the satisfaction of two solutions is higher than the predefined value or the satisfaction is the same, the comparison between the two solutions will be made according to the objective function value. According to the α-level comparison method, the ordinary comparison is replaced by the α-level comparison, thereby converting the general constrained minimization problem into a minimization problem under the α-level comparison.

[0063] The present invention has the following beneficial effects compared with the prior art: The present invention proposes an improved adaptive α-level comparison differential evolution algorithm for solving the collaborative heterogeneous constrained task allocation of multiple UAVs. The proposed method adopts a composite coding method, which can contain the timing allocation scheme of complex heterogeneous tasks within a limited dimension. The present invention adopts the Dubins curve model to calculate the flight path length of the UAV between task points, which greatly improves the accuracy of task allocation evaluation compared with the Euclidean distance model. The algorithm uses an improved adaptive α-level comparison to handle the complex spatiotemporal constraints in the heterogeneous task allocation scenario, and improves the crossover operation and selection operation of the benchmark differential evolution algorithm to greatly improve the algorithm's optimization ability, thereby planning a more efficient task allocation scheme. Based on the above discussion, the proposed improved adaptive α-level comparison differential evolution algorithm (aDELC) has good optimization performance, high convergence accuracy, and strong robustness. It is an efficient allocation method in solving the collaborative heterogeneous constrained task allocation scenario of multiple UAVs. BRIEF DESCRIPTION OF THE DRAWINGS

[0064] Figure 1 It is a schematic diagram of the six types of Dubin curve trajectories.

[0065] Figure 2 It is a schematic diagram of the UAV’s execution point target path.

[0066] Figure 3 It is a schematic diagram of the target path of the UAV execution line.

[0067] Figure 4 It is a schematic diagram of the UAV’s execution surface target path.

[0068] Figure 5It is the flow chart of the improved adaptive α-level comparison differential evolution algorithm.

[0069] Figure 6 This is the result graph of heterogeneous task allocation based on the aDELC algorithm.

[0070] Figure 7 Heterogeneous task distribution scenario data table in battlefield environment

[0071] Figure 8 It is the optimal result of heterogeneous task allocation based on the aDELC algorithm.

[0072] Fig. 9 It is a statistical data table of heterogeneous task allocation based on the aDELC algorithm.

[0073] Fig.10 It is the evolution curve of average task execution time. DETAILED DESCRIPTION

[0074] The present invention is further described below in conjunction with a specific implementation process. A constrained optimization method for solving heterogeneous task allocation for multiple UAVs is provided.

[0075] Preparation: Assume that the flight speed of each drone in the drone cluster is v = 1m / s. The parameter values ​​are as follows: Maximum number of iterations NC max =500; time window T for point target win =150; total time threshold T sum =450; Number of differential evolution populations N pop =10; scaling factor F(t) is a random number between [0.2,0.9], crossover factor C r , cross factor upper and lower limits C rmax =0.9, C rmin = 0.6. The coordinates of the starting point and the flight direction of the drone are as follows Figure 7 As shown, the number of drones M = 3, and the number of tasks N = 8. The starting point of the drone cluster {S 1 ,S 2 ,…S M}、Target point{T 1 ,T 2 ,…T N}, the type of each target and the turning radius and reconnaissance radius of the drone are all known. The maximum number of iterations of the algorithm NC max .

[0076] Step 1: Calculate the dimension D = 2*N of the population individuals according to the target number N;

[0077] The dimension D = 2*N of the population individuals is calculated based on the number of targets N. Each target occupies two dimensions. The scope and meaning of each dimension are determined by the type of heterogeneous targets: an individual in the population is X(t) = {(x 1 ,x 2 ),(x 3 ,x 4 ),…(x D-1 ,x D )} Starting from the first dimension, each two-dimensional data represents the distribution of a target, and has different meanings depending on the target type. For line targets: (x 2i-1 ,x 2i ) After scaling, the integer part of the first dimension data represents the corresponding drone ID, and the size of the decimal part describes the sequence of the task performed by the drone. The second dimension data is scaled to determine the direction of its entry into the line target; for the surface target: (x 2i-1 ,x 2i )After scaling, the first-dimensional data is the same as the line target, and the second-dimensional data is scaled to determine its direction of entry into the surface target. What is different from the line target is that the surface target has four entry directions while the line target has two.

[0078] The encoding method described in step one is introduced in detail as follows: In order to solve the problem of allocating three types of heterogeneous reconnaissance targets, the present invention adopts a composite encoding method. Each target will be represented by two-dimensional data, and different information will be represented according to different target types. For point targets, the two-dimensional data respectively represent the numbers of the two drones that target the point target; for line targets, the first dimension of data represents the number of the drone, and the second dimension determines the direction in which the drone enters the line target; similarly, the two-dimensional data representing the surface target respectively represent the number of the drone and the direction in which it enters the surface target. There are two entry directions for line targets and four different entry directions for surface targets. This choice will further affect the exit position and exit direction of the target, and its schematic diagram is attached. Figure 3 and attached Figure 4 , E1 to E4 in the figure represent their entry numbers. The specific allocation is shown in the following table:

[0079]

[0080]

[0081] Step 2: Initialize the number of iterations t = 1, initialize the differential evolution population, and the number of populations is N pop , the individual dimension of the differential evolution population is D, and N is randomly generated pop The multi-UAV task allocation scheme is represented by the differential evolution population as a vector group X(t) = {X 1 (t),…,X i(t),...,X Npop (t)}, where X i (t) = {(x 1 ,x 2 ),(x 3 ,x 4 ),…(x D-1 ,x D )} represents the i-th individual, which corresponds to the i-th complete task allocation plan, i = 1, 2…, N pop , x j (t) The meaning is determined by the type and specific dimension of the corresponding target, j = 1, 2, ..., D;

[0082] Step 3: Obtain the mission planning sequence of the individual drone cluster through decoding operation as P = {[P 1,1 ,P 1,2 ,…,P 1,n1 ],[P 2,1 ,P 2,2 ,…,P 3,n2 ],…[P M,1 ,P M,2 ,…,P M,nM ]}, where [P 1,1 ,P 1,2 ,…,P 1,n1 ] represents the task allocation result of UAV 1; if P 1,1 =4, indicating that the first mission executed by UAV 1 is numbered 4; subscript n1 indicates the total number of targets assigned to UAV 1 in this mission, and the decoding operation refers to formula (5);

[0083] The decoding calculation formula for the task execution order described in step 3 is as follows:

[0084]

[0085] Each individual obtains the execution task list of each drone according to the decoding formula (5), and then sorts the data representing the task in the list from large to small according to the decimals, and finally obtains the order in which the drone executes the task.

[0086] An example is used to illustrate the encoding and decoding process used in this paper: the number of drones M = 2, the number of trees N = 4, and the task types are point target, point target, line target, and surface target. Attached figure (5) describes the decoding process of the variables in detail. After scaling and decoding, we get the corresponding information represented by each dimension. The execution tasks of each drone are sorted from small to large according to the size of the decimal. This sorting is the execution order of the drones performing reconnaissance tasks. At the same time, the entry directions of line targets and surface targets are also analyzed during the decoding process, which are represented by red fonts in the figure. Therefore, the final task execution order is decoded as drone 1 performing tasks 1-3 (2) -2 respectively, and the entry direction of line target No. 3 is 2. Similarly, the task execution order of drone 2 is 2-1-4 (1), and the entry direction of surface target No. 4 is 1. It can be represented by the vector P = {[1,3,2], [2,1,4]}.

[0087] Step 4: Calculate the execution time f of the task allocation scheme corresponding to the differentially evolved individuals T (X), satisfaction level μ(X) formula is as follows:

[0088]

[0089] The above formula L ij It represents the path length that drone i flies when executing the jth target, which includes two parts: the path length from the previous target or the starting point to the target and the flight time required to execute the task. The flight time required to execute three heterogeneous targets is different. In order to improve the calculation reliability of the flight path length of drones in the task allocation scenario, the present invention abandons the Euclidean distance used in task allocation under normal circumstances, and calculates the path length value strictly according to the Dubins curve;

[0090] The satisfaction of a solution is obtained by combining the satisfaction of all constraints according to the following formula:

[0091] μ(X i )=min{μ g1 (X i ),μ g2 (X i ),μ h (X i )} (2)

[0092]

[0093] In the above formula Indicates that the inequality constraints in the task allocation scenario are calculated with reference to formula (8) and formula (9), μ h The equality constraint in the task allocation scenario is calculated according to formula (10). The constraint satisfaction level depends on the parameter b 1 ,b2 ,b 3 The correct settings of the three parameters depend on the specific problem. In the present invention, the median of the corresponding constraint satisfaction in the initial population is used, because the median of the randomly initialized constraint satisfaction can show the appropriate degree of constraint violation.

[0094] The UAV trajectory model based on the Dubins curve described in step 4 is introduced in detail as follows:

[0095] In the heterogeneous task allocation problem considered, multiple UAVs need to conduct reconnaissance on heterogeneous targets under kinematic constraints. In this paper, to simplify the calculation, we make the following assumptions:

[0096] (1) The speed of each drone is constant;

[0097] (2) The UAV performs a reconnaissance mission at a given altitude;

[0098] (3) Drones fly at different altitudes to avoid collisions between drones;

[0099] (4) The fuel carried by each UAV is limited by the final mission time. In addition, the Dubins aircraft model is used to describe the motion constraints of each UAV as follows:

[0100]

[0101] Where (x, y) is the horizontal position in the Cartesian inertial reference frame; θ is the heading angle; v is the flight speed; r k is the minimum turning radius; c is the control input. c>0 means the drone turns left, c<0 means the drone turns right, and c=0 means the drone goes straight. c=±1 means the drone turns with the minimum turning radius. Therefore, the vector (x, y, θ) can be used to represent the state of the drone in the battlefield environment.

[0102] The Dubins path is used to generate the connection path between targets. Considering the requirements of each target for the terminal heading in the modeling scenario of this article, the Dubins curve in this article will take into account the heading constraints of the terminal.

[0103] The Dubin curve with terminal heading constraint refers to the shortest path for a machine with the minimum turning radius to move from the initial state to the next state. Its shortest path contains two segments with radius r k An arc and a straight line, and it must be one of the following path sets {RSL,RSR,LSR,LSL,RLR,LRL}, where R represents a clockwise arc, L represents a counterclockwise arc, and S represents a straight line segment. Figure 1 The six path situations mentioned above are shown in FIG.

[0104] Step 4 Lij It represents the length of the path that UAV i flies when executing the jth target. The specific calculation is as follows:

[0105]

[0106] in represents the flight distance of UAV i from the last target point or starting point to perform task j, It represents the flight distance required to execute mission j. The flight distance varies depending on the target type. To execute a point target mission, the drone needs to fly around the target point. The flight distance of a line target needs to be determined by the length of the line target, and the flight distance of a surface target will be determined by the length and width of the surface target. Its basic reconnaissance model is as follows: Figure 2 , 3 , as shown in Figure 4.

[0107] The equality constraints and inequality constraints described in step 4 are described in detail as follows:

[0108] Time window constraint: The exploration of point targets requires two drones to conduct reconnaissance. The arrival time of the two drones is t i , t j A certain time window must be satisfied. The inequality constraint is described as g 1 .

[0109] g 1 =|t i -t j |-T win (8)

[0110] Total time threshold constraint: There is also a limit on the completion time of the entire task. This constraint can be understood as the range constraint of the drone due to the limited fuel it carries. This inequality constraint is described as g 2 , the inequality constraint calculation is given by formula (9).

[0111] g 2 =|F|-T sum (9)

[0112] Redundancy constraint: From the task allocation scenario setting, it can be seen that the allocation between drones and targets must meet certain redundancy requirements. Each drone must be assigned to a target, and the number of drones assigned to each target must meet the requirements of its task type. For example, a line target requires two drones while a surface target only requires one drone. To determine whether there is redundancy or lack of the number of drones required for each task, the equality constraint is described as h:

[0113]

[0114] where c ijis the decision variable representing the allocation relationship between UAV i and task j, c ij If it is 1, it means that UAV i performs the j task, N limj represents the number of drones required for task j, where point target N lim =2, Line target and area target N lim =1; M represents the total number of UAVs performing reconnaissance missions, and N represents the total number of reconnaissance missions.

[0115] Step 5: Update the adaptive security level α, the formula is as follows:

[0116]

[0117] In the above formula, t is the current iteration number, NC max is the maximum number of iterations, β(t) is the convergence coefficient of the control safety level α(t), and its value is adaptively changed according to the proportion of relatively feasible solutions in the present invention, and is given by formula (12); μ 0 The value is the median of the constraint satisfaction of all individuals in the initial population.

[0118] The adaptive convergence of the security level α(t) of the t-th iteration described in step 5 is described in detail as follows:

[0119]

[0120] In the present invention, in order to make the security level α slow down or accelerate the convergence according to the evolution of the population and improve the global search ability and local search ability, the parameter β is t Parameters and the relative feasible solution ratio k of the population flag Establish a connection, the relatively feasible solution refers to individual X i The satisfaction level μ(X i )≥α, the infeasible solution refers to the individual X i The satisfaction level μ(X i )<α. t k is the key coefficient for controlling the convergence rate of security level α, and its value will directly determine the increment of α(t) compared with α(t-1). flag is the individual satisfaction level μ(X i ) is greater than the safety level α(t) in the population as shown in formula (13).

[0121] Step 6: Calculate the population individual X by the mutation operation of the differential evolution algorithm i The variation vector V of (t) i (t), refer to formula (14);

[0122] The mutation operation of the differential evolution algorithm described in step 6 is shown in the following formula:

[0123]

[0124] Among them, r 1 、r 2 、r 3 is a random index, but it must be ensured that they are different from each other. Fc∈[0,2] is the scaling factor that controls the mutation vector. If a dimension of the mutation vector exceeds the value range, it is randomly transformed to the range, as shown in formula (15); where V i,j Represents the j-th dimension of the i-th mutation vector.

[0125] V i,j =rand(0,1),if V i,j <0orV i,j >1 (15)

[0126] Step 7: According to the mapping rule, the mutation vector dimension that exceeds the boundary is mapped to the value range through random transformation, refer to formula (15);

[0127] Step 8: Mix the target vector and the mutation vector according to the crossover operation of the differential evolution algorithm to generate each individual X i The corresponding experimental vector U i (t), refer to formula (16);

[0128] The crossover operation of the differential evolution algorithm described in step 8 is described in detail as follows:

[0129]

[0130] Among them C r (t) represents the crossover factor, r n The random number index that changes in each round ensures that there is at least one crossover operation in each round, otherwise the population will not change each time. rmax With C rmin represents the upper and lower bounds of the crossover factor, and C r The size of (t) determines the probability of crossover, and the selection of its size has a huge impact on the evolution rate of the population. In this invention, the feasibility ratio is continued to be used to determine the crossover constant k. flag Make adaptive changes. When the proportion of relatively feasible solutions in the population is large, it means that there are fewer violations of constraints in the current population, and the probability of convergence to the optimal solution is greater. In this case, performing smaller mutations and performing detailed dimensional searches during the crossover operation is more conducive to fine convergence to the optimal solution. On the contrary, if the proportion of relatively feasible solutions in the current population is small, it means that the current population has serious violations of constraints, and the evolutionary step length of the population and the crossover probability of the mutation vector should tend to be larger, which is more conducive to jumping out of infeasible solutions.

[0131] Step 9: Calculate the fitness value f of the experimental vector U(t) by formula (1) and formula (2) T (U(t)) and satisfaction level μ(U(t));

[0132] Step 10: Perform the selection operation of the differential evolution algorithm through the α-level comparison method to form a new differential evolution population, and select the optimal individual X through the α-level comparison method best (t), refer to formula (17);

[0133] The selection operation of the differential evolution algorithm under the α level comparison described in step 10 is described in detail as follows:

[0134]

[0135] Assume that f(X 1 )、f(X 2 ) and μ(X 1 ), μ(X 1 ) are X 1 and X 1 The fitness value and satisfaction of individual X 1 ,X 2 The alpha level comparison between is defined as follows:

[0136]

[0137]

[0138] According to the α-level comparison, only when the satisfaction of two solutions is higher than the predefined value or the satisfaction is the same, the comparison between the two solutions will be made according to the objective function value. According to the α-level comparison method, the ordinary comparison is replaced by the α-level comparison, thereby converting the general constrained minimization problem into a minimization problem under the α-level comparison.

[0139] Step 11: Obtain the relative feasible solution ratio k based on the satisfaction level of individuals in the population and the current adaptive safety level α(t) flag , and update the crossover factor C according to this value r (t) and control parameter β t ; Refer to formula (11), formula (12) and formula (13).

[0140] Step 11: Crossover Factor C r (t) and control parameter β t The update of will be calculated according to formula (16) and formula (13).

[0141] Step 12: Update the iteration period t=t+1 and determine whether t>NC is satisfiedmax If it is not satisfied, return to step 4 for iteration; if it is satisfied, the operation is considered to be completed and the global optimal solution X is output. best (t);

[0142] The selection of the best individual in step 12 is also performed in accordance with the α-level comparison method.

[0143] Step 13: The global optimal solution X best (t) Decoding; outputting the obtained task allocation plan.

[0144] A specific example is given below.

[0145] by Figure 3 The following environment is used as an example:

[0146] Preparation: Assume that the flight speed of each drone in the drone cluster is v = 1m / s. The parameter values ​​are as follows: Maximum number of iterations NC max =500; time window T for point target win =150; total time threshold T sum =450; Number of differential evolution populations N pop =10; the scaling factor F(t) is a random number between [0.2,0.9], and the upper and lower limits of the crossover factor C rmax =0.9, C rmin = 0.6. The coordinates of the starting point and the flight direction of the drone are as follows Figure 7 As shown, the number of drones M = 3, and the number of missions N = 8.

[0147] Step 1: Calculate the dimension of the population individuals D = 2*N = 16 according to the target number N;

[0148] Step 2: Set the number of iterations t = 1, initialize the differential evolution population, the individual dimension of the differential evolution population is D = 2*N, and randomly generate N pop = 10 groups of multi-UAV task allocation schemes, the differential evolution population is represented by a vector group X(t) = {X 1 (t),X 2 (t),...,X Npop (t)}, where X i (t) = {(x 1 ,x 2 ),(x 3 ,x 4 ),…(x D-1 ,x D )} represents the i-th individual, that is, corresponding to the i-th complete task allocation plan, i = 1, 2…, 10, X j(t) The meaning j=1, 2, ..., 16 is determined by the type and specific dimension of the corresponding target and is calculated by the decoding formula (5);

[0149] Step 3: Through the decoding operation, the UAV swarm mission planning sequence of the individual representative is obtained as P = {[P 1,1 ,P 1,2 ,…,P 1,n1 ],[P 2,1 ,P 2,2 ,…,P 3,n2 ],…[P M,1 ,P M,2 ,…,P M,nM ]}, where [P 1,1 ,P 1,2 ,…,P 1,n1 ] represents the task planning and allocation result of UAV 1; if P 1,1 =4, indicating that the first mission executed by UAV 1 is numbered 4; the subscript n1 indicates the total number of targets assigned to UAV 1 in this mission;

[0150] Step 4: Calculate the task execution time f of the task allocation scheme corresponding to the differentially evolved individuals T (X), satisfaction level μ(X) formula is as follows:

[0151]

[0152] The above formula L ij It represents the path length that drone i flies when executing the jth target, which includes two parts: the path length from the previous target or the starting point to the target and the flight time required to execute the task. The flight time required to execute three heterogeneous targets is different. In order to improve the calculation reliability of the flight path length of drones in the task allocation scenario, the present invention abandons the Euclidean distance used in task allocation under normal circumstances, and calculates the path length value strictly according to the Dubins curve;

[0153] The satisfaction of a solution is obtained by combining the satisfaction of all constraints according to the following formula:

[0154]

[0155]

[0156] In the above formula, μ g1 μ g2 Indicates that the inequality constraints in the task allocation scenario are calculated with reference to formula (8) and formula (9), μ h The equality constraint in the task allocation scenario is calculated according to formula (10). The constraint satisfaction level depends on the parameter b 1 ,b2 ,b 3 The correct settings of the three parameters depend on the specific problem. In the present invention, the median of the corresponding constraint satisfaction in the initial population is used, because the median of the randomly initialized constraint satisfaction can show the appropriate degree of constraint violation.

[0157] Step 5: Update the adaptive security level α, the formula is as follows:

[0158]

[0159] In the above formula, t is the current iteration number, NC max =500 is the maximum number of iterations, β(t) is the convergence coefficient of the control safety level α(t), and its value is adaptively changed according to the proportion of relatively feasible solutions in the present invention, and is given by formula (12); μ 0 The value is the median of the constraint satisfaction of all individuals in the initial population.

[0160] Step 6: Calculate the population individual X by the mutation operation of the differential evolution algorithm i The variation vector V of (t) i (t), refer to formula (14);

[0161] Step 7: According to the mapping rule, the mutation vector dimension that exceeds the boundary is mapped to the value range through random transformation, refer to formula (15);

[0162] Step 8: Mix the target vector and the mutation vector according to the crossover operation of the differential evolution algorithm to generate the experimental vector U i (t), refer to formula (16);

[0163] Step 9: Calculate the fitness value f of the experimental vector U(t) by formula (1) and formula (2) T (U(t)) and satisfaction level μ(U(t));

[0164] Step 10: Perform the selection operation of the differential evolution algorithm through the α-level comparison method to form a new differential evolution population, and select the optimal individual X through the α-level comparison method best (t), refer to formula (17);

[0165] Step 11: Obtain the relative feasible solution ratio k based on the satisfaction level of individuals in the population and the current adaptive safety level α(t) flag , and update the crossover factor C according to this value r (t) and control parameter β t ; Refer to formula (11), formula (12) and formula (13).

[0166] Step 12: Update the iteration period t=t+1 and determine whether t>NC is satisfied max If it is not satisfied, return to step 4 for iteration; if it is satisfied, the operation is considered to be completed and the global optimal solution X is output. best (t);

[0167] Step 13: The global optimal solution X best (t) decoding; outputting the obtained task allocation plan;

[0168] Fig. 9 shows the numerical results of 50 independent assignments of the aDELC algorithm, Fig.10 The evolution curve of the average time consumption of the task allocation scheme planned by the aDELC algorithm for 50 independent operations. The heterogeneous task allocation scheme planned by the aDELC algorithm realizes the redundancy requirements of the tasks for UAVs; Fig. 9 The data statistics of the algorithm in 50 independent operations are shown, where the second column is the shortest task plan time allocated in 50 optimizations, the third column is the average time of the task plan allocated in 50 optimizations, the fourth column is the median time of the task plan allocated in 50 optimizations, the fifth column is the longest time of the task plan allocated in 50 optimizations, the sixth column is the standard deviation of the time of the task plan allocated in 50 optimizations, and the seventh column is the success rate of the task plan allocated in 50 optimizations that does not violate the constraints. Fig. 9 and Fig.10 It can be seen that the use of the aDELC algorithm for heterogeneous task allocation can plan a reasonable allocation plan that meets the task redundancy requirements within the allowed time window and total time threshold constraints, and the success rate of the algorithm is maintained at 100% with good robustness. Therefore, the adaptive differential evolution algorithm based on α-level comparison proposed in the present invention is an effective multi-UAV heterogeneous task allocation method, which can stably plan a better UAV task execution allocation plan and ensure that the UAV cluster meets the overall task time threshold and the time window constraints of special targets.

Claims

1. A constrained optimization method for solving heterogeneous task allocation of multiple UAVs, characterized in that: The steps include: Preparation: Assume that the flight speed of each drone in the drone cluster is v = 1m / s, and the parameter values ​​are as follows: Maximum number of iterations NC max =500; time window T for point target win =150; total time threshold T sum =450; Number of differential evolution populations N pop =10; The scaling factor F(t) is a random number between [0.2, 0.9], and the crossover factor C r , cross factor upper and lower limits C rmax =0.9, C rmin = 0.6, the coordinates of the starting point and the flight direction of the drone are shown in Figure 7, the number of drones M = 3, the number of tasks N = 8, and the starting point of the drone cluster {S1, S2, ... S M }、Target point {T1,T2,…T N }, the type of each target and the turning radius of the drone are known, and the maximum number of iterations of the algorithm is NC max , Step 1: Calculate the dimension D = 2*N of the population individuals according to the target number N; The dimension D = 2*N of the population individuals is calculated based on the number of targets N. Each target occupies two dimensions. The scope and meaning of each dimension are determined by the type of heterogeneous targets: an individual in the population is X(t) = {(x1, x2), (x3, x4), … (x D-1 ,x D )} Starting from the first dimension, every two dimensions represent the distribution of a target, and have different meanings depending on the target type. For line targets: (x 2i-1 ,x 2i ) After scaling, the integer part of the first dimension data represents the corresponding drone ID, and the size of the decimal part describes the sequence of the task performed by the drone. The second dimension data is scaled to determine the direction of its entry into the line target; for the surface target: (x 2i-1 ,x 2i ) After scaling, the first dimension data is the same as the line target, and the second dimension data is scaled to determine the direction in which it enters the surface target. The difference between the line target and the surface target is that the entrance direction of the surface target is four, while the line target has two. Step 2: Initialize the number of iterations t = 1, initialize the differential evolution population, and the number of populations is N pop , the individual dimension of the differential evolution population is D, and N is randomly generated pop The multi-UAV task allocation scheme is represented by the differential evolution population as a vector group X(t) = {X1(t),…,X i (t),...,X Npop (t)}, where X i (t)={(x1,x2),(x3,x4),…(x D-1 ,x D )} represents the i-th individual, which corresponds to the i-th complete task allocation plan, i = 1, 2…, N pop , x j (t) The meaning is determined by the type and specific dimension of the corresponding target, j = 1, 2, ..., D; Step 3: Obtain the mission planning sequence of the individual drone cluster through decoding operation as P = {[P 1,1 ,P 1,2 ,…,P 1,n1 ],[P 2,1 ,P 2,2 ,…,P 3,n2 ],…[P M,1 ,P M,2 ,…,P M,nM ]}, where [P 1,1 ,P 1,2 ,…,P 1,n1 ] represents the task allocation result of UAV 1; if P 1,1 =4, indicating that the first mission executed by UAV 1 is numbered 4; the subscript n1 indicates the total number of targets assigned to UAV 1 in this mission; Step 4: Calculate the task execution time f of the task allocation scheme for the Durbin four-curve corresponding to the differentially evolved individuals T (X), satisfaction level μ(X); Step 5: According to the comparison of the satisfaction level and safety level α of the individuals in the population, the proportion of the relatively feasible solutions of the population is obtained, and the adaptive safety level α is updated; Step 6: Calculate the population individual X by the mutation operation of the differential evolution algorithm i The variation vector V of (t) i (t); Step 7: According to the mapping rule, the mutation vector dimension that exceeds the boundary is mapped to the value range through random transformation; Step 8: Mix the target vector and the mutation vector according to the crossover operation of the differential evolution algorithm to generate each individual X i The corresponding experimental vector U i (t); Step 9: Calculate the fitness value f of the experimental vector U(t) T (U(t)) and satisfaction level μ(U(t)); Step 10: Perform the selection operation of the differential evolution algorithm through the α-level comparison method to form a new differential evolution population, and select the optimal individual X through the α-level comparison method best (t); Step 11: Obtain the relative feasible solution ratio based on the satisfaction level of individuals in the population and the current adaptive safety level α, and update the cross factor C based on this value r (t) and control parameter β t ; Step 12: Let t = t + 1, and determine whether t> NC max If it is not satisfied, return to step 4 for iteration; if it is satisfied, the operation is considered to be completed and the global optimal solution X is output. best (t); Step 13: The global optimal solution X best (t) Decoding; outputting the obtained task allocation plan.

2. A constrained optimization method for solving heterogeneous task allocation of multiple UAVs according to claim 1, characterized in that: The encoding method described in step 1 is introduced in detail as follows: In order to solve the problem of allocating three types of heterogeneous reconnaissance targets, the present invention adopts a composite encoding method. Each target will be represented by two-dimensional data, and different information will be represented according to different target types. For point targets, the two-dimensional data respectively represent the labels of two drones targeting the point target; for line targets, the first dimension data represents the label of the drone, and the second dimension determines the direction of the drone entering the line target; similarly, the two-dimensional data representing the surface target respectively represent the label of the drone and the direction of entering the surface target. The entry direction for the line target includes two, and the entry direction for the surface target includes four different entry directions. This selection will further affect the exit position and exit direction for the target. The specific allocation is shown in the following table: The decoding calculation formula for the task execution order described in step 3 is as follows: Each individual obtains the execution task list of each drone according to the decoding formula (1), and then sorts the data representing the task in the list from large to small according to the decimals, and finally obtains the order in which the drone executes the task.

3. The constrained optimization method for solving heterogeneous task allocation of multiple UAVs according to claim 2 is characterized by: The UAV trajectory model based on the Dubins curve described in step 4 is introduced in detail as follows: In the heterogeneous task allocation problem considered, multiple UAVs need to conduct reconnaissance on heterogeneous targets under kinematic constraints. To simplify the calculation, we make the following assumptions: (1) The speed of each drone is constant; (2) The UAV performs a reconnaissance mission at a given altitude; (3) Drones fly at different altitudes to avoid collisions between drones; (4) The fuel carried by each UAV is limited by the final mission time. In addition, the Dubins aircraft model is used to describe the motion constraints of each UAV as follows: Where (x, y) is the horizontal position in the Cartesian inertial reference frame; θ is the heading angle; v is the flight speed; r k is the minimum turning radius; c is the control input, c>0 means the drone turns left, c<0 means the drone turns right, c=0 means the drone goes straight, and c=±1 means the drone turns with the minimum turning radius. Therefore, the vector (x, y, θ) can be used to represent the state of the drone in the battlefield environment. The Dubins path is used to generate the connection path between targets. Considering the requirements of each target for the terminal heading in the modeling scenario of this paper, the Dubins curve in this paper will take into account the heading constraint of the terminal. The Dubin curve with terminal heading constraint refers to the shortest path for a machine with the minimum turning radius to move from the initial state to the next state. Its shortest path contains two segments with radius r k An arc and a straight line, and it must be one of the following path sets {RSL,RSR,LSR,LSL,RLR,LRL}, where R represents a clockwise arc, L represents a counterclockwise arc, and S represents a straight line segment.

4. The constrained optimization method for solving heterogeneous task allocation of multiple UAVs according to claim 3 is characterized by: Calculate the task execution time f of the task allocation scheme corresponding to the differentially evolved individuals T (X), satisfaction level μ(X), which is calculated as follows: f T (X)=max{f i / v},i=1,2,...M The above formula L ij It represents the path length that drone i flies when executing the jth target, which includes two parts: the path length from the previous target or the starting point to the target and the flight time required to execute the task. The flight time required to execute three heterogeneous targets is different. In order to improve the calculation reliability of the flight path length of drones in the task allocation scenario, the present invention abandons the Euclidean distance used in task allocation under normal circumstances, and calculates the path length value strictly according to the Dubins curve; in represents the flight distance of UAV i from the last target point or starting point to perform task j, It represents the flight distance required to perform mission j. Different target types have different flight distances. To perform a point target mission, the drone needs to fly around the target point. The flight distance of a line target needs to be determined by the length of the line target. The flight distance of a surface target will be determined by the length and width of the surface target. The basic reconnaissance model is shown in Figures 2, 3, and 4. The equality constraints and inequality constraints described in step 4 are described in detail as follows: Time window constraint: The exploration of point targets requires two drones to conduct reconnaissance. The arrival time of the two drones is t i , t j A certain time window must be satisfied. This inequality constraint is described as g1. g1=|t i -t j |-T win (8) Total time threshold constraint: There is also a limit on the completion time of the entire task. This constraint can be understood as the range constraint of the drone due to the limited fuel it carries. The inequality constraint is described as g2. The inequality constraint calculation is given by formula (9): g2=|F|-T sum (9) Redundancy constraint: From the task allocation scenario setting, it can be seen that the allocation between drones and targets must meet certain redundancy requirements. Each drone must be assigned to one target. At the same time, the number of drones assigned to each target must meet the requirements of its task type. For example, a line target requires two drones while a surface target only requires one drone. It is necessary to determine whether there is redundancy or lack in the number of drones required for each task. The equality constraint is described as h: where c ij is the decision variable representing the allocation relationship between UAV i and task j, c ij If it is 1, it means that UAV i performs the j task, N limj represents the number of drones required for task j, where point target N lim =2, Line target and area target N lim =1; M represents the total number of UAVs performing reconnaissance missions, and N represents the total number of reconnaissance missions.

5. The constrained optimization method for solving heterogeneous task allocation of multiple UAVs according to claim 1 is characterized by: Safety level α t Adaptive convergence rule for : To make the safety level α t It can slow down or accelerate convergence according to the evolution of the population, improve global and local search capabilities, and adjust the parameter β t Parameters and the relative feasible solution ratio k of the population flag Establish a connection, the relatively feasible solution refers to individual X i The satisfaction level μ(X i )≥α, the infeasible solution refers to the individual X i The satisfaction level μ(X i )<α, β t To control the safety level α t The key coefficient of the convergence rate, its value will directly determine α t Compared with the increment of the parent generation, the individual satisfaction level μ in the current population is greater than the safety level α t The proportion of the number of species in a population.

6. The constrained optimization method for solving heterogeneous task allocation of multiple UAVs according to claim 5 is characterized by: Improvements on adaptive crossover factor and selection operation in differential evolution algorithm: The crossover operation of the differential evolution algorithm described in step 8 is described in detail as follows: Cr(t)=C rmax -k flag ·(C rmax -C rmin ) (16) Among them C r (t) represents the crossover factor, r n The random number index that changes in each round ensures that there is at least one crossover operation in each round, otherwise the population will not change each time. rmax With C rmin represents the upper and lower bounds of the crossover factor, and C r The size of (t) determines the probability of crossover, and the selection of its size has a huge impact on the evolution rate of the population. In this invention, the feasibility ratio is continued to be used to determine the crossover constant k. flag Make adaptive changes. When the proportion of relatively feasible solutions in the population is large, it means that there are fewer violations of constraints in the current population, and the population may converge to the optimal solution from a larger possibility. Then, when performing crossover operations, performing smaller mutations and performing detailed dimensional searches are more conducive to fine convergence to the optimal solution. On the contrary, if the proportion of relatively feasible solutions in the current population is small, it means that the current population has serious violations of constraints. Then the population's evolutionary step size and the crossover probability of the mutation vector should tend to be larger, which is more conducive to jumping out of infeasible solutions. The selection operation of the differential evolution algorithm under the α level comparison described in step 10 is described in detail as follows: Assuming f(X1), f(X2) and μ(X1), μ(X2) are the fitness values ​​and satisfaction values ​​of X1 and X2 respectively, for any α∈[0,1], the α level comparison between individuals X1 and X2 is defined as follows: According to the α-level comparison, only when the satisfaction of two solutions is higher than the predefined value or the satisfaction is the same, the comparison between the two solutions will be made according to the objective function value. According to the α-level comparison method, the ordinary comparison is replaced by the α-level comparison, thereby converting the general constrained minimization problem into a minimization problem under the α-level comparison.

Images