A three-dimensional path planning method for UUV particle swarm based on multi-constraint optimization

Abstract

The invention aims to solve the problems of traditional particle swarm methods being sensitive to initialization conditions, difficult to adjust parameters, easy to fall into local optimal solutions and slow convergence speed in path planning, and discloses a UUV particle swarm three-dimensional path planning method with multi-constraint optimization, which specifically includes: on the basis of the traditional particle swarm algorithm, using pre-optimization to initialize some particles, and generating an initial population between a starting point and an end point straight line; introducing a simulated annealing criterion as a local search mechanism during optimization, and adjusting it according to the simulated annealing principle during position update; when simulated annealing cannot update an individual, the search space is increased through a mutation operation; when the particle updates its speed and position, its upper and lower bounds are restricted and corrected to avoid the particle from running off or stagnating due to excessively large or small speed, and to avoid the particle position from deviating from the map due to deviation; finally, some individuals are allowed to perform local search to further improve the convergence performance.

Description

Technical Field

[0001] The present invention relates to the field of unmanned underwater vehicle (UUV) path planning, and in particular to a multi-constraint optimized UUV particle swarm three-dimensional path planning method. Background Art

[0002] With the development of technology, the world pays more and more attention to the application of UUVs, including marine exploration: seabed topography mapping, marine biological surveys and marine resource exploration; marine environmental monitoring: UUVs can be used for real-time monitoring and data collection of marine environmental parameters to help scientists better understand the changes and evolution of marine ecosystems; seabed resource development: UUVs can be used to assist in the development and survey of seabed mineral resources, and provide important data support for the development of seabed resources; marine safety and rescue: UUVs can perform seabed searches and target positioning. In addition, UUVs can also be used for marine accident rescue, such as monitoring and cleaning of oil pollution incidents; marine scientific research: UUVs are used in marine science During the research, field observations and data collection in the fields of marine biology, geology, geophysics, etc. can be carried out, providing scientists with valuable marine research data; Underwater cultural heritage protection: UUVs are equipped with cameras, laser scanners and other equipment. UUVs can obtain high-definition images and three-dimensional models of underwater cultural relics, helping archaeologists to protect and study these precious cultural heritages. In addition to the above application areas, UUVs can also play a role in military aspects, marine engineering construction, marine energy development, marine traffic management, etc. The basis for the realization of these applications is the path planning of UUVs, which can safely reach the designated location within a limited time and consume limited energy.

[0003] Compared with the path planning in the two-dimensional plane, the path planning in the three-dimensional environment where UUV is located is more difficult. The main reason is the complexity of the underwater environment. Complex environmental conditions such as fish schools, islands and reefs, and ocean currents have brought considerable difficulties to the autonomous application of UUV. UUV itself also has many limitations, such as energy limitations, communication limitations, and compressive strength limitations, which limit the further development of UUV. Planning a reasonable route can save energy to the greatest extent, shorten the voyage, avoid danger, and improve the success rate of the mission. However, both traditional algorithms and emerging intelligent algorithms have certain limitations, such as low path search efficiency or high search cost. The quality of planning is directly related to the efficiency and cost of UUV execution of commands. At the same time, most algorithms are applicable to two-dimensional space, ignoring the impact of height on path planning, while the actual marine environment is extremely complex, and topography such as seamounts and ridges is inevitable in the process of executing tasks. This makes the traditional path planning method unable to meet the obstacle avoidance requirements in complex three-dimensional space environments. Therefore, path planning for UUV in a three-dimensional environment is inevitable.

[0004] Commonly used path planning methods include graph theory algorithms: Dijkstra algorithm, A* algorithm, etc.; search algorithms: depth-first search (DFS), breadth-first search (BFS), etc.; genetic algorithms; reinforcement learning algorithms: Q-learning, deep reinforcement learning, etc.; other heuristic algorithms: taboo search algorithm, particle swarm optimization algorithm, etc. Compared with other algorithms, particle swarm optimization has the advantages of being simple and easy to implement, having strong global search capabilities, and good robustness in a three-dimensional environment. However, the application of particle swarm optimization in the field of path planning has the following problems:

[0005] (1) Sensitive to parameter settings: The parameter settings in the particle swarm algorithm have a great impact on the performance and convergence speed of the algorithm. For example, parameters such as the number of particles, inertia weight, and acceleration factor need to be properly adjusted to achieve the best results. Inappropriate parameter selection may cause the algorithm to converge too slowly or fail to find a good solution.

[0006] (2) It is easy to fall into the local optimal solution: The particle swarm algorithm relies on information exchange and collaboration between particles. This collaboration mechanism may cause the algorithm to linger near the local optimal solution and it is difficult to escape the constraints of the local optimal solution. Especially in complex problem spaces, the particle swarm algorithm may not be able to fully search the entire solution space.

[0007] At the same time, compared with the present invention, the following methods have the following differences:

[0008] The difference between the improved method provided by the paper "A Particle Swarm Optimization Algorithm Based on Simplex Search" and the method described in the present invention is as follows: the paper does not improve and optimize the initialization of the particle swarm, and the initial population is not of high quality; the paper does not limit the speed and position of particles during the search, and the particles are prone to deviation and stagnation;

[0009] The improved method provided by the paper "Adaptive Viscous Particle Swarm Algorithm Based on Simulated Annealing Mechanism" is different from the method described in the present invention as follows: the paper does not improve the initialization of the particle swarm, and the initial population is not of high quality; the paper does not perform mutation after the simulated annealing criterion, and the particle exploration is not high; the paper does not perform limit correction on particle speed and position, and particles are prone to deviation;

[0010] The improved method provided by patent CN 117075471 A "A two-stage multi-target global path planning method for an intelligent unmanned ship" differs from the method described in the present invention as follows: the patent does not optimize the initial population of the particle swarm; the patent does not perform mutation particle exploration after the simulated annealing search, and the exploration is low; the patent does not weight multiple indicators, which is not in line with the actual situation.

[0011] In order to solve the problems that particle swarms are sensitive to parameter settings and prone to local optimal solutions, the present invention proposes a three-dimensional path planning method for UUV particle swarms with multi-constraint optimization, improves the population initialization, weights multiple constraints to better meet the actual UUV navigation to reduce energy consumption, and performs local search to help particles jump out of the local optimal solution. Summary of the invention

[0012] In order to solve the above-mentioned problems in path planning, the present invention proposes the following technical solutions: a multi-constraint optimized UUV particle swarm three-dimensional path planning method is designed. On the basis of the traditional particle swarm algorithm, the particle swarm initialization is improved, and some particles are initialized with high quality by pre-optimization, and then the initial population is generated between the starting point and the end point for further optimization, thereby improving the quality of the initial solution; in addition, in the optimization process, the simulated annealing criterion is introduced as a local search mechanism to avoid falling into the local optimal solution, and the position is adjusted according to the simulated annealing principle when the position is updated; when the simulated annealing cannot update the individual, the search space is increased through a mutation operation; and when the particle updates the speed and position, the upper and lower bounds of the speed and position are restricted and corrected to avoid the deviation or stagnation of the particle due to excessive or insufficient speed, and also to avoid the deviation of the particle position from the map; finally, some individuals are allowed to perform local search with a certain probability to further improve the convergence performance. Specifically, the following steps are included:

[0013] Step 1:

[0014] A three-dimensional seabed map is created under the UUV navigation environment, and the position coordinates of the starting point and the end point are set, and the offset of the moving direction and the dimension of the solution are generated.

[0015] Step 2:

[0016] Initialize the particle swarm parameters and initialize the population pop size , dimension dim, gene value pop max and pop min , total number of iterations, learning factor c 1 and c 2 , inertia weight w and speed range v max and v min .

[0017] Step 3:

[0018] Initialize the particle swarm population.

[0019] Step 4:

[0020] Calculate the initial fitness value, first linear interpolation and coordinate reshaping: linear interpolation of the input variables according to the range of the starting point and the end point to obtain the coordinates cor, the formula is as follows:

[0021] cor i,1 =cor i,1 *(end(1)-start(1))+start(1) (1)

[0022] cor i,2 =cor i,2 *size(2) (2)

[0023] cor i,3 =cor i,3 *size(3) (3)

[0024] Among them, cor i,j is the value of the ith point on the jth coordinate axis; start(1) is the position of the starting point of the map on the x-axis; end(1) is the position of the end point of the map on the x-axis; size(2) is the length of the map on the y-axis; size(3) is the length of the map on the z-axis.

[0025] Then normalization and mapping are performed: the first dimension x coordinate of the coordinate cor is normalized and mapped to the range of [0,1], linear mapping is performed according to the range of the starting point and the end point, the normalized x coordinate is restored to the actual coordinate value, and the second dimension y coordinate and the third dimension z coordinate of the coordinate cor are scaled according to the scaling factor.

[0026] Then perform B-spline interpolation: perform B-spline interpolation on the coordinates cor to obtain the coordinates cor of the interpolated point all , calculate the velocity v based on the coefficients of B-spline interpolation all , calculate the acceleration a based on the velocity difference and the interpolated breakpoints all , the calculation formula is as follows:

[0027] cor i =[x i ,y i ,z i ],x i ∈R,y i ∈R,z i ∈R (4)

[0028] cor={cor 1 ,cor 2 ,...,cor n} (5)

[0029] cor′=pointN(x′,cor) (6)

[0030]

[0031] cor″′=[start,cor i ″,end] (8)

[0032] cor all ″′,v all ″′,a all ″′=B-splineinterpolation(cor″′) (9)

[0033] cor all,i =[x all,i ,y all,i ,z all,i ],i=1,2,...,m (10)

[0034] Among them, x′ is the parameter vector; cor is the original coordinate point; cor′ is the adjusted coordinate point; cor i ″ is the coordinate point after scaling the coordinate point; cor″′ is the coordinate point after adding the start point and the end point at the beginning and the end of the coordinate point; cor all ″′ is the position information of all points obtained by interpolation; v all ″′ is the velocity information of all points obtained by interpolation; a all ″′ is the acceleration information of all points obtained by interpolation; cor all are the position coordinates of all points after interpolation.

[0035] Then calculate the curvature: Calculate the curvature k based on the acceleration and velocity, the formula is as follows:

[0036]

[0037] Among them, k i is the curvature of the i-th point on the curve; v i is the velocity at the i-th point on the curve after interpolation; a i is the acceleration at the i-th point on the curve after interpolation.

[0038] Continue with constraint processing: set the coordinates cor all The z coordinate of the site map is interpolated in two dimensions to obtain the interpolated height z interp , the penalty term of the statistical constraints, including that the insertion height cannot be lower than the site map and the coordinates cannot exceed the boundary, is as follows:

[0039]

[0040] Where m is the number of points on the interpolated path; z interpis the height information of the interpolation point on the map; restraint(1) is the non-collision constraint; restraint(2) is the non-exceeding boundary constraint; restraint(3) is the non-exceeding site size constraint; all constraints are combined into the total constraint res, which is the sum of constraints restraint(1), restraint(2) and restraint(3).

[0041] The objective function calculation after optimizing the four objectives is as follows: calculate the total distance dis, that is, the sum of the Euclidean distances between all points, calculate the total height height, that is, the sum of the z coordinates of all points, and calculate the curvature mean k , and the penalty term res of the constraint condition, combine the above calculation results, and use the weighted coefficient to get the objective function value fitness, the formula is as follows:

[0042]

[0043] fitness=0.5*dis+0.2*height+0.3*mean k +res*10 3 (18)

[0044] Where dis is the total length of the path; height is the total height of the UUV during navigation; mean k is the average curvature of the path; res is the sum of all constraints; fitness is the function value of multi-objective optimization.

[0045] Step 5:

[0046] Find the initial extreme value, fitness best is the optimal fitness value, indicating the best fitness value currently found, and initializes the best fitness best is positive infinity, index best is the optimal individual index, indicating the position of the individual corresponding to the optimal fitness value in the population, and initializing index best is 0, traverse each individual in the population, calculate its fitness value, and compare it with the current optimal fitness value fitness best Compare, if the individual's fitness value is less than the current optimal fitness value fitness best , then update fitness best is the fitness value of the individual, and the index of the individual is recorded as index best , and finally according to index best Find the corresponding z bestis the extreme position of the population, indicating the position of the individual with the best fitness value in the entire population, and the corresponding It is the extreme fitness value of the group, which means the fitness value of the individual with the best fitness value in the entire population.

[0047] Step 6:

[0048] Set the simulated annealing temperature drop process, where the starting value is 200 and the ending value is 0, indicating that the temperature will drop from the starting value to 0. The formula is as follows:

[0049]

[0050] Among them, temp i is the i-th temperature value; iterations is the total number of temperature values ​​that need to be generated.

[0051] Since the initial temperature is 200, when i=0, temp 0 =200; when i=iterations-1, temp iterations-1 =0; when i takes any value, temp i It will change between 200 and 0. As the number of iterations increases, the value of i keeps changing, thus achieving a continuous decrease in temperature.

[0052] Step 7:

[0053] Record the historical optimal solution, define and initialize the vector History to store the historical optimal solution of each iteration, and convert the optimal solution of the current iteration to Add to the History vector.

[0054] Step 8:

[0055] Update speed, the formula is as follows:

[0056] V i,j =wV i,j +c 1 r 1 (pbest i,j -x i,j )+c 2 r 2 (gbest j -x i,j ) (20)

[0057] Among them, V i,j is the velocity of the ith particle in the jth dimension; w is the inertia weight; c 1 and c 2 is the acceleration constant; r 1 and r 2is a random number between [0,1]; pbest i,j is the position of the extreme value of the i-th individual in the j-th dimension; gbest j is the position of the group extreme value in the jth dimension; x i,j is the position of the ith particle in the jth dimension.

[0058] In order to ensure that the speed range is within v min and v max The speed that exceeds the range is corrected, and the formula is as follows:

[0059]

[0060] Among them, v j,i is the velocity of particle j in the i-th dimension; v min is the lower bound of velocity; v max is the upper bound of the speed.

[0061] Step 9:

[0062] According to the updated speed, update the particle's position. The formula is as follows:

[0063] x j,i =x j,i +v j,i (twenty two)

[0064] Among them, x j,i is the position of particle j in the i-th dimension; v j,i is the velocity of particle j in the i-th dimension.

[0065] Correct the new position to ensure that it is within the specified range. The formula is as follows:

[0066]

[0067] Among them, x j,i is the position of particle j in the i-th dimension; ub i is the upper bound on the i-th dimension; lb i is the lower bound on the ith dimension.

[0068] Step 10:

[0069] Update and mutate the population.

[0070] Step 11:

[0071] The random number is used to determine whether to perform local search optimization. If the random number is less than 0.5, the process proceeds to step 12; otherwise, the process proceeds to step 13.

[0072] Step 12:

[0073] Conduct further local search optimization.

[0074] Step 13:

[0075] Record the optimal value of each iteration.

[0076] Step 14:

[0077] Determine whether the number of iterations of the entire population has reached 100. If so, proceed to step 15; otherwise, return to step 8.

[0078] Step 15:

[0079] Output the optimal value.

[0080] The present invention has the following beneficial effects:

[0081] (1) The method of the present invention optimizes the UUV navigation path length, navigation heave height and navigation turning range by linear weighting. The path generated by the present invention has smaller fluctuations in heave height between -0.15 and +0.11, while the path generated by the traditional particle swarm algorithm fluctuates between -0.31 and +0.62, indicating that the method of the present invention is more stable in navigation altitude, more conducive to reducing the frequent surfacing and diving of UUV, and thus reducing the energy consumption of UUV. The path generated by the present invention has smaller fluctuations in turning range, which is reflected by the curvature change, with a maximum value of 0.31, while the maximum curvature change of the path generated by the traditional particle swarm algorithm is 3.78, which is 91.7% less than the maximum curvature of the traditional particle swarm algorithm, indicating that the method of the present invention has smaller fluctuations in turning range, and the generated path is smoother and more meets the navigation angle constraint.

[0082] (2) Based on the traditional particle swarm algorithm initialization population, the method of the present invention uses pre-optimization to perform high-quality initialization on some particles, and initializes a population between the starting point and the end point straight line, which increases the diversity of the population and helps to better explore the search space. Since the initialization of the present invention is improved, the initial function value of the present invention is 3.16e4 and the initial function value of the traditional particle swarm algorithm is 1.02e5, which is 69% less than the initial function value of the traditional particle swarm algorithm, which shows that the initial population of the particle swarm of the method of the present invention is better.

[0083] (3) The method of the present invention introduces a simulated annealing criterion on the basis of traditional particle swarm. When updating individual solutions, a certain probability is adopted to accept a worse solution according to the change of fitness value and current temperature to avoid falling into a local optimal solution. When simulated annealing cannot update the individual, a mutation operation is performed to increase the search space, thereby improving the global search capability of the algorithm and increasing the diversity and exploratory nature of the algorithm.

[0084] (4) When updating the speed and position, the method of the present invention limits the upper and lower bounds of the speed and position, and makes corrections to avoid the speed being too large or too small, which causes the particles to deviate or stagnate, and also avoids the deviation of the particle position from the map. In addition, a part of individuals are selected in the group for local search, which can converge to the local optimal solution and particle swarm more quickly. The traditional particle swarm algorithm has an obvious "circling" path, which makes the path length obviously longer, while the path generated by the method of the present invention is better in smoothness and rationality. BRIEF DESCRIPTION OF THE DRAWINGS

[0085] Figure 1 This is the main flow chart of a multi-constraint optimized UUV particle swarm three-dimensional path planning method;

[0086] Figure 2 Flowchart of particle swarm initialization in the main flow chart;

[0087] Figure 3 A flow chart for updating and mutating the population in the main flow chart;

[0088] Figure 4 Flowchart for local search optimization in the main flow chart;

[0089] Figure 5 It is a front view of the formed path diagram;

[0090] Figure 6 A side view of the formed path diagram;

[0091] Figure 7 A top view of the path is formed;

[0092] Figure 8 It is the change diagram of the convergence curve;

[0093] Fig. 9 is the height change diagram;

[0094] Fig.10 This is the curvature variation diagram. DETAILED DESCRIPTION

[0095] Figure 1 The main flow chart of the multi-constraint optimized UUV particle swarm three-dimensional path planning method of the present invention includes the following steps:

[0096] Step 1:

[0097] Create a grid coordinate matrix, obtain the grid data including the vector of the number of rows, columns and height data, merge the size information with the vector consisting of the upper integer value of the maximum height, set the starting and ending coordinates, the starting coordinates [1,17,4], the ending coordinates [92,78,4], set the number of nodes to be inserted to 10, create an empty matrix direct as the moving direction matrix, this matrix is ​​used to store the direction of movement in three-dimensional space, each row represents the offset in one direction, use a triple loop to traverse the moving direction, generate the offset of the moving direction by taking the value of -1, 0 or 1 in each dimension, and add it to the matrix direct, delete all rows in direct where all elements are 0, store the result, set the dimension of the solution to the number of nodes multiplied by 3, and now the three-dimensional seabed map is established and the map parameters and the starting and ending coordinates are set.

[0098] Step 2:

[0099] Initialize the particle swarm parameters and set the population size pop size =80 is used to determine the size of the particle group, that is, the number of particles; dim=30 is set to determine the dimension of each particle, that is, the dimension of the solution to the problem; the maximum value of the gene pop is set max =1 and the minimum value of gene pop min = 0 is used to determine the value range of the solution in each dimension; set the total number of iterations to 100; set the learning factor c 1 =1.457 and c 2 =1.457 is used to adjust the particle speed update, affecting the degree of the particle's trade-off between individual experience and group experience; setting the inertia weight w=0.7 is used to balance the impact of the particle's individual experience and group experience on the speed, determining the impact of the particle's previous speed on the current speed; setting the speed range v min = -0.5 and v max =0.5 is used to limit the range of particle speed changes to ensure that the speed does not exceed the set range.

[0100] Step 3:

[0101] The particle swarm population initialization in the claims is described in detail step by step, and the flow chart is shown in Figure 2 :

[0102] Step 3.1:

[0103] Defines a size of pop size*dim's zero matrix pop is the population, a two-dimensional array that stores the generated solution vectors and is used to save the initialized population. Next, the parameter NN is assigned a value of 2, where NN represents the number of elite individuals selected. The value 2 indicates that the first two solution vectors in the population array pop already exist. An empty vector a=[] is initialized to store the generated solutions, and then a vector is generated from pop in each loop through three cycles. min to pop max The generated partial solutions are appended to the vector a by horizontal splicing, so that the partial solutions generated in each cycle are added to the end of a, and the complete solution is gradually constructed. Finally, the complete solution a generated by three cycles is assigned to the first individual of the population, and the solution is used as the first individual of the initial population.

[0104] Step 3.2:

[0105] Initialize the open set open to store the nodes to be traversed, add the starting point to open, initialize f as an array to record the distance information of each node, each element represents three values ​​of the node: the distance to the starting point, the distance to the end point and the total distance. For the starting point, set its distance to the starting point to 0, set the estimated distance to the end point to the Euclidean distance from the starting point to the end point, and the total distance to the Euclidean distance from the starting point to the end point. Initialize the closed set close set , used to store the nodes that have been traversed. By adding nodes to the closed set, we can avoid repeatedly traversing the same nodes. Initialize the parent set parent to store their parent node information. After the search is completed, the shortest path can be determined by backtracking the parent node information.

[0106] Step 3.3:

[0107] Select the node with the smallest f value from the open set open as the current node last is the node currently being traversed, and add the current node to the closed set close set middle.

[0108] Step 3.4:

[0109] Determine whether the current node last is the end point. If so, proceed to step 3.13; otherwise, proceed to step 3.5.

[0110] Step 3.5:

[0111] Traverse the nodes around the current node last and determine whether the node is feasible based on the conditions, including whether it exceeds the lower boundary, whether it exceeds the upper boundary, whether it is an obstacle, whether it is in the closed set close set If the node is feasible, go to step 3.6; otherwise, set the function value of the node to infinity and continue traversing the surrounding nodes.

[0112] Step 3.6:

[0113] Determine whether the node is the end point. If so, proceed to step 3.8; otherwise, calculate the function value f of the node. 1 and f 2 and f 12 , the formula is as follows:

[0114] f 12 =f 1 +f 2 (1)

[0115] Among them, f 1 is the distance from the starting point to the current node; f 2 The estimated distance from the current node to the end point.

[0116] And proceed to step 3.7.

[0117] Step 3.7:

[0118] Determine whether the node is open. If so, proceed to step 3.9; otherwise, proceed to step 3.8.

[0119] Step 3.8:

[0120] Add the node to the open set open and go to step 3.10.

[0121] Step 3.9:

[0122] Determine whether the function value of the node is better. If so, go to step 3.10; otherwise, go to step 3.11.

[0123] Step 3.10:

[0124] Update the corresponding function value f and parent node parent information and go to step 3.11.

[0125] Step 3.11:

[0126] Delete the current node from the open set open and go to step 3.12.

[0127] Step 3.12:

[0128] Determine whether the traversal loop is completed. If so, go to step 3.13; otherwise, return to step 3.3.

[0129] Step 3.13:

[0130] Start backtracking: Add the endpoint to the path: Add the endpoint to the path route, which is used to store the final shortest path; Find the parent node: Find the endpoint in the closed set close set, and get its parent node information; add the parent node to the path: add the parent node information to the path route; until the parent node is the starting point, the backtracking ends, the path route is reversed, the shortest forward path is obtained, and the final path route is output as a population of the initial particle swarm.

[0131] Step 3.14:

[0132] The remaining solution vectors are generated through a loop starting from NN+1. In each loop, a random number vector of size 1*dim is generated in the interval [0,1). Next, the element values ​​of the random number vector are mapped to the specified range [pop min ,pop max ], the purpose of this linear transformation is to limit the element values ​​of the randomly generated solution vector to pop min and pop max Finally, the generated solution vector is assigned to the population array pop as the remaining population of the particle swarm.

[0133] Step 4:

[0134] Calculate the initial fitness value, first linear interpolation and coordinate reshaping: linear interpolation of the input variables according to the range of the starting point and the end point to obtain the coordinates cor, the formula is as follows:

[0135] cor i,1 =cor i,1 *(end(1)-start(1))+start(1) (2)

[0136] cor i,2 =cor i,2 *size(2) (3)

[0137] cor i,3 =cor i,3 *size(3) (4)

[0138] Among them, cor i,j is the value of the ith point on the jth coordinate axis; start(1) is the position of the starting point of the map on the x-axis; end(1) is the position of the end point of the map on the x-axis; size(2) is the length of the map on the y-axis; size(3) is the length of the map on the z-axis.

[0139] Then normalization and mapping are performed: the first dimension x coordinate of the coordinate cor is normalized and mapped to the range of [0,1], linear mapping is performed according to the range of the starting point and the end point, the normalized x coordinate is restored to the actual coordinate value, and the second dimension y coordinate and the third dimension z coordinate of the coordinate cor are scaled according to the scaling factor.

[0140] Then perform B-spline interpolation: perform B-spline interpolation on the coordinates cor to obtain the coordinates cor of the interpolated point all , calculate the velocity v based on the coefficients of B-spline interpolation all , calculate the acceleration a based on the velocity difference and the interpolated breakpoints all , the calculation formula is as follows:

[0141] cor i =[x i ,y i ,z i ],x i ∈R,y i ∈R,z i ∈R (5)

[0142] cor={cor 1 ,cor 2 ,...,cor n} (6)

[0143] cor′=pointN(x′,cor) (7)

[0144]

[0145] cor″′=[start,cor i ″,end] (9)

[0146] cor all ″′,v all ″′,a all ″′=B-splineinterpolation(cor″′) (10)

[0147] cor all,i =[x all,i ,y all,i ,z all,i ],i=1,2,...,m (11)

[0148] Among them, x′ is the parameter vector; cor is the original coordinate point; cor′ is the adjusted coordinate point; cor i ″ is the coordinate point after scaling the coordinate point; cor″′ is the coordinate point after adding the start point and the end point at the beginning and the end of the coordinate point; cor all″′ is the position information of all points obtained by interpolation; v all ″′ is the velocity information of all points obtained by interpolation; a all ″′ is the acceleration information of all points obtained by interpolation; cor all are the position coordinates of all points after interpolation.

[0149] Then calculate the curvature: Calculate the curvature k based on the acceleration and velocity, the formula is as follows:

[0150]

[0151] Among them, k i is the curvature of the i-th point on the curve; v i is the velocity at the i-th point on the curve after interpolation; a i is the acceleration at the i-th point on the curve after interpolation.

[0152] Continue with constraint processing: set the coordinates cor all The z coordinate of the site map is interpolated in two dimensions to obtain the interpolated height z interp , the penalty term of the statistical constraints, including that the insertion height cannot be lower than the site map and the coordinates cannot exceed the boundary, is as follows:

[0153]

[0154] Where m is the number of points on the interpolated path; z interp is the height information of the interpolation point on the map; restraint(1) is the non-collision constraint; restraint(2) is the non-exceeding boundary constraint; restraint(3) is the non-exceeding site size constraint; all constraints are combined into the total constraint res, which is the sum of constraints restraint(1), restraint(2) and restraint(3).

[0155] The objective function calculation after optimizing the four objectives is as follows: calculate the total distance dis, that is, the sum of the Euclidean distances between all points, calculate the total height height, that is, the sum of the z coordinates of all points, and calculate the curvature mean k , and the penalty term res of the constraint condition, combine the above calculation results, and use the weighted coefficient to get the objective function value fitness, the formula is as follows:

[0156]

[0157] fitness=0.5*dis+0.2*height+0.3*mean k +res*10 3 (19)

[0158] Where dis is the total length of the path; height is the total height of the UUV during navigation; mean k is the average curvature of the path; res is the sum of all constraints; fitness is the function value of multi-objective optimization.

[0159] Step 5:

[0160] Find the initial extreme value, fitness best is the optimal fitness value, indicating the best fitness value currently found, and initializes the best fitness best is positive infinity, index best is the optimal individual index, indicating the position of the individual corresponding to the optimal fitness value in the population, and initializing index best is 0, traverse each individual in the population, calculate its fitness value, and compare it with the current optimal fitness value fitness best Compare, if the individual's fitness value is less than the current optimal fitness value fitness best , then update fitness best is the fitness value of the individual, and the index of the individual is recorded as index best , and finally according to index best Find the corresponding z best is the group extreme value position and the corresponding It is the extreme fitness value of the group, which means the fitness value of the individual with the best fitness value in the entire population.

[0161] Step 6:

[0162] Set the simulated annealing temperature drop process, where the starting value is 200 and the ending value is 0, indicating that the temperature will drop from the starting value to 0. The formula is as follows:

[0163]

[0164] Among them, temp i is the i-th temperature value; iterations is the total number of temperature values ​​that need to be generated.

[0165] Since the initial temperature is 200, when i=0, temp 0 =200; when i=iterations-1, temp iterations-1 =0; when i takes any value, temp i It will change between 200 and 0. As the number of iterations increases, the value of i keeps changing, thus achieving a continuous decrease in temperature.

[0166] Step 7:

[0167] Record the historical optimal solution, define and initialize the vector History to store the historical optimal solution of each iteration, and convert the optimal solution of the current iteration to Add to the History vector.

[0168] Step 8:

[0169] Update speed, the formula is as follows:

[0170] V i,j =wV i,j +c 1 r 1 (pbest i,j -x i,j )+c 2 r 2 (gbest j -x i,j ) (twenty one)

[0171] Among them, V i,j is the velocity of the ith particle in the jth dimension; w is the inertia weight; c 1 and c 2 is the acceleration constant; r 1 and r 2 is a random number between [0,1]; pbest i,j is the position of the extreme value of the i-th individual in the j-th dimension; gbest j is the position of the group extreme value in the jth dimension; x i,j is the position of the ith particle in the jth dimension.

[0172] In order to ensure that the speed range is within v min and v max The speed that exceeds the range is corrected, and the formula is as follows:

[0173]

[0174] Among them, v j,i is the velocity of particle j in the i-th dimension; v min is the lower bound of velocity; v max is the upper bound of the speed.

[0175] Step 9:

[0176] According to the updated speed, update the particle's position. The formula is as follows:

[0177] x j,i =x j,i +v j,i (twenty three)

[0178] Among them, x j,i is the position of particle j in the i-th dimension; v j,i is the velocity of particle j in the i-th dimension.

[0179] Correct the new position to ensure that it is within the specified range. The formula is as follows:

[0180]

[0181] Among them, x j,i is the position of particle j in the i-th dimension; ub i is the upper bound on the i-th dimension; lb i is the lower bound on the ith dimension.

[0182] Step 10:

[0183] A detailed step-by-step description of the population update mutation in the claims is provided, and the flow chart is shown in Figure 3 :

[0184] Step 10.1:

[0185] Calculate the fitness value of the current individual j And the fitness value fit1 of the new individual, calculate the change in fitness value, the formula is as follows:

[0186] df = fit1-fitness j (25)

[0187] Among them, df is the fitness value difference, that is, the difference between the fitness value of the new individual and the fitness value of the old individual; fit1 is the fitness value of the new individual, which is calculated by the objective function; fitness j is the fitness value of the jth individual, calculated by the objective function.

[0188] Step 10.2:

[0189] Determine whether df is less than 0. If so, proceed to step 10.3; otherwise, proceed to step 10.4.

[0190] Step 10.3:

[0191] This means that the fitness value of the new individual is better, so the new individual is accepted, and the position of individual j is updated to the new individual position x 1 , and update the fitness value of individual j to the new fitness value fit1, which can ensure that the individuals in the population always maintain excellent fitness values, and update the individual position and fitness value. The formula is as follows:

[0192]

[0193] Among them, pop j is the jth individual in the population; x 1 is the newly generated individual position, that is, the updated individual position; fitness j is the fitness value of the jth individual; fit1 is the fitness value of the new individual.

[0194] Then proceed to step 10.6.

[0195] Step 10.4:

[0196] If the condition of step 10.2 is not met, it means that the fitness value of the new individual is worse than the fitness value of the original individual. Generate a random number r to determine whether it satisfies the simulated annealing formula, that is, compare r with The specific formula is as follows:

[0197]

[0198] Among them, r is a random number with a value range of [0,1); df is the fitness value difference, that is, the difference between the new individual fitness value and the old individual fitness value; tempiter is the simulated annealing temperature corresponding to the current iteration number; dr is the value of r and The difference.

[0199] If dr is less than 0, it means that the annealing formula is satisfied, and then return to step 10.3; otherwise, proceed to step 10.5.

[0200] Step 10.5:

[0201] If the annealing formula is not satisfied, the individual is mutated. First, the mutation rate mp is defined as follows:

[0202]

[0203] Among them, pop size is the population size, i.e. the number of individuals.

[0204] Set the distribution index etam=10 to control the variation intensity.

[0205] The formula for setting the scale is as follows:

[0206] N=size(pop,1) (29)

[0207] D=size(pop,2) (30)

[0208] Among them, size(pop,1) is the number of rows of the population matrix, that is, the number of individuals; size(pop,2) is the number of columns of the population matrix, that is, the number of dimensions of an individual; N is the number of individuals in the population; D is the number of dimensions representing each individual.

[0209] Normalization is performed, and the formula is as follows:

[0210]

[0211] Among them, δ 1 is the normalized value of the individual relative to the minimum boundary; δ 2 is the normalized value of the individual relative to the maximum boundary; x 2 Represents the value of a dimension in an individual; pop min Represents the minimum boundary value; pop max Represents the maximum boundary value.

[0212] δ 1 The value range of is 0 to 1, reflecting the position of the individual relative to the boundary on this dimension. 1 When δ is 0, the individual is equal to the minimum boundary value in this dimension; 1 When δ is 1, the individual is equal to the maximum boundary value in this dimension; 1 When it is 0.5, the individual is located in the middle of the boundary on this dimension.

[0213] For the selection of mutation position: pos is a vector of length D, each element of which is a random number uniformly distributed between 0 and 1. It is used to determine which dimensions will be mutated. The formula is as follows:

[0214]

[0215] Among them, rand(1,D) is a random vector of length D, which represents a random number uniformly sampled from the interval [0,1].

[0216] For generating random numbers: μ is a vector of length D, each element of which is a random number uniformly distributed between 0 and 1. It is used to calculate the variation amplitude. The formula is as follows:

[0217] μ=rand(1,D) (34)

[0218] Among them, rand(1,D) is a random vector of length D, which represents a random number uniformly sampled from the interval [0,1].

[0219] Then calculate the variation range of some dimensions in the population, for x 2≤0.5, first exclude the dimensions with random number μ>0.5, and determine the index of the dimension that meets the condition. The formula is as follows:

[0220]

[0221] Among them, pos i is the i-th element of the mutation position vector pos; μ i is the i-th element of the random vector μ.

[0222] Next, select the elements that meet the conditions from the random vector μ to obtain the random vector μ 1 , the formula is as follows:

[0223] μ 1 =μ(index) (36)

[0224] Among them, index is the index of the mutation position that meets the requirements.

[0225] The formula for calculating the variation range δ is as follows:

[0226]

[0227] Among them, δ 1 is a proportional factor vector used to calculate the variation amplitude; η m is the distribution index, which represents the parameter that controls the amplitude variation.

[0228] Add the variation amplitude δ to the corresponding dimension of individual x to obtain the mutated individual x′, the formula is as follows:

[0229] x′=x+δ·(pop max -pop min ) (38)

[0230] For x 2 > 0.5, first exclude the dimensions with random number μ≤0.5, and determine the index of the dimension that meets the condition. The formula is as follows:

[0231]

[0232] Among them, pos i is the i-th element of the mutation position vector pos; μ i is the i-th element of the random vector μ.

[0233] Next, select the elements that meet the conditions from the random vector μ to obtain the random vector μ 1 , the formula is as follows:

[0234] μ 1 =μ(index) (40)

[0235] Among them, index is the index of the mutation position that meets the requirements.

[0236] The formula for calculating the variation range δ is as follows:

[0237]

[0238] Among them, δ 2 is a proportional factor vector used to calculate the variation amplitude; η m is the distribution index, which represents the parameter that controls the amplitude variation.

[0239] Add the variation amplitude δ to the corresponding dimension of individual x to obtain the mutated individual x′, the formula is as follows:

[0240] x′=x+δ·(pop max -pop min ) (42)

[0241] Finally, the position of the individuals generated by the above mutation is corrected, and the individuals whose value is less than the lower limit pop min Correct the elements: first find the position of the elements that are less than the lower limit, and correct the elements that exceed the lower limit to the randomly generated ones in the pop min and pop max The value between is as follows:

[0242] ret(index)=pop min +0.01·(pop max -pop min )·rand(1,Σ(index)) (43)

[0243] Among them, ∑(index) is the number of true values ​​in the statistical vector index, which is used to determine the number of generated random numbers; ret is the corrected individual position vector.

[0244] For values ​​greater than the upper limit pop max Correct the elements: find the position of the elements that are greater than the upper limit, and correct the elements that exceed the upper limit to the randomly generated ones in the pop min and pop max The value between is as follows:

[0245] ret(index)=pop max -0.01·(pop max -pop min )·rand(1,∑(index)) (44)

[0246] Among them, Σ(index) is the number of true values ​​in the statistical vector index, which is used to determine the number of generated random numbers; ret is the corrected individual position vector.

[0247] Then proceed to step 10.6;

[0248] Step 10.6:

[0249] Update the individual optimal solution and the global optimal solution. For each particle, determine whether its current fitness value is less than the individual optimal solution fitness value of the particle. If so, update the individual optimal solution position of the particle to the current position, and update the individual optimal solution fitness value to the current fitness value. For each particle, determine whether its current fitness value is less than the global optimal solution fitness value. If so, update the global optimal solution position to the current position, and update the global optimal solution fitness value to the current fitness value. Sort the fitness values ​​of the current population. After sorting, the individual with the lowest fitness is at the beginning of the array, and the individual with the highest fitness is at the end of the array, and then enter step 11.

[0250] Step 11:

[0251] The random number is used to determine whether to perform local search optimization. If the random number is less than 0.5, the process goes to step 12; otherwise, the process goes to step 13.

[0252] Step 12:

[0253] A detailed step-by-step description of the local search optimization in the claims is provided, and the flowchart is shown in Figure 4 :

[0254] Step 12.1:

[0255] Select the best and worst ranked first N individuals (N is 1 / 10 of the population, rounded up) and merge them into a new population. Select the indexes of the first N and last N individuals from the sorted indexes and save the indexes of these individuals in ch index middle.

[0256] Step 12.2:

[0257] Initialize the number of iterations to 20, set the expansion factor α to an array of length 20, and the value of each element is 1.5. Set the scaling factor β to an array of length 20, and the value of each element is 0.5. In each iteration, the expansion operation and the scaling operation will use the same expansion factor and scaling factor. The expansion operation expands the position of the worst individual outward. The expansion factor α determines the degree of expansion. α is set to 1.5, indicating that each expansion operation will increase the position of the worst individual by 1.5 times. The scaling operation is to generate a new point and replace it with the worst individual. The scaling factor β determines the position of the generated point. β is set to 0.5, indicating that the generated point is between the best individual and the worst individual, and random perturbations are used to increase the diversity of the search.

[0258] Step 12.3:

[0259] Arrange the fitness values ​​of the new population in descending order of fitness value.

[0260] Step 12.4:

[0261] Calculate the mean x4 of individuals other than the worst individual as the new reference point. The formula is as follows:

[0262]

[0263] Among them, x4 is the average value of all individuals in the population except the worst individual in each dimension; N is the number of individuals in the population; pop' i is the value of the i-th individual in the population.

[0264] Calculate the reflection point x5 of the worst individual about x4. The formula for calculating x5 is as follows:

[0265] x5 j =2·x4 j -pop' 1,j (46)

[0266] Where j is the index of the dimension, ranging from 1 to D; x5 j is the value of the individual in the newly generated population in the jth dimension; x4 j is the value of the individual in the original population in the jth dimension; pop' 1,j is the value of the first individual in the population in all dimensions.

[0267] And correct x5 according to the correction formula in step 10.5, and calculate the fitness fit5 according to the new reflection point x5 and the formula in step 4.

[0268] Step 12.5:

[0269] If the fitness value of x5 is better than the fitness value of the current worst individual, then the expansion operation is performed: the expansion point x6 is calculated, and the formula is as follows:

[0270] x6=x4+α·(x4-pop' 1,j ) (47)

[0271] Among them, x6 is the expanded point, which is a vector; x4 is the mean point, which is a vector, representing the mean of all individuals in the population except the worst individual; pop' 1,j is the value of the first individual in the population in all dimensions.

[0272] And correct x6 according to the correction formula in step 10.5, calculate the fitness fit6 according to the new reflection point x6 and the formula in step 4, if the fitness value of the expansion point x6 is better than x5, replace x6 with the worst individual; otherwise, replace x5 with the worst individual.

[0273] Step 12.6:

[0274] If the fitness value of x5 is between the current best individual and the second to last individual, a scaling operation is performed: the scaling point x7 is calculated using the following formula:

[0275] x7=x4+β·rand(1,dim)⊙(x4-pop' 1,j ) (48)

[0276] Among them, x7 is a vector, representing the expanded point; x4 is a vector, representing the mean point; rand(1,dim) is a random vector with dim dimension, where dim represents the dimension of the vector; pop' 1,j is the value of the first individual in the population in all dimensions.

[0277] And correct x7 according to the correction formula in step 10.5, replace x7 with the worst individual, and calculate the fitness fit7 based on the new reflection point x7 and the formula in step 4.

[0278] Step 12.7:

[0279] If the fitness value of x5 is worse than the worst individual, then select a point x8 and perform an update operation to calculate x8. The formula is as follows:

[0280] x8=x4+β·rand(1,dim)⊙(pop' 1,j -x4) (49)

[0281] Among them, x8 is a vector; x4 is a vector, representing the mean point; rand(1,dim) is a random vector with dim dimension, where dim represents the dimension of the vector; pop' 1,j is the value of the first individual in the population in all dimensions.

[0282] And correct x8 according to the correction formula in step 10.5, calculate the fitness fit8 according to the new reflection point x8 and the formula in step 4, if the fitness value of x8 is better than the worst individual, replace x8 with the worst individual, otherwise, generate a point x9 between the best individual and the worst individual, and replace it with the worst individual.

[0283] Step 12.8:

[0284] Sort the entire population in ascending order according to the fitness value, and determine whether the number of simplex iterations has reached 20. If so, return the optimal individual as z best , that is, the optimization result; otherwise, return to step 12.3.

[0285] Step 13:

[0286] The best value of each iteration is recorded.

[0287] Step 14:

[0288] Determine whether the number of iterations of the entire population has reached 100. If so, proceed to step 15; otherwise, return to step 8.

[0289] Step 15:

[0290] Output the optimal value.

[0291] The method of the present invention is simulated by Matlab software. Figure 5 , Figure 6 and Figure 7 It is a path simulation diagram generated by the UUV from the starting point S[1,17,4] to the end point F[92,78,4] under the three-dimensional seabed map, where the units of the x-axis, y-axis and z-axis are km, the solid line is the path generated by the present invention, and the dotted line is the path generated by the traditional particle swarm algorithm. Figure 5 To generate the front view of the path, Figure 6 To generate a side view of the path, Figure 7 To generate a top view of the path, it can be concluded from the figure that the path generated by the traditional particle swarm has an obvious "circling" path, while the path generated by the method of the present invention is more reasonable, shorter and smoother.

[0292] Figure 8is a graph of the convergence curve change. It can be concluded from the graph that the initial objective function value of the present invention after the improved population initialization is 31626.8, while the initial objective function value of the traditional particle swarm algorithm after initialization is as high as 102831.0, which is 3.2 times that of the present invention; the objective function value of the method of the present invention is reduced from the initial 31626.8 to 468.0 at the 4th iteration, while the traditional particle swarm algorithm reduces the objective function value from the initial 102831.0 to 1892.7 at the 9th iteration, and the objective function value of the traditional particle swarm algorithm is 478.2 at the 64th iteration, which is close to the present invention. The effect of the fourth iteration is shown in FIG. 1 . The efficiency of the objective function of the present invention in reducing the objective function is 16 times that of the traditional particle swarm algorithm. When the present invention iterates 73 times, the objective function value begins to be basically stable at 267.0, while the objective function value of the traditional particle swarm algorithm tends to be stable at 379.7 only when it iterates 98 times. When the present invention iterates 96 times, the objective function value reaches a minimum of 262.0, while when the traditional particle swarm algorithm iterates 98 times, the objective function value reaches a minimum of 379.7. The minimum objective function value generated by the traditional particle swarm algorithm is 1.4 times the minimum objective function value generated by the present invention.

[0293] Fig. 9 The path generated by the present invention has smaller fluctuations in the heave height between -0.15km and +0.11km, while the path generated by the traditional particle swarm algorithm fluctuates between -0.31km and +0.62km, which is 72.1% less than the traditional particle swarm algorithm. This shows that the method of the present invention is more stable in the navigation altitude, which is more conducive to reducing the frequent surfacing and diving of UUV, thereby reducing the energy consumption of UUV. Fig.10 The curvature change diagram shows that the path generated by the present invention has a smaller fluctuation in the turning amplitude, and its maximum value is 0.31m as shown by the curvature change. -1 , while the maximum change of path curvature produced by traditional particle swarm algorithm is 3.78m -1 , while the maximum change of path curvature produced by traditional particle swarm algorithm is 3.78m -1 , compared with the traditional particle swarm algorithm, the maximum curvature is reduced by 91.7%, indicating that the method of the present invention has smaller fluctuations in the turning amplitude, and the resulting path is smoother and better meets the turning angle constraints during navigation.

Claims

1. A multi-constraint optimized UUV particle swarm three-dimensional path planning method, characterized in that: The following steps are involved: Step 1: Establish a three-dimensional seabed map in the UUV navigation environment, set the position coordinates of the starting point and the end point, and generate the offset of the moving direction and the dimension of the solution; Step 2: Initialize the particle swarm parameters, that is, initialize the population pop size , dimension dim, gene value pop max and pop min , total number of iterations, learning factors c1 and c2, inertia weight w and speed range v max and v min ; Step 3: Initialize the population of particle swarm; Step 4: Calculate the initial fitness value, first linear interpolation and coordinate reshaping: linear interpolation of the input variables according to the range of the starting point and the end point to obtain the coordinates cor, the formula is as follows: cor i,1 =cor i,1 *(end(1)-start(1))+start(1) (1) cor i,2 =cor i,2 *size(2) (2) cor i,3 =cor i,3 *size(3) (3) Among them, cor i,j is the value of the i-th point on the j-th coordinate axis; start(1) is the position of the starting point of the map on the x-axis; end(1) is the position of the end point of the map on the x-axis; size(2) is the length of the map on the y-axis; size(3) is the length of the map on the z-axis; Then normalize and map: normalize the first dimension x coordinate of the coordinate cor and map it to the range of [0,1], perform linear mapping according to the range of the starting point and the end point, restore the normalized x coordinate to the actual coordinate value, and scale the second dimension y coordinate and the third dimension z coordinate of the coordinate cor according to the scaling factor; Then perform B-spline interpolation: perform B-spline interpolation on the coordinates cor to obtain the coordinates cor of the interpolated point all , calculate the velocity v based on the coefficients of B-spline interpolation all , calculate the acceleration a based on the velocity difference and the interpolated breakpoints all , the calculation formula is as follows: cor i =[x i ,y i ,z i ],x i ∈R,y i ∈R,z i ∈R (4) color={color1,color2,...,color n } (5) cor′=pointN(x′,cor) (6) cor″′=[start,cor i ″,end] (8) cor all ″′,v all ″′,a all ″′=B-splineinterpolation(cor″′) (9) cor all,i =[x all,i ,y all,i ,z all,i ],i=1,2,...,m (10) Among them, x′ is the parameter vector; cor is the original coordinate point; cor′ is the adjusted coordinate point; cor i ″ is the coordinate point after scaling the coordinate point; cor″′ is the coordinate point after adding the start point and the end point at the beginning and the end of the coordinate point; cor all ″′ is the position information of all points obtained by interpolation; v all ″′ is the velocity information of all points obtained by interpolation; a all ″′ is the acceleration information of all points obtained by interpolation; cor all are the position coordinates of all points after interpolation; Then calculate the curvature: Calculate the curvature k based on the acceleration and velocity, the formula is as follows: Among them, k i is the curvature of the i-th point on the curve; v i is the velocity at the i-th point on the curve after interpolation; a i is the acceleration at the i-th point on the curve after interpolation; Continue with constraint processing: set the coordinates cor all The z coordinate of the site map is interpolated in two dimensions to obtain the interpolated height z interp , the penalty term of the statistical constraints, including that the insertion height cannot be lower than the site map and the coordinates cannot exceed the boundary, is as follows: Where m is the number of points on the interpolated path; z interp It is the height information of the interpolation point on the map; restraint(1) is a non-collision constraint; restraint(2) is a non-exceeding boundary constraint; restraint(3) is a To ensure that the site size cannot be exceeded; combine all constraints into the total constraint res, take constraint restraint(1), The sum of restraint(2) and restraint(3); The objective function calculation after optimizing the four objectives is as follows: calculate the total distance dis, that is, the sum of the Euclidean distances between all points, calculate the total height height, that is, the sum of the z coordinates of all points, and calculate the curvature mean k , and the penalty term res of the constraint condition, combine the above calculation results, and use the weighted coefficient to get the objective function value fitness, the formula is as follows: fitness=0.5*dis+0.2*height+0.3*mean k +res*10 3 (18) Where dis is the total length of the path; height is the total height of the UUV during navigation; mean k is the average curvature of the path; res is the sum of all constraints; fitness is the function value of multi-objective optimization; Step 5: Find the initial extreme value, fitness best is the optimal fitness value, indicating the best fitness value currently found, and initializes the best fitness best is positive infinity, index best is the optimal individual index, indicating the position of the individual corresponding to the optimal fitness value in the population, and initializing index best is 0, traverse each individual in the population, calculate its fitness value, and compare it with the current optimal fitness value fitness best Compare, if the individual's fitness value is less than the current optimal fitness value fitness best , then update fitness best is the fitness value of the individual, and the index of the individual is recorded as index best , and finally according to index best Find the corresponding z best is the group extreme value position and the corresponding is the extreme fitness value of the group; Step 6: Set the simulated annealing temperature drop process, where the starting value is 200 and the ending value is 0, indicating that the temperature will drop from the starting value to 0. The formula is as follows: Among them, temp i is the i-th temperature value; iterations is the total number of temperature values ​​that need to be generated; Since the initial temperature is 200, when i=0, temp0=200; when i=iterations-1, temp iterations-1 =0; when i takes any value, temp i will change between 200 and 0. As the number of iterations increases, the value of i keeps changing, thus achieving a continuous decrease in temperature. Step 7: Record the historical optimal solution, define and initialize the vector History to store the historical optimal solution of each iteration, and convert the optimal solution of the current iteration into fitness zbest Add to the History vector; Step 8: Update speed, the formula is as follows: V i,j =wV i,j +c1r1(pbest i,j -x i,j )+c2r2(gbest j -x i,j ) (20) Among them, V i,j is the velocity of the ith particle in the jth dimension; w is the inertia weight; c1 and c2 are acceleration constants; r1 and r2 are random numbers between [0,1]; pbest i,j is the position of the extreme value of the i-th individual in the j-th dimension; gbest j is the position of the group extreme value in the jth dimension; x i,j is the position of the i-th particle in the j-th dimension; In order to ensure that the speed range is within v min and v max The speed that exceeds the range is corrected, and the formula is as follows: Among them, v j,i is the velocity of particle j in the i-th dimension; v min is the lower bound of velocity; v max is the upper bound of the velocity; Step 9: According to the updated speed, update the particle's position. The formula is as follows: x j,i =x j,i +v j,i (22) Among them, x j,i is the position of particle j in the i-th dimension; v j,i is the velocity of particle j in the i-th dimension; Correct the new position to ensure that it is within the specified range. The formula is as follows: Among them, x j,i is the position of particle j in the i-th dimension; ub i is the upper bound on the i-th dimension; lb i is the lower bound on the i-th dimension; Step 10: Update and mutate the population; Step 11: Determine whether to perform local search optimization by using a random number. If the random number is less than 0.5, proceed to step 12; otherwise, proceed to step 13. Step 12: Conduct further local search optimization; Step 13: Record the optimal value of each iteration; Step 14: Determine whether the number of iterations of the entire population has been reached, if yes, go to step 15; otherwise, go back to step 8; Step 15: Output the optimal value.

Images