# Multistrategy Integrated Marine Predator Algorithm Applied to 3D Surface WSN Coverage Optimization

Authors: Zhendong Wang; Hang Xiao; Shuxin Yang; Junling Wang; Soroosh Mahmoodi
Journal: Wireless Communications and Mobile Computing
 (2022)

---

## Abstract

Achieving maximum network coverage with a limited number of sensor nodes is key to node deployment of wireless sensor network (WSN). This paper proposes an improved marine predator algorithm (IMPA) for 3D surface wireless sensor network deployment. A population evolution strategy based on random opposition-based learning and differential evolution operator is proposed to enrich the population diversity and improve the global search capability of the algorithm. The grouping idea of the Shuffled Frog Leaping Algorithm (SFLA) is then introduced. A local search strategy based on the SFLA is proposed to replace the FADs effect of MPA and enhance the ability of the algorithm to escape from the local optimum. A quasireflected opposition-based learning strategy is also presented to improve the optimization accuracy, accelerate the convergence speed of the algorithm, and improve the quality of the solution. Fifteen benchmark functions are selected for testing. The results are compared with seven different algorithms. The results show that the improved algorithm has excellent optimization performances. Finally, the IMPA is applied to optimize WSN coverage on 3D surfaces. The experimental results show that the proposed IMPA has good terrain adaptation and optimal deployment capabilities. It can improve the coverage of the network, reduce the deployment cost, and extend the network life cycle.


---

## Body

## 1. Introduction

Wireless sensor network (WSN) is a type of network composed of sensor nodes deployed in the monitoring area through self-organization and multihop, which can sense, process, and transmit the information of the perceived region in real-time. WSN has been widely applied to military action, space exploration, environmental science, smart home, fine agriculture, emergency rescue, and other fields [1–5]. The coverage problem is the first fundamental issue faced by WSN configuration. The sensor nodes may be arbitrarily distributed in the configuration area, which reflects the status of a certain area of the WSN being monitored and tracked. And it involves communication quality, network life, and other issues. So, it directly affects the service quality of the entire network [6]. The sensor has a very limited battery due to the size of the sensor itself. Especially in the harsh deployment area, it is not feasible to replenish the energy of the sensors by manual means. Under certain coverage conditions, it is necessary to reasonably utilize the energy of the nodes to extend the network life cycle. As a result, achieving maximum network coverage with a limited number of sensor nodes becomes a key within the sensor node deployment.At present, the coverage optimization algorithms for traditional 2D planes are very well developed. However, they cannot meet realistic and complex environments [7–12]. Compared with the 2D deployment environment, the 3D sensor node deployment is more complex and the network maintenance is more difficult. In addition, the existing 2D coverage methods cannot be directly applied in the 3D environment, otherwise a large number of coverage holes will be generated. Thus, the problem of WSN node deployment based on the 3D environments has become a vital research [13–15]. WSN node deployment for 3D environments can be divided into two categories. One category is node deployment in 3D space, such as smart cities [16], smart factories [17], and underwater wireless networks [18]. The other category is node deployment of WSNs on 3D surfaces, such as node deployment on mountains [19]. The WSN coverage of 3D surfaces as a particular coverage scenario for the 3D environments has been little studied in this aspect. However, there are still many scholars who have proposed appropriate coverage methods.Swarm intelligence algorithm is an effective method for solving many practical problems [20–22], and it has been used by many scholars to solve the 3D surface coverage optimization problem. The marine predator algorithm (MPA) as a swarm intelligence algorithm was proposed by Faramarzi et al. [23] in 2020. The MPA was proposed to simulate the predatory behavior of marine organisms. Their research was inspired by the different hunting methods of ocean life and the speed of movement between predator and prey. Although the MPA has significant advantages in optimization problems and engineering applications, it still has some disadvantages such as the low solution accuracy, the slow convergence, and the poor local exploitation capability. In order to maximize the 3D surface WSN coverage and prolong the network life cycle, aiming at the defects of the MPA, this paper proposes a multistrategy integrated improved marine predator algorithm called IMPA. In the case of fewer sensor nodes, it can effectively improve the network coverage and extend the network life cycle.This algorithm introduces a random opposition-based learning strategy and differential evolution operator to enrich the population diversity and expand the population search range. A grouping idea of the Shuffled Frog Leaping Algorithm (SFLA) is proposed to improve the exploitation ability of the algorithm. The global information and local search within the group are employed to guide the worst individual to approach the optimal individual. It can accelerate the convergence of the algorithm and reduce the risk of falling into the local optimum issues. To further speed up the convergence of the algorithm, a quasireflected opposition-based learning (QOBL) mechanism is proposed. And the boundary strategy is improved to enhance the quality of the solution. The algorithm is compared with the MSIMPA [24], the MPA, the Enhanced Sparrow Search Algorithm (ESSA) [25], the Improved Grey Wolf Optimizer (IGWO) [26], the Whale Optimization Algorithm (WOA) [27], the Particle Swarm Optimization (PSO) [28], and the Differential Evolution Algorithm (DE) [29] on 15 benchmark functions. The comparisons verify a strong optimization performance of the proposed algorithm. Finally, the proposed algorithm is applied to the WSN coverage optimization of two 3D terrains. The simulation results show that the IMPA can effectively improve the coverage of nodes, save the deployment cost of the network, and improve the network life cycle, which further verifies the practicability and effectiveness of the algorithm. In this study, the main contributions can be summarized as follows:
(1)
Most researchers calculate coverage using a 0-1 perception model, which is overly idealized in 3D environment. Therefore, this paper proposes a probabilistic coverage model suitable for 3D surfaces(2)
A population evolution strategy based on random opposition-based learning and differential evolution operator is proposed to enrich the population diversity while improving the global search ability of the algorithm(3)
A local search strategy based on SFLA is proposed to replace the local search mechanism of the MPA to enhance the algorithm’s local exploitation capability and the ability to jump out of the local optimum(4)
A novel 3D surface WSN coverage optimization technique based on multistrategy integrated improved marine predator algorithm is presented. Simulation experiments are conducted in two 3D terrains and the method has higher coverage and longer life cycle compared to other four algorithmsThe rest of the paper is organized as follows. Section2 presents the current research status and related work on 3D surface coverage and MPA algorithms. Section 3 introduces the relevant models for WSN node deployment on 3D surfaces. Section 4 describes the proposed IMPA algorithm in detail. Section 5 compares and analyzes the algorithms on the benchmark test functions. Section 6 presents the application of IMPA for WSN coverage in two 3D terrains. Finally, Section 7 concludes and proposes future research.


## 2. Related Work

At present, the research on WSN coverage control in 3D environment is a hot topic. For the WSN coverage optimization problem on 3D surfaces, some scholars use traditional techniques to improve coverage and ensure network service quality. For example, Anand et al. proposed a node coverage method applicable to 3D surfaces based on the terrain segmentation method of the Voronoi diagram [30]. Their simulation results demonstrated improvement in the algorithm’s efficacy. A 3D deployment method based on the virtual force algorithm was proposed [31]. The simulation results show that the algorithm can be adapted to all types of surfaces and achieve full coverage while maintaining network connectivity. Kim [32] introduced a mobile sensor network for the coverage problem in mountainous terrain and proposed a 3D coverage method by considering the physical characteristics such as slopes in mountainous terrain. Simulation experimental results show that the method is significantly improved in terms of network coverage and coverage time. Full coverage of the monitoring area can be achieved by deploying an appropriate number of sensors. However, the method is not applicable in large-scale node deployment scenarios. These traditional techniques have provided important solutions and valuable research results for solving the coverage optimization problems of WSN, but there are some shortcomings. For example, the structure of some algorithms is too complex, which makes it difficult to meet the service requirements of users in practical deployment.Since the swarm intelligence optimization algorithm has shown excellent results in various optimization problems and engineering applications, it opens up a new path for the research of coverage optimization of WSN, many researchers have utilized the swarm intelligence algorithm to optimize the WSN coverage of 3D surfaces [33–35]. For example, Li et al. [34] proposed a Modified Parallel Tunicate Swarm Algorithm (MPTSA) based on the improved traditional parallel strategy and applied it to the 3D surface coverage optimization. Their simulations and experimental results show the MPTSA can effectively improve the coverage. However, this method does not consider the lifetime of the network. Literature [36] proposed an innovative method for determining the 3D surface perception blind area for the 3D spatial target coverage problem. A differential evolution algorithm was employed to solve the 3D surface nodes deployment problem to realize full coverage of the monitoring area with the least number of nodes. However, the proposed method for judging the 3D perception blind area ignored the limitations where both sensing nodes and monitoring points are located on the terrain surface. Wang and Xie [37] proposed a more comprehensive method for determining the 3D surface perception blind area based on the mentioned limitations of perceptual blind area determination. They also innovatively enhanced the grid method commonly used in 2D WSN coverage to make it applicable to the coverage calculation of 3D surfaces. Moreover, an Enhanced Grey Wolf Optimizer (EGWO) using a hybrid inner and outer layer location update strategy was proposed and employed to the deployment optimization of 3D surfaces. The experimental results show that the EGWO has excellent deployment ability. However, the coverage perception model of this work is a binary perception model, which is too ideal in complex 3D terrain.Currently, the research on MPA focuses on improving and applying the algorithm itself. Houssein et al. [38] proposed a MPAOBL-GWO algorithm incorporating the opposition-based learning and the Grey Wolf Optimizer (GWO). The GWO was applied as the local search method of MPA. It improved the exploitation capability of the algorithm and avoided falling into the local optimum issues. However, the convergence speed of the algorithm still needs to be improved. Abdel-Basset et al. [39] proposed an improved marine predator algorithm (IMPA) by using the improved population strategy based on the population’s average fitness. This strategy enhanced the performance and convergence speed of the algorithm. However, the risk of falling into the local optimum searching for individuals while reducing the population diversity was increased. Literature [40] introduced a differential evolution operator into the MPA. It proposed an enhanced marine predator algorithm (EMPA) for identifying static and dynamic photovoltaic model parameters. The EMPA maintained the diversity of the population in the search process and increased the probability of the algorithm away from the local optimum. However, the enhancement of the optimization accuracy of the improved algorithm was very limited. In literature [24], a multistrategy improved marine predator algorithm (MSIMPA) was proposed. It employed the chaotic mapping and opposition-based learning strategy to improve the population’s generation quality within the initial stage. Moreover, a grouping learning strategy to crossover the population operations was proposed to improve the algorithm’s optimization accuracy and convergence speed. The simulation results show that the algorithm could effectively improve the deployment coverage of the 2D WSN.Based on the above related works, some achievements have been achieved in the study of WSN coverage optimization for 3D surfaces and the improvement of the MPA algorithm. However, there are some shortcomings in these studies. For example, the coverage perception model is a binary model, which is too idealized in practical applications. In addition, since the 3D surface is more complex compared to the 2D plane, the coverage of some algorithms still needs to be further improved. Therefore, a lot of research is still needed for coverage optimization of 3D surface WSN to achieve higher coverage while obtaining a longer network life cycle using fewer sensor nodes. To address the above issues, this paper proposes a coverage optimization method based on multistrategy integrated improved marine predator algorithm.


## 3. 3D Surface WSN Coverage Model Description

### 3.1. Perceptual Model of Sensor Nodes

In WSN, suppose that the nodes of a wireless sensor are a setS=s1,s2,⋯,sn where the 3D coordinate of sj is denoted by xj,yj,zj that j=1,2,⋯,n. The set of the monitored nodes within the monitoring area is M=m1,m2,⋯,mn, where the 3D coordinate of mk is denoted by xk,yk,zk that k=1,2,⋯,n. The sensing and communication radiuses of each node are Rs and Rc, respectively. The Euclidean distance between node sj and node mk can be defined by
(1)dissj,mk=xj−xk2+yj−yk2+zj−zk2.The common 0-1 perception model is mainly used in the WSN coverage of the 2D environment. In the 3D surface environment, it is necessary to establish a coverage model that is closer to the real world. Therefore, this paper adopts the probabilistic perception model. The probability of the monitoring nodemk can be perceived by the node sj in the following:
(2)Pcovsj,mk=11+α∗dissj,mkβ,dissj,mk≤Rs,0,dissj,mk>Rs,where α and β are type parameters related to the physical characteristics of the sensor. Usually, the value of β is from 1 to 4 and α is an adjustable parameter.Therefore, the joint perception probability of all nodes in the setS to monitor node mk can be evaluated by
(3)PCsall,mk=1−∏i=1n1−Pcovsj,mk,where sall denotes all the nodes of the sensor within the monitoring range.


### 3.2. Judgment of 3D Perception Blind Area

There are often many obstacles in the deployment of sensor nodes in the real world. These obstacles will block the sensing signals of the sensor nodes, causing a substantial reduction in the effective monitoring area of the nodes.As shown in Figure1, even though nodes A, B, and C are all located in the perception range of node S, the monitoring points A and B are not sensed by the node S because the node signals are blocked by the obstacle (the shaded area in Figure 1). The lack of sensing causes the formation of a particular area for the 3D surface node deployment, which is called the 3D perception blind area [41].Figure 1
Schematic diagram of the 3D perception blind area.For instance, the general method of the 3D perception blind area which is monitored from pointB can be determined by the following procedure. Assuming that the coordinate of node S is xs,ys,zs, the coordinate of point B is xb,yb,zb, the spatial surface equation is z=fx,y, and point B is in the sensing range of node S. The line segment SB is obtained by connecting both points S and B, and then the equation of the straight spatial line SB can be created. As shown in Equation (4), the simultaneous equations of the spatial line SB and the spatial surface z. It can obtain the number η of solutions to Equation (4) under conditions x∈xs,xb and y∈ys,yb.
(4)x−xsxb−xs=y−ysyb−ys=z−zszb−zsz=fx,y.The determination of the 3D perception blind area between nodeS and point B depends on the value of η, which can be divided into the following two situations.Ifη=0, there is no perception blind area between node S and point B. It defines that the surface does not intersect with the line segment SB, which means the point B is covered by the node S.Ifη≥1, there is a perception blind area between node S and point B, which means that the point B is not covered by the node S. Figure 1 intuitively shows that the surface intersects with the line SB, so the node S cannot monitor the point B.The procedure mentioned above is an effective method to determine the existence of blind areas. However, there are some limitations. For instance, Figure1 shows that the line segment AS has no intersection with the surface, but there is an obstacle between the point A and the node S. Therefore, the node S cannot cover the point A. For this situation, a determination rule is added in the following.Suppose the midpoint of the line segmentAS is Q1xq,yq,zq. If the value for the z-axis of the point Q1 is smaller than the value of the function corresponding to the surface equation z=fx,y, then there would be a perceptual blind area between both points, which means node S cannot monitor the point A. The condition of the method is presented in the following:
(5)PcovS,A=0ifzq≤fxq,yq,1otherwise.It summarizes that a monitoring point should be perceived by the nodeS when the following three conditions are satisfied. First, the evaluated value of Equation (2) not equal to 0; second, Equation (4) has no solution, and finally, the evaluated value of Equation (5) becomes 1.


### 3.3. Coverage Description

It is not easy to divide the monitoring region into small regions of equal areas when solving the problem of 3D WSN coverage. Besides, it is not also necessary to do this division. A surface division based on the grid method was proposed [37]. The method gave a reasonable error condition for the number of grid divisions. The basic idea is to project the surface vertically into a two-dimensional plane. Then, the projection region is divided into several small grids with equal areas. The area of the small surface corresponding to each small grid can be calculated by the calculus idea and usage of the surface with the integral surface formula. The summation of the areas of the small surfaces corresponding to all small grids would be the total area of the whole surface. The calculation method of the surface area is shown in
(6)z=fx,y,SF=∬Dxy1+dzdx2+dzdy2dxdy,where SF denotes the area of the surface and Dxy is the integration region consisting of x and y.In this paper, the grid method is employed to calculate the coverage of the 3D surface. First, the projection area of the surfaceF mapped to the 2D plane is divided into l∗w grids, where l represents the number of grids in the x-axis direction, and w represents the number of grids in the y-axis direction. That is to say, the target area is divided into several cube grids of equal size. Consider the set SF=s1F,s2F,⋯,snF represents the area of the small surface, where n=l∗w. The projection region is divided into several small grids of equal area and locates the monitoring point at the center of the small surface corresponding to the small grid. If the monitoring point can be perceived by the node S, it means that the node S could cover the small surface where the monitoring point is located. The final coverage area can be obtained by accumulating the areas of all the small surfaces which can be monitored. Therefore, the coverage CF of the 3D surface can be described as follows:
(7)CF=1SF∑x=1l∑y=1wPCsall,mw⋅x−1+y×sw⋅x−1+yF,where sall denotes all the sensors deployed in the target area, the subscript wx−1+y is used to calculate the index of the grid, that is, the index of the monitoring point. Equation (7) is the proposed objective function which is optimized by the algorithm presented in this research.


## 3.1. Perceptual Model of Sensor Nodes

In WSN, suppose that the nodes of a wireless sensor are a setS=s1,s2,⋯,sn where the 3D coordinate of sj is denoted by xj,yj,zj that j=1,2,⋯,n. The set of the monitored nodes within the monitoring area is M=m1,m2,⋯,mn, where the 3D coordinate of mk is denoted by xk,yk,zk that k=1,2,⋯,n. The sensing and communication radiuses of each node are Rs and Rc, respectively. The Euclidean distance between node sj and node mk can be defined by
(1)dissj,mk=xj−xk2+yj−yk2+zj−zk2.The common 0-1 perception model is mainly used in the WSN coverage of the 2D environment. In the 3D surface environment, it is necessary to establish a coverage model that is closer to the real world. Therefore, this paper adopts the probabilistic perception model. The probability of the monitoring nodemk can be perceived by the node sj in the following:
(2)Pcovsj,mk=11+α∗dissj,mkβ,dissj,mk≤Rs,0,dissj,mk>Rs,where α and β are type parameters related to the physical characteristics of the sensor. Usually, the value of β is from 1 to 4 and α is an adjustable parameter.Therefore, the joint perception probability of all nodes in the setS to monitor node mk can be evaluated by
(3)PCsall,mk=1−∏i=1n1−Pcovsj,mk,where sall denotes all the nodes of the sensor within the monitoring range.


## 3.2. Judgment of 3D Perception Blind Area

There are often many obstacles in the deployment of sensor nodes in the real world. These obstacles will block the sensing signals of the sensor nodes, causing a substantial reduction in the effective monitoring area of the nodes.As shown in Figure1, even though nodes A, B, and C are all located in the perception range of node S, the monitoring points A and B are not sensed by the node S because the node signals are blocked by the obstacle (the shaded area in Figure 1). The lack of sensing causes the formation of a particular area for the 3D surface node deployment, which is called the 3D perception blind area [41].Figure 1
Schematic diagram of the 3D perception blind area.For instance, the general method of the 3D perception blind area which is monitored from pointB can be determined by the following procedure. Assuming that the coordinate of node S is xs,ys,zs, the coordinate of point B is xb,yb,zb, the spatial surface equation is z=fx,y, and point B is in the sensing range of node S. The line segment SB is obtained by connecting both points S and B, and then the equation of the straight spatial line SB can be created. As shown in Equation (4), the simultaneous equations of the spatial line SB and the spatial surface z. It can obtain the number η of solutions to Equation (4) under conditions x∈xs,xb and y∈ys,yb.
(4)x−xsxb−xs=y−ysyb−ys=z−zszb−zsz=fx,y.The determination of the 3D perception blind area between nodeS and point B depends on the value of η, which can be divided into the following two situations.Ifη=0, there is no perception blind area between node S and point B. It defines that the surface does not intersect with the line segment SB, which means the point B is covered by the node S.Ifη≥1, there is a perception blind area between node S and point B, which means that the point B is not covered by the node S. Figure 1 intuitively shows that the surface intersects with the line SB, so the node S cannot monitor the point B.The procedure mentioned above is an effective method to determine the existence of blind areas. However, there are some limitations. For instance, Figure1 shows that the line segment AS has no intersection with the surface, but there is an obstacle between the point A and the node S. Therefore, the node S cannot cover the point A. For this situation, a determination rule is added in the following.Suppose the midpoint of the line segmentAS is Q1xq,yq,zq. If the value for the z-axis of the point Q1 is smaller than the value of the function corresponding to the surface equation z=fx,y, then there would be a perceptual blind area between both points, which means node S cannot monitor the point A. The condition of the method is presented in the following:
(5)PcovS,A=0ifzq≤fxq,yq,1otherwise.It summarizes that a monitoring point should be perceived by the nodeS when the following three conditions are satisfied. First, the evaluated value of Equation (2) not equal to 0; second, Equation (4) has no solution, and finally, the evaluated value of Equation (5) becomes 1.


## 3.3. Coverage Description

It is not easy to divide the monitoring region into small regions of equal areas when solving the problem of 3D WSN coverage. Besides, it is not also necessary to do this division. A surface division based on the grid method was proposed [37]. The method gave a reasonable error condition for the number of grid divisions. The basic idea is to project the surface vertically into a two-dimensional plane. Then, the projection region is divided into several small grids with equal areas. The area of the small surface corresponding to each small grid can be calculated by the calculus idea and usage of the surface with the integral surface formula. The summation of the areas of the small surfaces corresponding to all small grids would be the total area of the whole surface. The calculation method of the surface area is shown in
(6)z=fx,y,SF=∬Dxy1+dzdx2+dzdy2dxdy,where SF denotes the area of the surface and Dxy is the integration region consisting of x and y.In this paper, the grid method is employed to calculate the coverage of the 3D surface. First, the projection area of the surfaceF mapped to the 2D plane is divided into l∗w grids, where l represents the number of grids in the x-axis direction, and w represents the number of grids in the y-axis direction. That is to say, the target area is divided into several cube grids of equal size. Consider the set SF=s1F,s2F,⋯,snF represents the area of the small surface, where n=l∗w. The projection region is divided into several small grids of equal area and locates the monitoring point at the center of the small surface corresponding to the small grid. If the monitoring point can be perceived by the node S, it means that the node S could cover the small surface where the monitoring point is located. The final coverage area can be obtained by accumulating the areas of all the small surfaces which can be monitored. Therefore, the coverage CF of the 3D surface can be described as follows:
(7)CF=1SF∑x=1l∑y=1wPCsall,mw⋅x−1+y×sw⋅x−1+yF,where sall denotes all the sensors deployed in the target area, the subscript wx−1+y is used to calculate the index of the grid, that is, the index of the monitoring point. Equation (7) is the proposed objective function which is optimized by the algorithm presented in this research.


## 4. Improved Marine Predator Algorithm

### 4.1. Marine Predator Algorithm

The marine predator algorithm [23] is a novel metaheuristic algorithm that simulates the foraging strategy of marine organisms. It updates the position of predator and prey by switching between the Lévy flight and the Brownian motion. Both predator and prey are solution candidates, since each prey in the ecosystem is a potential predator. Like most metaheuristic algorithms, MPA generates the initial population randomly in the search space. The mathematical model is described as follows:
(8)X0=randXmax−Xmin+Xmin,where rand is a random number that obeys a uniform distribution of 0,1 and Xmin and Xmax denote the lower and upper limits of the search space, respectively.In MPA, the location of the prey can be represented by a matrix. The corresponding mathematical expression is described as follows:(9)Prey=X1,1⋯X1,d⋮⋱⋮Xn,1⋯Xn,dn×d,where n denotes the number of populations and d is the dimension of the search space. Xi,j denotes the current position of the i-th prey in the j-th dimension. The individual with the best fitness value is called the top predator. It is used to construct the Elite matrix, whose size is the same as the Prey matrix. The mathematical expression of the Elite matrix is presented in the following:
(10)Elite=X1,1I⋯X1,dI⋮⋱⋮Xn,1I⋯Xn,dIn×d,where XI represents the top predator vector. The Elite matrix is obtained by replicating this vector n times. The searching procedure of the MPA algorithm can be separated into three stages in accordance with the speed variation of marine biological predation, which is described as follows.


#### 4.1.1. High-Speed Ratio Phase

This phase occurs in the first third of the number of iterations, when the population is in the exploration phase. At this time, the speed of the predator is greater than that of the prey, and the mathematical model of prey position update is described as follows:(11)Whilet≤13tmax,i=1,2,⋯,n,(12)stepsizei=RB⊗Elitei−RB⊗Preyi,(13)Preyi=Preyi+P×R⊗stepsizei,where t and tmax denote the current and the maximum number of iterations, respectively. RB is a vector of random numbers obeying the normal distribution used to represent the Brownian motion. The symbol ⨂ denotes the entry-wise multiplication, R is a random number between 0 and 1, and P=0.5 is a constant.


#### 4.1.2. Unit-Speed Ratio Phase

This phase occurs between one-third and two-thirds of the number of iterations and gradually indicates the population’s transition from exploration to exploitation. At this phase, the predator and the prey movement speed are equal. In this stage, the first half of the population is used for the exploration, and the movement strategy of the prey is the Lévy motion. The remaining half of the population would then be used for the exploitation where the movement strategy of the predator would be the Brownian motion. The mathematical model of this phase is described as follows:(14)While13tmax<t<23tmax,i=1,2,⋯,n/2,(15)stepsizei=RL⊗Elitei−RL⊗Preyi,(16)Preyi=Preyi+P×R⊗stepsizei,(17)While13tmax<t<23tmax,i=n/2,⋯,n,(18)stepsizei=RB⊗RB⊗Elitei−Preyi,(19)CF=1−ttmax2t/tmax,(20)Preyi=Elitei+P×CF⊗stepsizei,where RL is a vector of random numbers based on the Lévy distribution and CF is an adaptive parameter that controls the movement step of the predator.


#### 4.1.3. Low-Speed Ratio Phase

A late iteration phase occurs in the last remaining third of the whole iteration. At this phase, the predator moves slower than the prey within a Lévy’s motion, which adopts the population into a local exploitation strategy. The mathematical model of this phase is described as follows:(21)Whilet≥23tmax,i=1,2,⋯,n,(22)stepsizei=RL⊗RL⊗Elitei−Preyi,(23)Preyi=Elitei+P×CF⊗stepsizei.


#### 4.1.4. Eddy Formation and FADs Effect

The formation of eddies or Fish Aggregation Devices (FADs) affects the foraging behavior of the marine predators. This strategy mainly improves searching for individuals to escape the local optimal solutions. Its mathematical model is described as follows:(24)Preyi=Preyi+CFXmin+RL⊗Xmax−Xmin⊗Ur≤Pf,Preyi+Pf1−r+rPreyr1−Preyr2r>Pf,where U is a binary vector containing 0 and 1 and Pf=0.2 represents the probability of influence on the optimization process. r is a random number in the range 0,1, and the subscripts r1 and r2 represent the random index values of the prey matrix. In addition, the MPA algorithm has memory storage to remember the history of the best-reached position, which is similar to an essential greedy strategy. During the algorithm iteration, the fitness value of each individual in the current iteration is compared with the optimal solution of the previous iteration. The previous solution will then be replaced if the current’s fitness is better.


### 4.2. Population Evolution Strategy Based on the Random Opposition-Based Learning and Differential Evolution Operator

The optimization process in each phase of MPA is highly depends on the prey matrix (described in section4.1). The MPA algorithm randomly initializes the population and then relies only on the Brownian or Lévy motion strategies to update the individual positions. The late iteration of the algorithm will reduce population diversity, narrow search range, reduce exploration ability, and tends to stagnate in the local optimum. Usually, in the optimization process of swarm intelligence algorithms, when the searched individual is close to the optimal solution, the population will quickly approach the direction of the optimal solution and accelerate the convergence of the algorithm. On the contrary, the convergence rate becomes slower.To solve this problem, opposition-based learning [42] (OBL) provides a search strategy that considers both the current solution and its corresponding inverse solution, expanding the search range and avoiding stagnation of the current solution. The research shows that the probability of the opposite solution of the population close to the global optimal solution is greater than that of the random solution. Therefore, many researchers have attempted to use OBL to improve both the quality of the population and the performance of the algorithm [43–48].To enhance the diversity of the population further, Long et al. [49] proposed a random opposition-based learning (ROL) strategy based on OBL, which could effectively improve the ability of the algorithm to jump out of the local optimum. Suppose a feasible solution of a d-dimensional space is X=x1,x2,⋯,xd, where x1,x2,⋯,xd∈R and xj∈lbj,ubj. The random opposite solution of X would be X^r=x^r,1,x^r,2,⋯,x^r,d and its calculation formula can be presented as follows:
(25)x^r,j=lbj+ubj−R1×xjj=1,2,⋯,d,where R1∈0,1 is a random number.Given the effectiveness of the ROL strategy for the improvement of algorithm performance and inspired by the mutation, crossover, and selection operators in the differential evolution algorithm (DE), this paper proposes a population evolution strategy based on ROL and differential evolution operator (RDES) to generate a higher quality population.First, the ROL population of the initial population is evaluated by Equation (25). Because the positions of all individuals in the ROL population may not be better than those of the initial individuals, both the initial and ROL populations are merged to expand the population searching range. Then, the differential evolution operator is considered based on the merged population to enrich the population’s diversity further. The specific operations are as follows:
(1)
Mutation Operation. Since the mutation operator in the DE has a significant influence on the algorithm, it is essential to choose an effective individual mutation method. This paper used the DE/current-to-best/1 mutation operator [29] to perform mutation operation on the merged population, as shown in the following formula:
(26)Vit=Xit+F⋅Xbestt−Xkt+Xmt−Xnt,where F is scaling factor, Xit denotes the i-th individual of the current generation, Xbestt denotes the optimal individual of the current population. Xkt, Xmt, and Xnt denote three different individuals randomly selected except the current individual, where the k,m,n∈1,2,⋯,2∗N are not equal to each other, and N represents the population size. Equation (26) reveals that the mutation operator uses the individuals’ local information and the population’s global optimal information to lead the population to find an optimal solution. It expands the diversity of the population and accelerates the searching speed of the algorithm simultaneously(2)
Crossover Operation. The original individual Xit and the new generated individual Vit are cross-operated after the differential mutation according to
(27)Ui,jt=Vi,jt,rand≤Pc,Xi,jt,otherwise,where Ui,jt represents the value of the j-th gene of the current i-th individual, rand is a random number uniformly distributed in 0,1, and Pc represents the crossover probability. The algorithm would keep the gene value of Vi,jt if rand<Pc, otherwise, it replaces the corresponding gene value with Xi,jt(3)
Selection of Operation. If the crossover operation results in a new individual with a better fitness value than the original individual, the original individual would be substituted. This process can be described as follows:
(28)Xit+1=Uit,fUit≤fXit,Xit,otherwise.Finally, the fitness values of all individuals in the merged population are calculated and sorted. The topN individuals with the best fitness value are selected to construct the prey matrix using the greedy strategy.Based on the above descriptions, the mainframe of the RDES is given in Algorithm1.Algorithm 1: Pseudocode of the RDES.
1:    Initialize the population withN individuals2:Fori≤Ndo3:     Using Eq. (25) to generate the ROL population4:End for5:    Merge initialization population and ROL population6:Fori≤2Ndo7:     Execute the mutation operation by Eq. (26)8:Ifrand<Pcthen9:      Execute the crossover operation by Eq. (27)10:End if11:    Using Eq. (28) to execute the selection operation12:End for13:    Select the topN fitness individuals from the mergedPopulation as the initial Prey matrix


### 4.3. Local Search Strategy Based on Improved SFLA

SFLA is a swarm intelligence algorithm proposed by Eusuff et al. [50]. They simulated the population’s mechanism of information sharing and communication during the frog foraging. The algorithm combines the advantages of the Memetic Algorithm (MA) [51] and the PSO [28]. Moreover, it benefits few control parameters, is an easy implementation, and has a strong global searching ability. Therefore, some researchers have hybridized SFLA with other algorithms to improve the optimization performance of the algorithm. For instance, in [52], an improved ant colony algorithm was combined with SFLA for solving the path planning problem, which accelerated the convergence of the algorithm. In [53], in order to improve the local search capability of the PSO, SFLA was integrated with PSO and used for the detection of sonar images.The MPA uses the FADs effect (Equation (24)) to support individuals’ searching for escape from the local optimum trap. This method lacks the guidance of population information to update the individuals’ position and limited the local development capability. Finding a more optimal solution becomes complicated once the population falls into the local optimum. The SFLA divides the whole population into several groups based on its unique grouping operator. The evolution of the individual is guided by the local optimal solution in the group and the global optimal solution of the population. It realizes the information sharing and communication between individuals. Therefore, the algorithm would not easily fall into the local optimum and would be more conducive to find the global optimal solution [22]. In this paper, by considering the grouping idea of SFLA, both MPA and SFLA algorithms are combined, and a local search algorithm based on the improved SFLA is proposed. The proposed local search algorithm replaces the FADs effect to improve the local development ability of MPA.The main idea of the SFLA is described in the following. First, all the population individuals are arranged in descending order according to their fitness values and divided intom groups. The first individual is assigned to the first group, the second to the second group, until the m-th to the m-th group. Them+1-th individual is then assigned to the first group, the m+2-th to the second group, and so on, until assigning the last one. Then, the positions of the worst individual Pw are updated according to Equations (29)–(31). The specific update method is as follows:
(29)Pnew=Pw+s×rand×Pb−Pw,(30)Pnew=Pw+s×rand×Pg−Pw,(31)Pnew=lb+rand×ub−lb,where Pnew denotes the updated position of the worst individual, Pb and Pw denotes the best and worst individuals in the group, randis a random number between 0,1, Pg represents the position of the globally optimal individual for the entire population, s is a constant representing the step size factor and that its value can be adjusted according to the problems, and ub and lb are the upper and lower bounds of the search space, respectively.Since the third location update method (as shown in Equation (31)) of the basic SFLA algorithm is relatively random, it lacks the guidance of population information, and consequently, the searching speed is slow. Therefore, this paper considers the usage of the local information of the group and the global information of the population to guide the algorithm for strengthening the ability of further search. A new location update method is proposed, as shown in
(32)Pnew=Pw+rand×Pg+Pb2−Pw.The mainframe of the improved local search strategy is shown in Algorithm2. First, the prey population is divided into m groups according to the grouping strategy which was introduced above. After the completion of grouping, the worst position of an individual would be updated using Equation (29). The original position is then replaced if the new position includes a better fitness value. Otherwise, the new position would be evaluated by Equation (30). If neither of the above two methods can generate a better position, the position of the worst individual in the group would be updated by Equation (32).Figure2 depicts the three updating methods for the position of the worst individual in each group. As it can be seen in Figure 2, Equation (29) means the best individual in the group guides the worst individual to move towards it. Equation (30) indicates that the worst individual approaches the global optimal individual, and Equation (32) represents that the population approaches the middle position between the global optimal and the intragroup optimal.Figure 2
Schematic diagram of the position update of the worst individual (the black dashed line represents the guidance direction of the population information, and the solid blue line represents the trajectory of the worst individual position update).The improved local search method transmits the information according to the population grouping. It effectively combines the global information of the population and the local search. It also improves the exploitation ability of the population. In addition, the position of the worst individuals in the group is only updated each time after grouping, which reduces the computation time and increases the convergence rate of the algorithm. The fitness value of the population is inevitably improved accordingly when the individuals of all groups are remixed into the population. Therefore, it leads the population towards the optimal solution.Algorithm 2: Pseudocode of the local search strategy.
1:  Initialization parameters: The number of groupsm; the number of iterations within the group tm;2:Fori=1:mdo3:Forj=1:tmdo4:     RecordPw,Pb, Pg5:     Update the worst Prey using Eq. (29)6:IffitnessPnew<fitnessPwthen7:Pw=Pnew8:Else9:        Update the worst Prey using Eq. (30)10:fitnessPnew<fitnessPwthen11:Pw=Pnew12:Else13:         Update the worst Prey using Eq. (32)14:End if15:End if16:End for17:End for


### 4.4. Dynamic Learning Mechanism

Ergezer et al. proposed a variant of OBL called quasireflected opposition-based learning (QOBL) [54]. They proved that the probability of the QOBL solution is closer to the optimal solution than the general OBL solution is 11/16. This paper proposes a QOBL strategy for the prey population to improve the optimization accuracy and accelerate the convergence speed of the algorithm further. The QOBL position (Preyiqr) of the prey in the generation can be described as follows:
(33)Preyiqrt=randPreyi,ci,(34)ci=minPreyi+maxPreyi2,where Preyi represents the position of the i-th prey and ci represents the midpoint of the dynamic boundary. Unlike the fixed boundary, the dynamic search boundary adjusts the maximum and minimum values in all dimensions of the current i-th individual as the upper and lower bounds of the search space. This method can further reduce the possibility of stagnating the population in the local optima, enhance the exploitation capability of the algorithm, and increase the convergence speed [55].Finally, the elite strategy is employed to retain the individuals who carry better fitness value. If the fitness value of the QOBL position’s fitness value is superior to the current individual, the new position would be assigned to the current individual. The strategy is modeled as follows:(35)Preyit+1=PreyiqrtiffPreyiqrt<fPreyit,Preyitotherwise.Figure3 shows the process of applying the dynamic learning mechanism based on QOBL to the MPA algorithm. First, the QOBL population of the prey population is evaluated. Then, the fitness values of both initial prey and the QOBL populations would be calculated separately. Finally, the individuals with better fitness values are retained for constructing new prey populations according to the Equation (35).Figure 3
Schematic diagram of population selection based on QOBL.


### 4.5. IMPA Algorithm

In this paper, a multistrategy integrated improved marine predator algorithm (IMPA) is proposed to increase both accuracy of solution and convergence speed and prevent falling into the local optimum of MPA. The pseudocode of IMPA is given in Algorithm3, and the main framework of the algorithm is demonstrated in Figure 4.Algorithm 3: The proposed IMPA.
Input: Number of population N, the maximum number of iterations tmax, the crossover probability Pc, the number of groups m, the number of iterations within the group tm, the dimensions of each individual d, and search range lb,ubOutput: Optimal fitness1:  Initialize the population of Preys withN individuals2:Whilet≤tmaxdo3:     Execute RDES to generate a high-quality population (Algorithm1)4:     Calculate the fitness and build the Elite matrix5:Ift≤tmax/3then6:      Update Prey using Eq. (13)7:Else iftmax/3<t<2∗tmax/38:        For the first half of the population (i=1,2,⋯,N/2)9:        Update Prey using Eq. (16)10:        For the remained half of the population (i=N/2,⋯,N)11:        Update Prey using Eq. (20)12:Else ift≥2∗tmax/313:        Update Prey using Eq. (23)14:End if15:     Complete memory saving and Elite update16:     Execute the local search strategy (Algorithm2)17:Fori=1:Ndo18:       Calculate the QOBL position according to Eq. (33)19:       Select the better individual according to Eq. (35)20:End for21:  Complete memory saving and Elite updating22:  Next iterationt=t+125:End whileFigure 4
IMPA flow chart.


### 4.6. Algorithm Validity Test

In this section, the Schwefel2.26 Function is selected to test the validity of the three improved strategies. The population distribution maps of IMPA and MPA algorithms are drawn. The initial and final distributions of both algorithms are shown in Figures5 and 6 while the population size and the maximum number of iterations are 50 and 20, respectively.Figure 5
Population distribution of MPA.
(a)
Initial population distribution of MPA(b)
Final population distribution of MPAFigure 6
Population distribution of IMPA.
(a)
Initial population distribution of IMPA(b)
Final population distribution of IMPAFigures5(b) and 6(b) show that the population distribution of MPA is relatively scattered after 20 iterations and the convergence speed is low. Moreover, most individuals gather near the local optimum, which means that the MPA would not easily jump out of the local optimum when dealing with multidimensional complex functions. In contrast, all individuals in the IMPA converge to the global optimum quickly, proving the effectiveness of the proposed strategies.


### 4.7. Time Complexity Analysis

Time complexity is an essential indicator for predicting the performance of an algorithm. Suppose thatN is the size of populations, D is the dimension of individuals, and T is the maximum number of iterations. The time complexity of the MPA is evaluated ON∙D∙T. The proposed IMPA algorithm would not change the structure of individual position updates in the three phases of the basic MPA algorithm. The analysis of the time complexity is explained in the following.The increased complexity of the ROL population while performing the RDES (Algorithm1) is ON, where this evaluation becomes O2∙N∙D while performing the DE evolution operation. While executing the local search strategy (Algorithm 2), the evaluated complexity is related to the number of the groups m and the number of iterations tm, where m≤N/2. Therefore, in the worst case, the increased complexity of this part would be ON. The increased complexity while computing the QOBL population would be ON. Finally, the evaluated time complexity of IMPA would be OTN+2∙N∙D+N+N, which results in ON∙D∙T. It concluded that the improved strategies proposed in this study would not increase the time complexity of the proposed IMPA.


### 4.8. Global Convergence Analysis of IMPA

The convergence analysis of the three phases within the optimization process for the basic MPA algorithm has been extensively discussed [23]. It approved that the MPA is a convergent algorithm. The conclusion of the mentioned study is reviewed in the following.Assuming thatPreyt is a random sampling sequence of MPA, the MPA eventually converges to an Elite matrix (globally optimal position) consisting of the top predators, which is limt⟶∞Preyt=Elite.The global convergence of the hybrid optimization algorithm was analyzed, and several sufficient conditions for the global convergence of the hybrid optimization algorithm were approved [56]. The approved relevant theorems of the mentioned research are represented below.Lemma 1.
fn is a monotonically nonincreasing sequence with a lower limited bound in which limt⟶∞fn=inffn [56].Lemma 2.
If an algorithm satisfies both following conditions (properties), the optimal sequence will converge globally according to a probability equal to 1 [56].Condition 3.
The greedy strategy is used to retain the elite while comparing with the historical optimal solution during the iteration. The mentioned strategy is represented in the following:(36)ft+1best=fX,fX<ftbest,ftbest,otherwise,where ftbest is the optimal fitness value found after t times of comparisons which retained to the moment t, and X is the solution used for t+1-th comparison withftbest.Condition 4.
The probability of transitioning from any nonglobal optimal point (might be set toX∗) to the corresponding level set LX∗=XfX<fX∗,X∈S is not 0, where S represents the entire search space.Theorem 5.
The proposed IMPA algorithm converges globally according to the probability equal to 1.
Assume that the objective functionf of the IMPA algorithm optimization is set to find the minimum value, fn is the random sampling sequence, ftbest is the optimal sequence of sampling points, and fgbest is the global optimal solution. As shown below, Theorem 5 can be proved by Lemma 1 and Lemma 2.Proof.
Because the three proposed improved strategies and the marine memory storage procedure in the IMPA algorithm all use the greedy strategy to retain the optimal solution of the population, the IMPA algorithm would satisfy Condition3.According to Lemma1, the limit of the optimal sequence ftbest exists. Let limt⟶∞ftbest=fnbest, and fnbest≠fgbest. Because fgbest is the global optimal solution, therefore fnbest>fgbest. Since limt⟶∞ftbest=infftbest=fnbest, it follows that for ∀t∈Z+, ftbest≥fnbest. Let X∗ be the sampling point corresponding to fnbest, so, fnbest=fX∗. X∗ is obviously a nonglobal optimal solution. Because the IMPA uses the same update operation as well as the MPA algorithm in the high-, unit-, and low-speed ratio phases, the convergence of the population sequence of IMPA would be consistent with that of the MPA and will eventually converge to the Elite matrix, limn⟶∞fn=Elite=fgbest. That is, the nonglobal optimal point X∗ can converge to the Elite matrix (global optimal point). Therefore, when the X assigns value to a sufficiently small neighborhood around the elite matrix, PLX∗=PXfX<fX∗,X∈S≠0, the Condition 4 is satisfied.In summary, the IMPA algorithm satisfies both conditions of Lemma2. Therefore, the IMPA algorithm converges globally according to the probability equal to 1.


## 4.1. Marine Predator Algorithm

The marine predator algorithm [23] is a novel metaheuristic algorithm that simulates the foraging strategy of marine organisms. It updates the position of predator and prey by switching between the Lévy flight and the Brownian motion. Both predator and prey are solution candidates, since each prey in the ecosystem is a potential predator. Like most metaheuristic algorithms, MPA generates the initial population randomly in the search space. The mathematical model is described as follows:
(8)X0=randXmax−Xmin+Xmin,where rand is a random number that obeys a uniform distribution of 0,1 and Xmin and Xmax denote the lower and upper limits of the search space, respectively.In MPA, the location of the prey can be represented by a matrix. The corresponding mathematical expression is described as follows:(9)Prey=X1,1⋯X1,d⋮⋱⋮Xn,1⋯Xn,dn×d,where n denotes the number of populations and d is the dimension of the search space. Xi,j denotes the current position of the i-th prey in the j-th dimension. The individual with the best fitness value is called the top predator. It is used to construct the Elite matrix, whose size is the same as the Prey matrix. The mathematical expression of the Elite matrix is presented in the following:
(10)Elite=X1,1I⋯X1,dI⋮⋱⋮Xn,1I⋯Xn,dIn×d,where XI represents the top predator vector. The Elite matrix is obtained by replicating this vector n times. The searching procedure of the MPA algorithm can be separated into three stages in accordance with the speed variation of marine biological predation, which is described as follows.


### 4.1.1. High-Speed Ratio Phase

This phase occurs in the first third of the number of iterations, when the population is in the exploration phase. At this time, the speed of the predator is greater than that of the prey, and the mathematical model of prey position update is described as follows:(11)Whilet≤13tmax,i=1,2,⋯,n,(12)stepsizei=RB⊗Elitei−RB⊗Preyi,(13)Preyi=Preyi+P×R⊗stepsizei,where t and tmax denote the current and the maximum number of iterations, respectively. RB is a vector of random numbers obeying the normal distribution used to represent the Brownian motion. The symbol ⨂ denotes the entry-wise multiplication, R is a random number between 0 and 1, and P=0.5 is a constant.


### 4.1.2. Unit-Speed Ratio Phase

This phase occurs between one-third and two-thirds of the number of iterations and gradually indicates the population’s transition from exploration to exploitation. At this phase, the predator and the prey movement speed are equal. In this stage, the first half of the population is used for the exploration, and the movement strategy of the prey is the Lévy motion. The remaining half of the population would then be used for the exploitation where the movement strategy of the predator would be the Brownian motion. The mathematical model of this phase is described as follows:(14)While13tmax<t<23tmax,i=1,2,⋯,n/2,(15)stepsizei=RL⊗Elitei−RL⊗Preyi,(16)Preyi=Preyi+P×R⊗stepsizei,(17)While13tmax<t<23tmax,i=n/2,⋯,n,(18)stepsizei=RB⊗RB⊗Elitei−Preyi,(19)CF=1−ttmax2t/tmax,(20)Preyi=Elitei+P×CF⊗stepsizei,where RL is a vector of random numbers based on the Lévy distribution and CF is an adaptive parameter that controls the movement step of the predator.


### 4.1.3. Low-Speed Ratio Phase

A late iteration phase occurs in the last remaining third of the whole iteration. At this phase, the predator moves slower than the prey within a Lévy’s motion, which adopts the population into a local exploitation strategy. The mathematical model of this phase is described as follows:(21)Whilet≥23tmax,i=1,2,⋯,n,(22)stepsizei=RL⊗RL⊗Elitei−Preyi,(23)Preyi=Elitei+P×CF⊗stepsizei.


### 4.1.4. Eddy Formation and FADs Effect

The formation of eddies or Fish Aggregation Devices (FADs) affects the foraging behavior of the marine predators. This strategy mainly improves searching for individuals to escape the local optimal solutions. Its mathematical model is described as follows:(24)Preyi=Preyi+CFXmin+RL⊗Xmax−Xmin⊗Ur≤Pf,Preyi+Pf1−r+rPreyr1−Preyr2r>Pf,where U is a binary vector containing 0 and 1 and Pf=0.2 represents the probability of influence on the optimization process. r is a random number in the range 0,1, and the subscripts r1 and r2 represent the random index values of the prey matrix. In addition, the MPA algorithm has memory storage to remember the history of the best-reached position, which is similar to an essential greedy strategy. During the algorithm iteration, the fitness value of each individual in the current iteration is compared with the optimal solution of the previous iteration. The previous solution will then be replaced if the current’s fitness is better.


## 4.1.1. High-Speed Ratio Phase

This phase occurs in the first third of the number of iterations, when the population is in the exploration phase. At this time, the speed of the predator is greater than that of the prey, and the mathematical model of prey position update is described as follows:(11)Whilet≤13tmax,i=1,2,⋯,n,(12)stepsizei=RB⊗Elitei−RB⊗Preyi,(13)Preyi=Preyi+P×R⊗stepsizei,where t and tmax denote the current and the maximum number of iterations, respectively. RB is a vector of random numbers obeying the normal distribution used to represent the Brownian motion. The symbol ⨂ denotes the entry-wise multiplication, R is a random number between 0 and 1, and P=0.5 is a constant.


## 4.1.2. Unit-Speed Ratio Phase

This phase occurs between one-third and two-thirds of the number of iterations and gradually indicates the population’s transition from exploration to exploitation. At this phase, the predator and the prey movement speed are equal. In this stage, the first half of the population is used for the exploration, and the movement strategy of the prey is the Lévy motion. The remaining half of the population would then be used for the exploitation where the movement strategy of the predator would be the Brownian motion. The mathematical model of this phase is described as follows:(14)While13tmax<t<23tmax,i=1,2,⋯,n/2,(15)stepsizei=RL⊗Elitei−RL⊗Preyi,(16)Preyi=Preyi+P×R⊗stepsizei,(17)While13tmax<t<23tmax,i=n/2,⋯,n,(18)stepsizei=RB⊗RB⊗Elitei−Preyi,(19)CF=1−ttmax2t/tmax,(20)Preyi=Elitei+P×CF⊗stepsizei,where RL is a vector of random numbers based on the Lévy distribution and CF is an adaptive parameter that controls the movement step of the predator.


## 4.1.3. Low-Speed Ratio Phase

A late iteration phase occurs in the last remaining third of the whole iteration. At this phase, the predator moves slower than the prey within a Lévy’s motion, which adopts the population into a local exploitation strategy. The mathematical model of this phase is described as follows:(21)Whilet≥23tmax,i=1,2,⋯,n,(22)stepsizei=RL⊗RL⊗Elitei−Preyi,(23)Preyi=Elitei+P×CF⊗stepsizei.


## 4.1.4. Eddy Formation and FADs Effect

The formation of eddies or Fish Aggregation Devices (FADs) affects the foraging behavior of the marine predators. This strategy mainly improves searching for individuals to escape the local optimal solutions. Its mathematical model is described as follows:(24)Preyi=Preyi+CFXmin+RL⊗Xmax−Xmin⊗Ur≤Pf,Preyi+Pf1−r+rPreyr1−Preyr2r>Pf,where U is a binary vector containing 0 and 1 and Pf=0.2 represents the probability of influence on the optimization process. r is a random number in the range 0,1, and the subscripts r1 and r2 represent the random index values of the prey matrix. In addition, the MPA algorithm has memory storage to remember the history of the best-reached position, which is similar to an essential greedy strategy. During the algorithm iteration, the fitness value of each individual in the current iteration is compared with the optimal solution of the previous iteration. The previous solution will then be replaced if the current’s fitness is better.


## 4.2. Population Evolution Strategy Based on the Random Opposition-Based Learning and Differential Evolution Operator

The optimization process in each phase of MPA is highly depends on the prey matrix (described in section4.1). The MPA algorithm randomly initializes the population and then relies only on the Brownian or Lévy motion strategies to update the individual positions. The late iteration of the algorithm will reduce population diversity, narrow search range, reduce exploration ability, and tends to stagnate in the local optimum. Usually, in the optimization process of swarm intelligence algorithms, when the searched individual is close to the optimal solution, the population will quickly approach the direction of the optimal solution and accelerate the convergence of the algorithm. On the contrary, the convergence rate becomes slower.To solve this problem, opposition-based learning [42] (OBL) provides a search strategy that considers both the current solution and its corresponding inverse solution, expanding the search range and avoiding stagnation of the current solution. The research shows that the probability of the opposite solution of the population close to the global optimal solution is greater than that of the random solution. Therefore, many researchers have attempted to use OBL to improve both the quality of the population and the performance of the algorithm [43–48].To enhance the diversity of the population further, Long et al. [49] proposed a random opposition-based learning (ROL) strategy based on OBL, which could effectively improve the ability of the algorithm to jump out of the local optimum. Suppose a feasible solution of a d-dimensional space is X=x1,x2,⋯,xd, where x1,x2,⋯,xd∈R and xj∈lbj,ubj. The random opposite solution of X would be X^r=x^r,1,x^r,2,⋯,x^r,d and its calculation formula can be presented as follows:
(25)x^r,j=lbj+ubj−R1×xjj=1,2,⋯,d,where R1∈0,1 is a random number.Given the effectiveness of the ROL strategy for the improvement of algorithm performance and inspired by the mutation, crossover, and selection operators in the differential evolution algorithm (DE), this paper proposes a population evolution strategy based on ROL and differential evolution operator (RDES) to generate a higher quality population.First, the ROL population of the initial population is evaluated by Equation (25). Because the positions of all individuals in the ROL population may not be better than those of the initial individuals, both the initial and ROL populations are merged to expand the population searching range. Then, the differential evolution operator is considered based on the merged population to enrich the population’s diversity further. The specific operations are as follows:
(1)
Mutation Operation. Since the mutation operator in the DE has a significant influence on the algorithm, it is essential to choose an effective individual mutation method. This paper used the DE/current-to-best/1 mutation operator [29] to perform mutation operation on the merged population, as shown in the following formula:
(26)Vit=Xit+F⋅Xbestt−Xkt+Xmt−Xnt,where F is scaling factor, Xit denotes the i-th individual of the current generation, Xbestt denotes the optimal individual of the current population. Xkt, Xmt, and Xnt denote three different individuals randomly selected except the current individual, where the k,m,n∈1,2,⋯,2∗N are not equal to each other, and N represents the population size. Equation (26) reveals that the mutation operator uses the individuals’ local information and the population’s global optimal information to lead the population to find an optimal solution. It expands the diversity of the population and accelerates the searching speed of the algorithm simultaneously(2)
Crossover Operation. The original individual Xit and the new generated individual Vit are cross-operated after the differential mutation according to
(27)Ui,jt=Vi,jt,rand≤Pc,Xi,jt,otherwise,where Ui,jt represents the value of the j-th gene of the current i-th individual, rand is a random number uniformly distributed in 0,1, and Pc represents the crossover probability. The algorithm would keep the gene value of Vi,jt if rand<Pc, otherwise, it replaces the corresponding gene value with Xi,jt(3)
Selection of Operation. If the crossover operation results in a new individual with a better fitness value than the original individual, the original individual would be substituted. This process can be described as follows:
(28)Xit+1=Uit,fUit≤fXit,Xit,otherwise.Finally, the fitness values of all individuals in the merged population are calculated and sorted. The topN individuals with the best fitness value are selected to construct the prey matrix using the greedy strategy.Based on the above descriptions, the mainframe of the RDES is given in Algorithm1.Algorithm 1: Pseudocode of the RDES.
1:    Initialize the population withN individuals2:Fori≤Ndo3:     Using Eq. (25) to generate the ROL population4:End for5:    Merge initialization population and ROL population6:Fori≤2Ndo7:     Execute the mutation operation by Eq. (26)8:Ifrand<Pcthen9:      Execute the crossover operation by Eq. (27)10:End if11:    Using Eq. (28) to execute the selection operation12:End for13:    Select the topN fitness individuals from the mergedPopulation as the initial Prey matrix


## 4.3. Local Search Strategy Based on Improved SFLA

SFLA is a swarm intelligence algorithm proposed by Eusuff et al. [50]. They simulated the population’s mechanism of information sharing and communication during the frog foraging. The algorithm combines the advantages of the Memetic Algorithm (MA) [51] and the PSO [28]. Moreover, it benefits few control parameters, is an easy implementation, and has a strong global searching ability. Therefore, some researchers have hybridized SFLA with other algorithms to improve the optimization performance of the algorithm. For instance, in [52], an improved ant colony algorithm was combined with SFLA for solving the path planning problem, which accelerated the convergence of the algorithm. In [53], in order to improve the local search capability of the PSO, SFLA was integrated with PSO and used for the detection of sonar images.The MPA uses the FADs effect (Equation (24)) to support individuals’ searching for escape from the local optimum trap. This method lacks the guidance of population information to update the individuals’ position and limited the local development capability. Finding a more optimal solution becomes complicated once the population falls into the local optimum. The SFLA divides the whole population into several groups based on its unique grouping operator. The evolution of the individual is guided by the local optimal solution in the group and the global optimal solution of the population. It realizes the information sharing and communication between individuals. Therefore, the algorithm would not easily fall into the local optimum and would be more conducive to find the global optimal solution [22]. In this paper, by considering the grouping idea of SFLA, both MPA and SFLA algorithms are combined, and a local search algorithm based on the improved SFLA is proposed. The proposed local search algorithm replaces the FADs effect to improve the local development ability of MPA.The main idea of the SFLA is described in the following. First, all the population individuals are arranged in descending order according to their fitness values and divided intom groups. The first individual is assigned to the first group, the second to the second group, until the m-th to the m-th group. Them+1-th individual is then assigned to the first group, the m+2-th to the second group, and so on, until assigning the last one. Then, the positions of the worst individual Pw are updated according to Equations (29)–(31). The specific update method is as follows:
(29)Pnew=Pw+s×rand×Pb−Pw,(30)Pnew=Pw+s×rand×Pg−Pw,(31)Pnew=lb+rand×ub−lb,where Pnew denotes the updated position of the worst individual, Pb and Pw denotes the best and worst individuals in the group, randis a random number between 0,1, Pg represents the position of the globally optimal individual for the entire population, s is a constant representing the step size factor and that its value can be adjusted according to the problems, and ub and lb are the upper and lower bounds of the search space, respectively.Since the third location update method (as shown in Equation (31)) of the basic SFLA algorithm is relatively random, it lacks the guidance of population information, and consequently, the searching speed is slow. Therefore, this paper considers the usage of the local information of the group and the global information of the population to guide the algorithm for strengthening the ability of further search. A new location update method is proposed, as shown in
(32)Pnew=Pw+rand×Pg+Pb2−Pw.The mainframe of the improved local search strategy is shown in Algorithm2. First, the prey population is divided into m groups according to the grouping strategy which was introduced above. After the completion of grouping, the worst position of an individual would be updated using Equation (29). The original position is then replaced if the new position includes a better fitness value. Otherwise, the new position would be evaluated by Equation (30). If neither of the above two methods can generate a better position, the position of the worst individual in the group would be updated by Equation (32).Figure2 depicts the three updating methods for the position of the worst individual in each group. As it can be seen in Figure 2, Equation (29) means the best individual in the group guides the worst individual to move towards it. Equation (30) indicates that the worst individual approaches the global optimal individual, and Equation (32) represents that the population approaches the middle position between the global optimal and the intragroup optimal.Figure 2
Schematic diagram of the position update of the worst individual (the black dashed line represents the guidance direction of the population information, and the solid blue line represents the trajectory of the worst individual position update).The improved local search method transmits the information according to the population grouping. It effectively combines the global information of the population and the local search. It also improves the exploitation ability of the population. In addition, the position of the worst individuals in the group is only updated each time after grouping, which reduces the computation time and increases the convergence rate of the algorithm. The fitness value of the population is inevitably improved accordingly when the individuals of all groups are remixed into the population. Therefore, it leads the population towards the optimal solution.Algorithm 2: Pseudocode of the local search strategy.
1:  Initialization parameters: The number of groupsm; the number of iterations within the group tm;2:Fori=1:mdo3:Forj=1:tmdo4:     RecordPw,Pb, Pg5:     Update the worst Prey using Eq. (29)6:IffitnessPnew<fitnessPwthen7:Pw=Pnew8:Else9:        Update the worst Prey using Eq. (30)10:fitnessPnew<fitnessPwthen11:Pw=Pnew12:Else13:         Update the worst Prey using Eq. (32)14:End if15:End if16:End for17:End for


## 4.4. Dynamic Learning Mechanism

Ergezer et al. proposed a variant of OBL called quasireflected opposition-based learning (QOBL) [54]. They proved that the probability of the QOBL solution is closer to the optimal solution than the general OBL solution is 11/16. This paper proposes a QOBL strategy for the prey population to improve the optimization accuracy and accelerate the convergence speed of the algorithm further. The QOBL position (Preyiqr) of the prey in the generation can be described as follows:
(33)Preyiqrt=randPreyi,ci,(34)ci=minPreyi+maxPreyi2,where Preyi represents the position of the i-th prey and ci represents the midpoint of the dynamic boundary. Unlike the fixed boundary, the dynamic search boundary adjusts the maximum and minimum values in all dimensions of the current i-th individual as the upper and lower bounds of the search space. This method can further reduce the possibility of stagnating the population in the local optima, enhance the exploitation capability of the algorithm, and increase the convergence speed [55].Finally, the elite strategy is employed to retain the individuals who carry better fitness value. If the fitness value of the QOBL position’s fitness value is superior to the current individual, the new position would be assigned to the current individual. The strategy is modeled as follows:(35)Preyit+1=PreyiqrtiffPreyiqrt<fPreyit,Preyitotherwise.Figure3 shows the process of applying the dynamic learning mechanism based on QOBL to the MPA algorithm. First, the QOBL population of the prey population is evaluated. Then, the fitness values of both initial prey and the QOBL populations would be calculated separately. Finally, the individuals with better fitness values are retained for constructing new prey populations according to the Equation (35).Figure 3
Schematic diagram of population selection based on QOBL.


## 4.5. IMPA Algorithm

In this paper, a multistrategy integrated improved marine predator algorithm (IMPA) is proposed to increase both accuracy of solution and convergence speed and prevent falling into the local optimum of MPA. The pseudocode of IMPA is given in Algorithm3, and the main framework of the algorithm is demonstrated in Figure 4.Algorithm 3: The proposed IMPA.
Input: Number of population N, the maximum number of iterations tmax, the crossover probability Pc, the number of groups m, the number of iterations within the group tm, the dimensions of each individual d, and search range lb,ubOutput: Optimal fitness1:  Initialize the population of Preys withN individuals2:Whilet≤tmaxdo3:     Execute RDES to generate a high-quality population (Algorithm1)4:     Calculate the fitness and build the Elite matrix5:Ift≤tmax/3then6:      Update Prey using Eq. (13)7:Else iftmax/3<t<2∗tmax/38:        For the first half of the population (i=1,2,⋯,N/2)9:        Update Prey using Eq. (16)10:        For the remained half of the population (i=N/2,⋯,N)11:        Update Prey using Eq. (20)12:Else ift≥2∗tmax/313:        Update Prey using Eq. (23)14:End if15:     Complete memory saving and Elite update16:     Execute the local search strategy (Algorithm2)17:Fori=1:Ndo18:       Calculate the QOBL position according to Eq. (33)19:       Select the better individual according to Eq. (35)20:End for21:  Complete memory saving and Elite updating22:  Next iterationt=t+125:End whileFigure 4
IMPA flow chart.


## 4.6. Algorithm Validity Test

In this section, the Schwefel2.26 Function is selected to test the validity of the three improved strategies. The population distribution maps of IMPA and MPA algorithms are drawn. The initial and final distributions of both algorithms are shown in Figures5 and 6 while the population size and the maximum number of iterations are 50 and 20, respectively.Figure 5
Population distribution of MPA.
(a)
Initial population distribution of MPA(b)
Final population distribution of MPAFigure 6
Population distribution of IMPA.
(a)
Initial population distribution of IMPA(b)
Final population distribution of IMPAFigures5(b) and 6(b) show that the population distribution of MPA is relatively scattered after 20 iterations and the convergence speed is low. Moreover, most individuals gather near the local optimum, which means that the MPA would not easily jump out of the local optimum when dealing with multidimensional complex functions. In contrast, all individuals in the IMPA converge to the global optimum quickly, proving the effectiveness of the proposed strategies.


## 4.7. Time Complexity Analysis

Time complexity is an essential indicator for predicting the performance of an algorithm. Suppose thatN is the size of populations, D is the dimension of individuals, and T is the maximum number of iterations. The time complexity of the MPA is evaluated ON∙D∙T. The proposed IMPA algorithm would not change the structure of individual position updates in the three phases of the basic MPA algorithm. The analysis of the time complexity is explained in the following.The increased complexity of the ROL population while performing the RDES (Algorithm1) is ON, where this evaluation becomes O2∙N∙D while performing the DE evolution operation. While executing the local search strategy (Algorithm 2), the evaluated complexity is related to the number of the groups m and the number of iterations tm, where m≤N/2. Therefore, in the worst case, the increased complexity of this part would be ON. The increased complexity while computing the QOBL population would be ON. Finally, the evaluated time complexity of IMPA would be OTN+2∙N∙D+N+N, which results in ON∙D∙T. It concluded that the improved strategies proposed in this study would not increase the time complexity of the proposed IMPA.


## 4.8. Global Convergence Analysis of IMPA

The convergence analysis of the three phases within the optimization process for the basic MPA algorithm has been extensively discussed [23]. It approved that the MPA is a convergent algorithm. The conclusion of the mentioned study is reviewed in the following.Assuming thatPreyt is a random sampling sequence of MPA, the MPA eventually converges to an Elite matrix (globally optimal position) consisting of the top predators, which is limt⟶∞Preyt=Elite.The global convergence of the hybrid optimization algorithm was analyzed, and several sufficient conditions for the global convergence of the hybrid optimization algorithm were approved [56]. The approved relevant theorems of the mentioned research are represented below.Lemma 1.
fn is a monotonically nonincreasing sequence with a lower limited bound in which limt⟶∞fn=inffn [56].Lemma 2.
If an algorithm satisfies both following conditions (properties), the optimal sequence will converge globally according to a probability equal to 1 [56].Condition 3.
The greedy strategy is used to retain the elite while comparing with the historical optimal solution during the iteration. The mentioned strategy is represented in the following:(36)ft+1best=fX,fX<ftbest,ftbest,otherwise,where ftbest is the optimal fitness value found after t times of comparisons which retained to the moment t, and X is the solution used for t+1-th comparison withftbest.Condition 4.
The probability of transitioning from any nonglobal optimal point (might be set toX∗) to the corresponding level set LX∗=XfX<fX∗,X∈S is not 0, where S represents the entire search space.Theorem 5.
The proposed IMPA algorithm converges globally according to the probability equal to 1.
Assume that the objective functionf of the IMPA algorithm optimization is set to find the minimum value, fn is the random sampling sequence, ftbest is the optimal sequence of sampling points, and fgbest is the global optimal solution. As shown below, Theorem 5 can be proved by Lemma 1 and Lemma 2.Proof.
Because the three proposed improved strategies and the marine memory storage procedure in the IMPA algorithm all use the greedy strategy to retain the optimal solution of the population, the IMPA algorithm would satisfy Condition3.According to Lemma1, the limit of the optimal sequence ftbest exists. Let limt⟶∞ftbest=fnbest, and fnbest≠fgbest. Because fgbest is the global optimal solution, therefore fnbest>fgbest. Since limt⟶∞ftbest=infftbest=fnbest, it follows that for ∀t∈Z+, ftbest≥fnbest. Let X∗ be the sampling point corresponding to fnbest, so, fnbest=fX∗. X∗ is obviously a nonglobal optimal solution. Because the IMPA uses the same update operation as well as the MPA algorithm in the high-, unit-, and low-speed ratio phases, the convergence of the population sequence of IMPA would be consistent with that of the MPA and will eventually converge to the Elite matrix, limn⟶∞fn=Elite=fgbest. That is, the nonglobal optimal point X∗ can converge to the Elite matrix (global optimal point). Therefore, when the X assigns value to a sufficiently small neighborhood around the elite matrix, PLX∗=PXfX<fX∗,X∈S≠0, the Condition 4 is satisfied.In summary, the IMPA algorithm satisfies both conditions of Lemma2. Therefore, the IMPA algorithm converges globally according to the probability equal to 1.


## 5. Function Test and Result Analysis

In this study, 15 benchmark functions are selected for testing to verify the performance of the proposed IMPA algorithm. The results are compared with the seven algorithms of MSIMPA, MPA, ESSA, IGWO, WOA, PSO, and DE. The information of all the 15 test functions are shown in Table1, where F1-F5 are complex unimodal functions, F6-F10 are multidimensional multimodal functions, and F11-F15 are fixed-dimensional multimodal functions. The population size of all algorithms and the maximum number of iterations are set to 30 and 500, respectively. The multidimensional functions (F1-F10) are tested in 30 and 100 dimensions, respectively. The main parameters of each algorithm are set as shown in Table 2. Each of the six algorithms is run 30 times independently. The best fitness value (Best), the average best fitness value (Ave), the standard deviation of the best fitness value (Std) results, and the average running time are recorded. To intuitively observe the performance difference of each algorithm, the algorithms are ranked according to the average best fitness value. If the average best fitness values are equal, the algorithms would then be ranked according to the best standard deviation value.Table 1
Benchmark test functions.
FunctionDimIntervalMinF1x=∑i=1nxi230/100[-100, 100]0F2x=∑i=1nxi+∏i=1nxi30/100[-10, 10]0F3x=∑i=1n−1100xi+1−xi22+xi−1230/100[-30, 30]0F4x=∑i=1nxi+0.5230/100[-100, 100]0F5x=maxixi,1≤i≤n30/100[-100, 100]0F6x=−20exp−0.21n∑i=1nxi2−exp1n∑i=1ncos2πxi+20+e30/100[-32, 32]0F7x=14000∑i=1nxi2−∏i=1ncosxii+130/100[-600, 600]0F8x=∑i=1nixi4+random0,130/100[-1.28, 1.28]0F9x=∑i=1n−xisinxi30/100[-500, 500]-418.98nF10x=π/n10sinπy1+∑i=1n−1yi−121+10sin2πyi+1+yn−12+∑i=1nuxi,10,100,4yi=1+xi+1/4uxi,a,k,m=kxi−amxi>a0−a<xi<ak−xi−amxi<−a30/100[-50, 50]0F11x=x1+2x2−72+2x1+x2−522[-10, 10]0F12x=∑i=17X−aiX−aiT+ci−14[0, 10]-10.4029F13x=−∑i=14ciexp−∑j=13aijxj−pij23[0, 1]-3.8628F14x=1.5−x1+x1x22+2.25−x1+x1x222+2.625−x1+x1x2322[-4.5, 4.5]0F15=sin23πx1+x1−121+sin23πx2+x2−121+sin22πx22[-10, 10]0Table 2
Parameter Settings for the comparing optimization algorithms.
AlgorithmParametersIMPAP=0.5,Pc=0.5, m=8, tm=3, s=2MSIMPAFADs=0.2, P=0.5MPAFADs=0.2, P=0.5ESSAST=0.8, SD=0.2IGWOa=2⟶0WOAa=2⟶0, b=1PSOω=0.5, c1=1.3, c2=2.8DEPc=0.8, F=0.5


### 5.1. Numerical Analysis

Tables3 and 4 revealed that the IMPA algorithm has a significant advantage over other algorithms in the average best fitness and standard deviation values within the operation of both 30 and 100 dimensions. The IMPA algorithm ranks the best optimization performance over the other seven algorithms. Table 3 revealed that the IMPA algorithm could evaluate the theoretical optimum for F1-F2, F4-F5, F7, and F11-F14, the algorithm shows a strong global search capability. Although the IMPA could not assess the theoretical optimum for the F6, F8, F9-F10, and F15, it achieves the highest optimization accuracy over the other seven algorithms. None of the six algorithms in F3 could evaluate the optimal value, but the ESSA achieves the highest optimization accuracy, exceeding five orders of magnitude. The optimization accuracy and stability of the MSIMPA algorithm, as well as the IMPA, show good optimization performance on the multiple functions. The optimization performance of the MPA in F13 ranks the best with good optimization ability and stability. However, the IMPA and MPA show a slight difference in the terms of standard deviation value. It indicates that the IMPA could not reduce the optimization capability of the algorithm. The IMPA can still maintain a high solution accuracy in the 100-dimensional test. For instance, all the six algorithms fall into the local optimum in F10. However, the average optimization accuracy of the IMPA is at least 24 times greater than the magnitude of other algorithms. It shows a stronger capability to escape the local optimal. Moreover, it proves the effectiveness of the local search strategy based on the improved SFLA.Table 3
Comparison of benchmark function optimization results (30-dimension and fixed-dimension).
FunctionAlgorithmBestAveStdRankRunningtime (s)F1IMPA00010.2197MSIMPA00010.4990MPA3.249e-255.6984e-236.8011e-2360.1168ESSA3.626e-2731.0578e-285.2877e-2840.0929IGWO3.6013e-302.0879e-283.8839e-2850.3297WOA7.8706e-851.2859e-724.9154e-7230.0208PSO8.3500e-083.9739e-021.4805e-0170.1687DE9.1636e-051.4571e-017.6655e-0180.0705F2IMPA00010.2553MSIMPA00010.5151MPA2.0622e-142.1332e-131.9902e-1360.1194ESSA05.0428e-172.3574e-1650.0940IGWO1.6695e-188.3521e-185.7087e-1840.3340WOA2.2386e-575.2527e-511.755e-5030.0226PSO2.9933e-037.9203e-021.6586e-0180.1731DE3.6882e-031.0604e-026.1972e-0370.0716F3IMPA2.2747e+012.3184e+012.0370e-0130.2810MSIMPA2.3646e+012.4331e+012.9350e-0140.6913MPA2.4583e+012.5351e+014.3584e-0160.1373ESSA2.2796e-192.4733e-054.5664e-0510.1161IGWO2.3720e+012.4335e+018.5325e-0150.3542WOA2.6976e+012.7968e+014.7406e-0170.0312PSO2.6555e-054.1230e-011.3376e+0020.1801DE2.6727e+017.7838e+015.1636e+0180.0808F4IMPA00010.2494MSIMPA1.2009e-097.5903e-095.2053e-0930.4986MPA1.4150e-081.2359e-046.4591e-0440.1213ESSA2.0396e-133.0063e-094.9969e-0920.1020IGWO3.3595e-053.7862e-028.6222e-0260.3391WOA8.0338e-023.8291e-012.2480e-0180.0215PSO9.3274e-063.0549e-011.2967e+0070.1721DE1.0879e-045.0759e-043.1912e-0450.0715F5IMPA00010.2440MSIMPA00010.4935MPA5.3530e-103.9982e-092.4312e-0940.1151ESSA02.7442e-141.4981e-1330.0934IGWO7.7868e-071.1347e-051.2372e-0550.3325WOA4.4745e+004.1640e+012.7544e+0180.0209PSO2.5094e-045.8164e-022.6117e-0160.1698DE5.7486e+001.5239e+016.1165e+0070.0707F6IMPA8.8818e-168.8818e-16010.2431MSIMPA8.8818e-168.8818e-16010.5607MPA1.7852e-131.3499e-128.3401e-1360.1242ESSA8.8818e-161.1428e-144.6812e-1440.1093IGWO4.3521e-146.4126e-141.2392e-1450.3395WOA8.8818e-164.6777e-152.4567e-1530.0254PSO4.1938e-044.3581e-021.2346e-0180.1841DE2.8525e-036.6528e-033.7812e-0370.0804F7IMPA00010.2934MSIMPA00010.7127MPA00010.1413ESSA00010.1215IGWO03.7034e-035.9281e-0360.3616WOA03.7007e-182.0270e-1750.0331PSO6.5533e-064.0149e+007.8970e+0080.1906DE3.7623e-047.4684e-032.1069e-0270.0865F8IMPA9.4248e-061.0125e-056.7652e-0510.4376MSIMPA3.4548e-067.6693e-057.7777e-0521.6014MPA3.6961e-041.2423e-035.8351e-0440.2326ESSA2.3136e-056.1201e-045.7285e-0430.1838IGWO9.4534e-042.5539e-031.3054e-0350.4462WOA2.4683e-052.5666e-032.9612e-0360.0752PSO5.9669e-034.8591e-023.7271e-0280.2317DE1.6948e-023.9809e-021.5271e-0270.1256F9IMPA-11977.2949-10958.0172374.261420.2557MSIMPA-11412.6515-10955.865297.957930.7244MPA-9904.4343-8861.0457465.264640.1415ESSA-11029.4511-8172.11761117.562760.1118IGWO-11351.2906-7855.77351753.144370.3589WOA-12567.8533-10576.35021797.365840.0332PSO-12569.4866-12066.06151024.538310.1825DE-10862.8222-8053.88891609.651980.0859F10IMPA1.5705e-321.5705e-325.5674e-4810.7041MSIMPA3.4698e-101.3188e-099.3011e-1031.3146MPA1.5317e-094.1000e-052.2287e-0450.4124ESSA4.9661e-152.4918e-104.1901e-1020.2997IGWO2.7700e-062.9326e-041.2520e-0360.6210WOA6.1646e-032.3879e-022.0913e-0270.1618PSO6.2733e-081.8325e-052.6361e-0540.3174DE2.1047e-052.1863e-012.9378e-0180.2135F11IMPA00010.1037MSIMPA00010.1860MPA09.1953e-253.0168e-2450.0553ESSA02.6295e-321.4403e-3140.0869IGWO1.4018e-271.5287e-225.675e-2260.2465WOA5.4848e-069.8198e-041.0474e-0370.0146PSO3.0552e-044.1497e-018.0658e-0180.1467DE00010.0560F12IMPA-10.4029-10.40291.3194e-1610.2139MSIMPA-10.4029-10.40295.9774e-1330.3531MPA-10.4029-10.40292.3540e-1140.1048ESSA-10.4029-10.40297.3759e-1620.1252IGWO-10.4029-10.40292.3253e-0850.3106WOA-10.4029-7.43263.0665e+0080.0338PSO-10.4028-10.31542.6727e-0160.1729DE-10.4029-10.00451.5274e+0070.0771F13IMPA-3.8628-3.86282.6662e-1540.1102MSIMPA-3.8628-3.86282.6684e-1550.2435MPA-3.8628-3.86282.3715e-1510.0765ESSA-3.8628-3.86282.4184e-1520.1033IGWO-3.8628-3.86282.4795e-1530.2701WOA-3.8628-3.85820.006503270.0233PSO-3.862-3.84654.6015e-0280.1550DE-3.8628-3.86282.7101e-1560.0651F14IMPA00010.1241MSIMPA03.5595e-301.3958e-2930.1955MPA3.4592e-283.0451e-221.5159e-2150.0588ESSA02.8995e-241.1012e-2340.0891IGWO7.4237e-217.1524e-171.4177e-1660.2510WOA1.4198e-131.0161e-012.6348e-0180.0161PSO1.2304e-041.3578e-039.7707e-0470.1463DE00010.0575F15IMPA1.3498e-311.3498e-316.6809e-4710.1196MSIMPA1.3498e-311.3498e-316.6809e-4710.1903MPA7.2938e-276.392e-231.3055e-2260.0567ESSA1.3498e-311.3498e-316.6809e-4710.0890IGWO7.8004e-306.2979e-262.7686e-2550.2485WOA2.2215e-098.1541e-061.6292e-0570.0154PSO1.2427e-071.6901e-032.9564e-0380.1472DE1.3498e-311.3498e-316.6809e-4710.0561Table 4
Comparison of benchmark function optimization results (100-dimension).
FunctionAlgorithmBestAveStdRankRunningtime (s)F1IMPA00010.4799MSIMPA00011.7578MPA2.0936e-202.0131e-191.8807e-1950.2724ESSA1.3162e-2157.4533e-314.0189e-3040.1457IGWO3.4070e-132.7008e-122.0152e-1260.5957WOA4.056e-821.6867e-718.8984e-7130.0481PSO2.5648e-063.4789e-011.4948e+0070.2561DE6.8945e+021.2783e+035.3839e+0280.1209F2IMPA02.5048e-297020.5490MSIMPA00011.9195MPA1.2682e-131.418e-111.0044e-1150.2828ESSA4.6552e-622.6374e-141.4083e-1340.1475IGWO7.0279e-091.5616e-085.995e-0960.6088WOA5.5313e-612.3061e-499.5159e-4930.0489PSO4.8680e-033.9462e-014.8301e-0170.2833DE1.2305e+011.9098e+013.7135e+0080.1235F3IMPA9.2984e+019.3322e+011.7465e-0130.6006MSIMPA9.4100e+019.4491e+011.7690e-0142.4819MPA9.5854e+019.7101e+017.8353e-0150.3002ESSA8.5565e-093.6629e-045.9162e-0410.1675IGWO9.4975e+019.7472e+011.1160e+0060.6292WOA9.7476e+019.8149e+012.5179e-0170.0606PSO6.1575e-047.9545e+014.3079e+0220.2824DE3.7089e+051.5227e+061.0087e+0680.1382F4IMPA00010.5204MSIMPA7.3420e-025.7584e-013.4773e-0131.6146MPA3.0751e+003.9660e+006.4026e-0140.2583ESSA8.1256e-116.5052e-071.3268e-0620.1420IGWO5.4614e+007.8873e+001.0375e+0060.5740WOA1.6505e+004.1532e+009.5339e-0150.0437PSO3.0226e-073.4222e+011.3481e+0270.2539DE5.8683e+021.0714e+033.8677e+0280.1172F5IMPA00010.4683MSIMPA00011.4869MPA6.3693e-082.25e-071.1712e-0740.2249ESSA4.8795e-814.7458e-152.3905e-1430.1178IGWO1.1614e+004.6770e+002.8552e+0060.5037WOA2.3009e+018.6982e+011.6871e+0170.0388PSO3.4590e-049.5632e-031.0224e-0250.2272DE8.3771e+018.7385e+011.4535e+0080.1050F6IMPA8.8818e-168.8818e-16010.5018MSIMPA8.8818e-168.8818e-16011.8198MPA1.2297e-115.5801e-114.2214e-1150.2448ESSA8.8818e-161.4803e-151.3467e-1530.1498IGWO6.1236e-081.7033e-076.0547e-0860.5292WOA8.8818e-164.6777e-152.6279e-1540.0486PSO3.1384e-045.8055e-021.8421e-0170.2642DE4.7890e+006.2446e+008.0832e-0180.1407F7IMPA00010.5311MSIMPA00012.6065MPA00010.2776ESSA00010.1679IGWO2.7289e-132.9180e-036.3530e-0360.5641WOA00010.0631PSO7.2014e-071.4439e+014.1840e+0180.2778DE5.0666e+001.3019e+016.3009e+0070.1553F8IMPA3.0053e-061.8797e-051.4073e-0410.9748MSIMPA2.2357e-067.5344e-056.8927e-0522.0507MPA3.8496e-041.8048e-035.9696e-0440.5977ESSA3.3627e-057.2115e-047.7153e-0430.3942IGWO0.00542031.5155e-024.5389e-0360.8735WOA2.6318e-055.3274e-035.2580e-0350.2172PSO4.9434e-036.7585e-024.1483e-0270.4266DE1.4769e+003.0316e+001.0210e+0080.2811F9IMPA-40708.7689-33742.4422302.829130.5864MSIMPA-35220.1185-33533.4499895.514741.9733MPA-26300.7559-24710.3869907.327550.2990ESSA-28921.6901-23487.29372711.124460.1623IGWO-26070.9887-16182.55656433.380070.5971WOA-41895.7857-34167.74925693.474420.0678PSO-41898.2886-40509.56263560.758910.2540DE-10149.0836-60561.65082047.619980.1443F10IMPA4.7116e-334.7116e-331.3918e-4811.8984MSIMPA2.1453e-036.4522e-033.3446e-0343.7580MPA3.1036e-024.6095e-028.6771e-0351.0844ESSA1.0014e-154.9800e-091.0483e-0820.7149IGWO1.3000e-012.2248e-015.3100e-0271.4428WOA1.9735e-025.3398e-022.6769e-0260.4637PSO6.4788e-092.0935e-056.6560e-0530.6519DE1.4390e+044.3132e+055.7432e+0580.5338


### 5.2. Convergence Analysis

In this paper, the average convergence curves of the algorithms are selected to evaluate the convergence ability of the optimization algorithms. Figure7 shows each algorithm’s 30-dimensional average convergence curves on the 15 benchmark test functions. The IMPA has excellent convergence speed and convergence accuracy on the remaining 14 test functions except for the F3. In F1, F2, F5, and F7, both IMPA and MSIMPA converge to the theoretical optimum quickly. In F12 and F13, the optimization precision and convergence speed of the IMPA improved a little over the other seven algorithms. It is worth mentioning that the DE can quickly converge to the global optimal value in F11, F14 and F15, which shows excellent optimization ability. In addition, it can be seen from Tables 3 and 4 that the running time of the proposed IMPA algorithm on most test functions is about twice that of MPA, but the optimization accuracy of IMPA has been greatly improved. Therefore, it is worth sacrificing a little time to obtain a higher quality of the solution. IMPA can solve the issues of low optimization accuracy and low convergence rate of MPA and shows significant improvements in the convergence speed and solution accuracy. These improvements are due to the implementation of the population evolution strategy and the dynamic learning mechanism, which effectively enhance the algorithm’s ability for the global searching while enriching the diversity of the population.Figure 7
Comparison of convergence curves.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)


### 5.3. The Wilcoxon Rank Sum Test

Only relying on the comparison of mean and standard deviation values, the superiority of the IMPA cannot be approved entirely [57]. Therefore, it needs to be tested at the statistical level. In this paper, to distinguish the significant differences in the optimization performance of the IMPA over the other seven algorithms, Wilcoxon statistical tests are performed. Test results of 30-dimension and fixed-dimension on the benchmark functions at the significant level of p=5% are obtained. The test results are shown in Table 5. When p<5%, the zero hypothesis considers a rejection, indicating that there is a significant difference between both algorithms. When p>5%, it indicates that the difference between both algorithms is not substantial. Since the best algorithm on a benchmark function cannot be compared with itself, the best algorithm on each benchmark function is marked as N/A, indicating that the performance is comparable and optimal between the two. “S” represents the saliency judgment result, and the symbols “+”, “-”, and “=,” respectively, indicate that the performance of IMPA is better than, worse than, and equivalent to the comparison algorithm. As shown in Table 5, most of the p values are less than 5%, indicating that the IMPA algorithm is significantly different in most of the test functions among the other seven algorithms. For F1, F2, F5-F7, F11, and F15, the performance of IMPA and MSIMPA are comparable, and all can obtain the best value. In general, the IMPA algorithm shows significant advantages over the other seven algorithms.Table 5
Thep value of the Wilcoxon rank sum test on 15 benchmark functions (30-dimension and fixed-dimension).
FunctionMSIMPAMPAESSAIGWOWOAPSODEp value Sp value Sp value Sp value Sp value Sp value Sp value SF1N/A=1.212e-12+1.212e-12+1.212e-12+1.212e-12+1.212e-12+1.212e-12+F2N/A=1.212e-12+4.574e-12+1.212e-12+1.212e-12+1.212e-12+1.212e-12+F33.020e-11+3.020e-11+3.020e-11+3.020e-11+3.020e-11+5.573e-10+3.020e-11+F41.212e-12+1.212e-12+1.212e-12+1.212e-12+1.212e-12+1.212e-12+1.212e-12+F5N/A=1.212e-12+4.574e-12+1.212e-12+1.212e-12+1.212e-12+1.212e-12+F6N/A=1.212e-12+4.191e-02+1.031e-12+1.163e-09+1.212e-12+1.212e-12+F7N/A=N/A=N/A=6.617e-04+3.337e-01-1.212e-12+1.212e-12+F84.515e-02+3.020e-11+3.081e-08+3.020e-11+7.739e-06+3.020e-11+3.020e-11+F99.117e-01-3.020e-11+1.777e-10+6.722e-10+9.941e-01-3.256e-07+3.020e-11+F101.212e-12+1.212e-12+1.212e-12+1.212e-12+1.212e-12+1.212e-12+1.212e-12+F11N/A=4.574e-12+3.337e-01-1.212e-12+1.212e-12+1.212e-12+N/A=F122.758e-11+1.449e-11+8.435e-06+1.449e-11+1.449e-11+1.136e-11+7.338e-03+F133.128e-01-9.452e-09+3.568e-07+2.119e-05+1.720e-12+1.720e-12+3.337e-03+F142.19e-06+1.212e-12+1.657e-11+1.212e-12+1.212e-12+1.212e-12+N/A=F15N/A=1.212e-12+N/A1.212e-12+1.212e-12+1.212e-12+N/A=+/=/-6/7/214/1/013/1/115/0/013/0/215/0/012/3/0


## 5.1. Numerical Analysis

Tables3 and 4 revealed that the IMPA algorithm has a significant advantage over other algorithms in the average best fitness and standard deviation values within the operation of both 30 and 100 dimensions. The IMPA algorithm ranks the best optimization performance over the other seven algorithms. Table 3 revealed that the IMPA algorithm could evaluate the theoretical optimum for F1-F2, F4-F5, F7, and F11-F14, the algorithm shows a strong global search capability. Although the IMPA could not assess the theoretical optimum for the F6, F8, F9-F10, and F15, it achieves the highest optimization accuracy over the other seven algorithms. None of the six algorithms in F3 could evaluate the optimal value, but the ESSA achieves the highest optimization accuracy, exceeding five orders of magnitude. The optimization accuracy and stability of the MSIMPA algorithm, as well as the IMPA, show good optimization performance on the multiple functions. The optimization performance of the MPA in F13 ranks the best with good optimization ability and stability. However, the IMPA and MPA show a slight difference in the terms of standard deviation value. It indicates that the IMPA could not reduce the optimization capability of the algorithm. The IMPA can still maintain a high solution accuracy in the 100-dimensional test. For instance, all the six algorithms fall into the local optimum in F10. However, the average optimization accuracy of the IMPA is at least 24 times greater than the magnitude of other algorithms. It shows a stronger capability to escape the local optimal. Moreover, it proves the effectiveness of the local search strategy based on the improved SFLA.Table 3
Comparison of benchmark function optimization results (30-dimension and fixed-dimension).
FunctionAlgorithmBestAveStdRankRunningtime (s)F1IMPA00010.2197MSIMPA00010.4990MPA3.249e-255.6984e-236.8011e-2360.1168ESSA3.626e-2731.0578e-285.2877e-2840.0929IGWO3.6013e-302.0879e-283.8839e-2850.3297WOA7.8706e-851.2859e-724.9154e-7230.0208PSO8.3500e-083.9739e-021.4805e-0170.1687DE9.1636e-051.4571e-017.6655e-0180.0705F2IMPA00010.2553MSIMPA00010.5151MPA2.0622e-142.1332e-131.9902e-1360.1194ESSA05.0428e-172.3574e-1650.0940IGWO1.6695e-188.3521e-185.7087e-1840.3340WOA2.2386e-575.2527e-511.755e-5030.0226PSO2.9933e-037.9203e-021.6586e-0180.1731DE3.6882e-031.0604e-026.1972e-0370.0716F3IMPA2.2747e+012.3184e+012.0370e-0130.2810MSIMPA2.3646e+012.4331e+012.9350e-0140.6913MPA2.4583e+012.5351e+014.3584e-0160.1373ESSA2.2796e-192.4733e-054.5664e-0510.1161IGWO2.3720e+012.4335e+018.5325e-0150.3542WOA2.6976e+012.7968e+014.7406e-0170.0312PSO2.6555e-054.1230e-011.3376e+0020.1801DE2.6727e+017.7838e+015.1636e+0180.0808F4IMPA00010.2494MSIMPA1.2009e-097.5903e-095.2053e-0930.4986MPA1.4150e-081.2359e-046.4591e-0440.1213ESSA2.0396e-133.0063e-094.9969e-0920.1020IGWO3.3595e-053.7862e-028.6222e-0260.3391WOA8.0338e-023.8291e-012.2480e-0180.0215PSO9.3274e-063.0549e-011.2967e+0070.1721DE1.0879e-045.0759e-043.1912e-0450.0715F5IMPA00010.2440MSIMPA00010.4935MPA5.3530e-103.9982e-092.4312e-0940.1151ESSA02.7442e-141.4981e-1330.0934IGWO7.7868e-071.1347e-051.2372e-0550.3325WOA4.4745e+004.1640e+012.7544e+0180.0209PSO2.5094e-045.8164e-022.6117e-0160.1698DE5.7486e+001.5239e+016.1165e+0070.0707F6IMPA8.8818e-168.8818e-16010.2431MSIMPA8.8818e-168.8818e-16010.5607MPA1.7852e-131.3499e-128.3401e-1360.1242ESSA8.8818e-161.1428e-144.6812e-1440.1093IGWO4.3521e-146.4126e-141.2392e-1450.3395WOA8.8818e-164.6777e-152.4567e-1530.0254PSO4.1938e-044.3581e-021.2346e-0180.1841DE2.8525e-036.6528e-033.7812e-0370.0804F7IMPA00010.2934MSIMPA00010.7127MPA00010.1413ESSA00010.1215IGWO03.7034e-035.9281e-0360.3616WOA03.7007e-182.0270e-1750.0331PSO6.5533e-064.0149e+007.8970e+0080.1906DE3.7623e-047.4684e-032.1069e-0270.0865F8IMPA9.4248e-061.0125e-056.7652e-0510.4376MSIMPA3.4548e-067.6693e-057.7777e-0521.6014MPA3.6961e-041.2423e-035.8351e-0440.2326ESSA2.3136e-056.1201e-045.7285e-0430.1838IGWO9.4534e-042.5539e-031.3054e-0350.4462WOA2.4683e-052.5666e-032.9612e-0360.0752PSO5.9669e-034.8591e-023.7271e-0280.2317DE1.6948e-023.9809e-021.5271e-0270.1256F9IMPA-11977.2949-10958.0172374.261420.2557MSIMPA-11412.6515-10955.865297.957930.7244MPA-9904.4343-8861.0457465.264640.1415ESSA-11029.4511-8172.11761117.562760.1118IGWO-11351.2906-7855.77351753.144370.3589WOA-12567.8533-10576.35021797.365840.0332PSO-12569.4866-12066.06151024.538310.1825DE-10862.8222-8053.88891609.651980.0859F10IMPA1.5705e-321.5705e-325.5674e-4810.7041MSIMPA3.4698e-101.3188e-099.3011e-1031.3146MPA1.5317e-094.1000e-052.2287e-0450.4124ESSA4.9661e-152.4918e-104.1901e-1020.2997IGWO2.7700e-062.9326e-041.2520e-0360.6210WOA6.1646e-032.3879e-022.0913e-0270.1618PSO6.2733e-081.8325e-052.6361e-0540.3174DE2.1047e-052.1863e-012.9378e-0180.2135F11IMPA00010.1037MSIMPA00010.1860MPA09.1953e-253.0168e-2450.0553ESSA02.6295e-321.4403e-3140.0869IGWO1.4018e-271.5287e-225.675e-2260.2465WOA5.4848e-069.8198e-041.0474e-0370.0146PSO3.0552e-044.1497e-018.0658e-0180.1467DE00010.0560F12IMPA-10.4029-10.40291.3194e-1610.2139MSIMPA-10.4029-10.40295.9774e-1330.3531MPA-10.4029-10.40292.3540e-1140.1048ESSA-10.4029-10.40297.3759e-1620.1252IGWO-10.4029-10.40292.3253e-0850.3106WOA-10.4029-7.43263.0665e+0080.0338PSO-10.4028-10.31542.6727e-0160.1729DE-10.4029-10.00451.5274e+0070.0771F13IMPA-3.8628-3.86282.6662e-1540.1102MSIMPA-3.8628-3.86282.6684e-1550.2435MPA-3.8628-3.86282.3715e-1510.0765ESSA-3.8628-3.86282.4184e-1520.1033IGWO-3.8628-3.86282.4795e-1530.2701WOA-3.8628-3.85820.006503270.0233PSO-3.862-3.84654.6015e-0280.1550DE-3.8628-3.86282.7101e-1560.0651F14IMPA00010.1241MSIMPA03.5595e-301.3958e-2930.1955MPA3.4592e-283.0451e-221.5159e-2150.0588ESSA02.8995e-241.1012e-2340.0891IGWO7.4237e-217.1524e-171.4177e-1660.2510WOA1.4198e-131.0161e-012.6348e-0180.0161PSO1.2304e-041.3578e-039.7707e-0470.1463DE00010.0575F15IMPA1.3498e-311.3498e-316.6809e-4710.1196MSIMPA1.3498e-311.3498e-316.6809e-4710.1903MPA7.2938e-276.392e-231.3055e-2260.0567ESSA1.3498e-311.3498e-316.6809e-4710.0890IGWO7.8004e-306.2979e-262.7686e-2550.2485WOA2.2215e-098.1541e-061.6292e-0570.0154PSO1.2427e-071.6901e-032.9564e-0380.1472DE1.3498e-311.3498e-316.6809e-4710.0561Table 4
Comparison of benchmark function optimization results (100-dimension).
FunctionAlgorithmBestAveStdRankRunningtime (s)F1IMPA00010.4799MSIMPA00011.7578MPA2.0936e-202.0131e-191.8807e-1950.2724ESSA1.3162e-2157.4533e-314.0189e-3040.1457IGWO3.4070e-132.7008e-122.0152e-1260.5957WOA4.056e-821.6867e-718.8984e-7130.0481PSO2.5648e-063.4789e-011.4948e+0070.2561DE6.8945e+021.2783e+035.3839e+0280.1209F2IMPA02.5048e-297020.5490MSIMPA00011.9195MPA1.2682e-131.418e-111.0044e-1150.2828ESSA4.6552e-622.6374e-141.4083e-1340.1475IGWO7.0279e-091.5616e-085.995e-0960.6088WOA5.5313e-612.3061e-499.5159e-4930.0489PSO4.8680e-033.9462e-014.8301e-0170.2833DE1.2305e+011.9098e+013.7135e+0080.1235F3IMPA9.2984e+019.3322e+011.7465e-0130.6006MSIMPA9.4100e+019.4491e+011.7690e-0142.4819MPA9.5854e+019.7101e+017.8353e-0150.3002ESSA8.5565e-093.6629e-045.9162e-0410.1675IGWO9.4975e+019.7472e+011.1160e+0060.6292WOA9.7476e+019.8149e+012.5179e-0170.0606PSO6.1575e-047.9545e+014.3079e+0220.2824DE3.7089e+051.5227e+061.0087e+0680.1382F4IMPA00010.5204MSIMPA7.3420e-025.7584e-013.4773e-0131.6146MPA3.0751e+003.9660e+006.4026e-0140.2583ESSA8.1256e-116.5052e-071.3268e-0620.1420IGWO5.4614e+007.8873e+001.0375e+0060.5740WOA1.6505e+004.1532e+009.5339e-0150.0437PSO3.0226e-073.4222e+011.3481e+0270.2539DE5.8683e+021.0714e+033.8677e+0280.1172F5IMPA00010.4683MSIMPA00011.4869MPA6.3693e-082.25e-071.1712e-0740.2249ESSA4.8795e-814.7458e-152.3905e-1430.1178IGWO1.1614e+004.6770e+002.8552e+0060.5037WOA2.3009e+018.6982e+011.6871e+0170.0388PSO3.4590e-049.5632e-031.0224e-0250.2272DE8.3771e+018.7385e+011.4535e+0080.1050F6IMPA8.8818e-168.8818e-16010.5018MSIMPA8.8818e-168.8818e-16011.8198MPA1.2297e-115.5801e-114.2214e-1150.2448ESSA8.8818e-161.4803e-151.3467e-1530.1498IGWO6.1236e-081.7033e-076.0547e-0860.5292WOA8.8818e-164.6777e-152.6279e-1540.0486PSO3.1384e-045.8055e-021.8421e-0170.2642DE4.7890e+006.2446e+008.0832e-0180.1407F7IMPA00010.5311MSIMPA00012.6065MPA00010.2776ESSA00010.1679IGWO2.7289e-132.9180e-036.3530e-0360.5641WOA00010.0631PSO7.2014e-071.4439e+014.1840e+0180.2778DE5.0666e+001.3019e+016.3009e+0070.1553F8IMPA3.0053e-061.8797e-051.4073e-0410.9748MSIMPA2.2357e-067.5344e-056.8927e-0522.0507MPA3.8496e-041.8048e-035.9696e-0440.5977ESSA3.3627e-057.2115e-047.7153e-0430.3942IGWO0.00542031.5155e-024.5389e-0360.8735WOA2.6318e-055.3274e-035.2580e-0350.2172PSO4.9434e-036.7585e-024.1483e-0270.4266DE1.4769e+003.0316e+001.0210e+0080.2811F9IMPA-40708.7689-33742.4422302.829130.5864MSIMPA-35220.1185-33533.4499895.514741.9733MPA-26300.7559-24710.3869907.327550.2990ESSA-28921.6901-23487.29372711.124460.1623IGWO-26070.9887-16182.55656433.380070.5971WOA-41895.7857-34167.74925693.474420.0678PSO-41898.2886-40509.56263560.758910.2540DE-10149.0836-60561.65082047.619980.1443F10IMPA4.7116e-334.7116e-331.3918e-4811.8984MSIMPA2.1453e-036.4522e-033.3446e-0343.7580MPA3.1036e-024.6095e-028.6771e-0351.0844ESSA1.0014e-154.9800e-091.0483e-0820.7149IGWO1.3000e-012.2248e-015.3100e-0271.4428WOA1.9735e-025.3398e-022.6769e-0260.4637PSO6.4788e-092.0935e-056.6560e-0530.6519DE1.4390e+044.3132e+055.7432e+0580.5338


## 5.2. Convergence Analysis

In this paper, the average convergence curves of the algorithms are selected to evaluate the convergence ability of the optimization algorithms. Figure7 shows each algorithm’s 30-dimensional average convergence curves on the 15 benchmark test functions. The IMPA has excellent convergence speed and convergence accuracy on the remaining 14 test functions except for the F3. In F1, F2, F5, and F7, both IMPA and MSIMPA converge to the theoretical optimum quickly. In F12 and F13, the optimization precision and convergence speed of the IMPA improved a little over the other seven algorithms. It is worth mentioning that the DE can quickly converge to the global optimal value in F11, F14 and F15, which shows excellent optimization ability. In addition, it can be seen from Tables 3 and 4 that the running time of the proposed IMPA algorithm on most test functions is about twice that of MPA, but the optimization accuracy of IMPA has been greatly improved. Therefore, it is worth sacrificing a little time to obtain a higher quality of the solution. IMPA can solve the issues of low optimization accuracy and low convergence rate of MPA and shows significant improvements in the convergence speed and solution accuracy. These improvements are due to the implementation of the population evolution strategy and the dynamic learning mechanism, which effectively enhance the algorithm’s ability for the global searching while enriching the diversity of the population.Figure 7
Comparison of convergence curves.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)


## 5.3. The Wilcoxon Rank Sum Test

Only relying on the comparison of mean and standard deviation values, the superiority of the IMPA cannot be approved entirely [57]. Therefore, it needs to be tested at the statistical level. In this paper, to distinguish the significant differences in the optimization performance of the IMPA over the other seven algorithms, Wilcoxon statistical tests are performed. Test results of 30-dimension and fixed-dimension on the benchmark functions at the significant level of p=5% are obtained. The test results are shown in Table 5. When p<5%, the zero hypothesis considers a rejection, indicating that there is a significant difference between both algorithms. When p>5%, it indicates that the difference between both algorithms is not substantial. Since the best algorithm on a benchmark function cannot be compared with itself, the best algorithm on each benchmark function is marked as N/A, indicating that the performance is comparable and optimal between the two. “S” represents the saliency judgment result, and the symbols “+”, “-”, and “=,” respectively, indicate that the performance of IMPA is better than, worse than, and equivalent to the comparison algorithm. As shown in Table 5, most of the p values are less than 5%, indicating that the IMPA algorithm is significantly different in most of the test functions among the other seven algorithms. For F1, F2, F5-F7, F11, and F15, the performance of IMPA and MSIMPA are comparable, and all can obtain the best value. In general, the IMPA algorithm shows significant advantages over the other seven algorithms.Table 5
Thep value of the Wilcoxon rank sum test on 15 benchmark functions (30-dimension and fixed-dimension).
FunctionMSIMPAMPAESSAIGWOWOAPSODEp value Sp value Sp value Sp value Sp value Sp value Sp value SF1N/A=1.212e-12+1.212e-12+1.212e-12+1.212e-12+1.212e-12+1.212e-12+F2N/A=1.212e-12+4.574e-12+1.212e-12+1.212e-12+1.212e-12+1.212e-12+F33.020e-11+3.020e-11+3.020e-11+3.020e-11+3.020e-11+5.573e-10+3.020e-11+F41.212e-12+1.212e-12+1.212e-12+1.212e-12+1.212e-12+1.212e-12+1.212e-12+F5N/A=1.212e-12+4.574e-12+1.212e-12+1.212e-12+1.212e-12+1.212e-12+F6N/A=1.212e-12+4.191e-02+1.031e-12+1.163e-09+1.212e-12+1.212e-12+F7N/A=N/A=N/A=6.617e-04+3.337e-01-1.212e-12+1.212e-12+F84.515e-02+3.020e-11+3.081e-08+3.020e-11+7.739e-06+3.020e-11+3.020e-11+F99.117e-01-3.020e-11+1.777e-10+6.722e-10+9.941e-01-3.256e-07+3.020e-11+F101.212e-12+1.212e-12+1.212e-12+1.212e-12+1.212e-12+1.212e-12+1.212e-12+F11N/A=4.574e-12+3.337e-01-1.212e-12+1.212e-12+1.212e-12+N/A=F122.758e-11+1.449e-11+8.435e-06+1.449e-11+1.449e-11+1.136e-11+7.338e-03+F133.128e-01-9.452e-09+3.568e-07+2.119e-05+1.720e-12+1.720e-12+3.337e-03+F142.19e-06+1.212e-12+1.657e-11+1.212e-12+1.212e-12+1.212e-12+N/A=F15N/A=1.212e-12+N/A1.212e-12+1.212e-12+1.212e-12+N/A=+/=/-6/7/214/1/013/1/115/0/013/0/215/0/012/3/0


## 6. Application of 3D Surface WSN Coverage Based on IMPA

In this section, simulation of the performed experiments is conducted under two kinds of 3D terrains named undulating and rugged terrains. The flow chart of the proposed IMPA algorithm applied to 3D surface WSN coverage optimization is shown in Figure8. The three-dimensional maps of both terrains are shown in Figure 9. In this research, to simplify the problem model, both terrains are simulated using the surface equation of the function. And the simulation results for IMPA, MSIMPA, MPA, ESSA, and WOA are compared.Figure 8
The flow chart of the IMPA applied in WSN.Figure 9
3D maps of two terrains.
(a)
Undulating terrain(b)
Rugged terrain


### 6.1. Undulating Terrain Simulation Experiment

The surface equation of the simulated undulating terrain isz=xe−x2−y2. The length and width of the mapped surface into the 2D plane are set to 4 m. The sensing radius of the sensor nodes is0.5m. The population size and iterations set to 30 and 1000, respectively. The experiments are independently run 20 times for each group to obtain the average result. The surface is divided into 400 sections using the grid division method, and deployed with 30,40,⋯,80 nodes. The simulation results are as shown in Table 6.Table 6
Comparison of average coverage of undulating terrain.
AlgorithmAverage coverage (%)30 nodes40 nodes50 nodes60 nodes70 nodes80 nodesWOA61.9469.4677.5879.0184.188.98ESSA64.0771.7679.4881.3984.8487.13MPA70.1580.7986.6190.2894.3794.8MSIMPA72.2181.3787.6791.6394.7795.33IMPA74.1983.290.0593.696.5996.81Table6 revealed that the average coverage increases with the increase of sensor nodes. Under the same conditions, the network coverage of sensor nodes optimally deployed by the IMPA algorithm proposed in this paper is the largest. The coverage optimization ability of WOA and ESSA is poor, and there is a big gap compared with the other three algorithms. When the number of nodes is 50, the average coverage of the IMPA algorithm is 3.44% greater than that of the MPA and 12.47% higher than the WOA. When the number of nodes is 80, the average coverage of the IMPA reaches to 96.81%. The average coverage of the IMPA does not increase significantly after 70 nodes. It is mainly due to the redundancy of nodes in the network, weakening the algorithm’s performance.Figures10 and 11 show the optimized deployment results of the MPA and IMPA at 50 nodes, respectively. As shown in Figure 10(b), in order to understand the distribution location and coverage areas of the nodes, we map the 3D coverage results into a 2D plane. The red dots indicate the deployment locations of the sensor nodes and the blue areas indicate the coverage areas. Figure 12 shows the coverage iteration curves of the five algorithms.Figure 10
MPA optimized deployment on undulating terrain (50 nodes).
(a)
MPA optimized deployment(b)
MPA optimized deployment of 2D viewFigure 11
IMPA optimized deployment on undulating terrain (50 nodes).
(a)
IMPA optimized deployment(b)
IMPA optimized deployment of 2D viewFigure 12
Convergence comparison in undulating terrain (50 nodes).Figures10(b) and 11(b)show that the sensor nodes deployed by the IMPA are more uniform, and the coverage area is larger. Figure 12 shows that the convergence rate of the IMPA is more rapid compared to the MPA. At about the 200-th generation, the coverage rate of the IMPA reaches to 88%, while the MPA algorithm is only about 72%, a difference of about 16% between the two. This shows that the proposed improved strategy can speed up the optimization speed of the algorithm.


### 6.2. Rugged Terrain Simulation Experiment

Simulations using more complex rugged terrain are conducted to further validate the deployment capabilities of the IMPA algorithm. The simulated surface equation of the rugged terrain isz=cosxsiny. The length and width of the mapped surface into the 2D plane are adjusted to 10 m. The sensing radius of the sensor nodes is 1 m. The population size and iterations are set to 50 and 300, respectively. The simulations are independently run 20 times for each group to obtain the average value. In this simulation, the number of sensor nodes is adjusted to 100 because the rugged terrain is complex, and the 3D perception blind area is larger. The surface is also divided into 400 sections. The results of the optimized deployment of the IMPA are shown in Figure 13.Figure 13
IMPA optimized deployment on rugged terrain.
(a)
IMPA optimized deployment(b)
IMPA optimized deployment of 2D viewFigure13 shows that the optimized nodes by the IMPA are uniformly deployed on the surface. The coverage area significantly becomes greater while comparing to the initial deployment. However, deploying nodes in the convex region of the 3D surface due to the existence of the 3D perception blind area would be difficult. It results in a relatively concentrated coverage hole which is shown in the blank area of Figure 13(b). Figure 14 illustrates the coverage convergence curves of both algorithms in rugged terrain. It shows that the convergence rate of the IMPA is very high within the early stages. Around the 200-th generation, the coverage quickly reaches to 87.5% where its final coverage rate is optimized to 92.57%. The results indicate that the final coverage for the IMPA is 6.82% greater than the MPA (85.75%). It shows that the IMPA algorithm integrating the three proposed improved strategies of this research could maintain good deployment performance within the complex rugged terrain. Moreover, the IMPA algorithm shows good terrain adaptability.Figure 14
Convergence comparison in rugged terrain.In summary, by comparing both simulation results for the different terrains, the IMPA could effectively improve the coverage while solving the coverage problem for the complex 3D surface WSN. The nodes could be deployed more uniformly and have more robust adaptability to the different terrains. Finally, it verifies the improvement and effectiveness of the proposed strategies.


### 6.3. Network Life Cycle Comparison

The control of energy consumption is very important to the network life of WSN and is directly related to the life cycle of the entire network. Under the same conditions, the smaller the energy consumption of WSN, the longer the life cycle of the network. The LEACH protocol is a low-energy routing algorithm commonly used in WSN. The basic idea is to distribute the energy load of the entire network evenly to each sensor node by randomly and cyclically selecting the cluster head node, thereby reducing the energy consumption of the network and improving the overall survival time of the network. In this paper, the LEACH protocol is used for data transmission, and the specific energy consumption model is shown in the literature [58]. Reference [59] uses the death time of the last node to represent the life cycle of the network. Therefore, this paper also uses the death time of the last node to test the network life cycle.The comparison of different algorithms in undulating terrain is shown in Figure15. It can be seen from the figure that the network lifetime increases as the number of deployed sensor nodes increases and the proposed IMPA algorithms all obtain the highest network life cycle. Interestingly, the MSIMPA is able to achieve a higher coverage than the MPA, but his network life cycle is shorter than MPA. The comparison of the network life cycle of rugged terrain is shown in Figure 16. It can be seen that the network life cycle optimized by IMPA is the highest, reaching 4333 rounds, followed by the WOA algorithm, but the network coverage rate optimized by WOA is not ideal. Overall, the WSN optimized by the proposed IMPA algorithm can achieve higher network coverage and longer life cycle.Figure 15
Comparison of the life cycle of the five algorithms in undulating terrain.Figure 16
Comparison of the life cycle of the five algorithms in rugged terrain.


## 6.1. Undulating Terrain Simulation Experiment

The surface equation of the simulated undulating terrain isz=xe−x2−y2. The length and width of the mapped surface into the 2D plane are set to 4 m. The sensing radius of the sensor nodes is0.5m. The population size and iterations set to 30 and 1000, respectively. The experiments are independently run 20 times for each group to obtain the average result. The surface is divided into 400 sections using the grid division method, and deployed with 30,40,⋯,80 nodes. The simulation results are as shown in Table 6.Table 6
Comparison of average coverage of undulating terrain.
AlgorithmAverage coverage (%)30 nodes40 nodes50 nodes60 nodes70 nodes80 nodesWOA61.9469.4677.5879.0184.188.98ESSA64.0771.7679.4881.3984.8487.13MPA70.1580.7986.6190.2894.3794.8MSIMPA72.2181.3787.6791.6394.7795.33IMPA74.1983.290.0593.696.5996.81Table6 revealed that the average coverage increases with the increase of sensor nodes. Under the same conditions, the network coverage of sensor nodes optimally deployed by the IMPA algorithm proposed in this paper is the largest. The coverage optimization ability of WOA and ESSA is poor, and there is a big gap compared with the other three algorithms. When the number of nodes is 50, the average coverage of the IMPA algorithm is 3.44% greater than that of the MPA and 12.47% higher than the WOA. When the number of nodes is 80, the average coverage of the IMPA reaches to 96.81%. The average coverage of the IMPA does not increase significantly after 70 nodes. It is mainly due to the redundancy of nodes in the network, weakening the algorithm’s performance.Figures10 and 11 show the optimized deployment results of the MPA and IMPA at 50 nodes, respectively. As shown in Figure 10(b), in order to understand the distribution location and coverage areas of the nodes, we map the 3D coverage results into a 2D plane. The red dots indicate the deployment locations of the sensor nodes and the blue areas indicate the coverage areas. Figure 12 shows the coverage iteration curves of the five algorithms.Figure 10
MPA optimized deployment on undulating terrain (50 nodes).
(a)
MPA optimized deployment(b)
MPA optimized deployment of 2D viewFigure 11
IMPA optimized deployment on undulating terrain (50 nodes).
(a)
IMPA optimized deployment(b)
IMPA optimized deployment of 2D viewFigure 12
Convergence comparison in undulating terrain (50 nodes).Figures10(b) and 11(b)show that the sensor nodes deployed by the IMPA are more uniform, and the coverage area is larger. Figure 12 shows that the convergence rate of the IMPA is more rapid compared to the MPA. At about the 200-th generation, the coverage rate of the IMPA reaches to 88%, while the MPA algorithm is only about 72%, a difference of about 16% between the two. This shows that the proposed improved strategy can speed up the optimization speed of the algorithm.


## 6.2. Rugged Terrain Simulation Experiment

Simulations using more complex rugged terrain are conducted to further validate the deployment capabilities of the IMPA algorithm. The simulated surface equation of the rugged terrain isz=cosxsiny. The length and width of the mapped surface into the 2D plane are adjusted to 10 m. The sensing radius of the sensor nodes is 1 m. The population size and iterations are set to 50 and 300, respectively. The simulations are independently run 20 times for each group to obtain the average value. In this simulation, the number of sensor nodes is adjusted to 100 because the rugged terrain is complex, and the 3D perception blind area is larger. The surface is also divided into 400 sections. The results of the optimized deployment of the IMPA are shown in Figure 13.Figure 13
IMPA optimized deployment on rugged terrain.
(a)
IMPA optimized deployment(b)
IMPA optimized deployment of 2D viewFigure13 shows that the optimized nodes by the IMPA are uniformly deployed on the surface. The coverage area significantly becomes greater while comparing to the initial deployment. However, deploying nodes in the convex region of the 3D surface due to the existence of the 3D perception blind area would be difficult. It results in a relatively concentrated coverage hole which is shown in the blank area of Figure 13(b). Figure 14 illustrates the coverage convergence curves of both algorithms in rugged terrain. It shows that the convergence rate of the IMPA is very high within the early stages. Around the 200-th generation, the coverage quickly reaches to 87.5% where its final coverage rate is optimized to 92.57%. The results indicate that the final coverage for the IMPA is 6.82% greater than the MPA (85.75%). It shows that the IMPA algorithm integrating the three proposed improved strategies of this research could maintain good deployment performance within the complex rugged terrain. Moreover, the IMPA algorithm shows good terrain adaptability.Figure 14
Convergence comparison in rugged terrain.In summary, by comparing both simulation results for the different terrains, the IMPA could effectively improve the coverage while solving the coverage problem for the complex 3D surface WSN. The nodes could be deployed more uniformly and have more robust adaptability to the different terrains. Finally, it verifies the improvement and effectiveness of the proposed strategies.


## 6.3. Network Life Cycle Comparison

The control of energy consumption is very important to the network life of WSN and is directly related to the life cycle of the entire network. Under the same conditions, the smaller the energy consumption of WSN, the longer the life cycle of the network. The LEACH protocol is a low-energy routing algorithm commonly used in WSN. The basic idea is to distribute the energy load of the entire network evenly to each sensor node by randomly and cyclically selecting the cluster head node, thereby reducing the energy consumption of the network and improving the overall survival time of the network. In this paper, the LEACH protocol is used for data transmission, and the specific energy consumption model is shown in the literature [58]. Reference [59] uses the death time of the last node to represent the life cycle of the network. Therefore, this paper also uses the death time of the last node to test the network life cycle.The comparison of different algorithms in undulating terrain is shown in Figure15. It can be seen from the figure that the network lifetime increases as the number of deployed sensor nodes increases and the proposed IMPA algorithms all obtain the highest network life cycle. Interestingly, the MSIMPA is able to achieve a higher coverage than the MPA, but his network life cycle is shorter than MPA. The comparison of the network life cycle of rugged terrain is shown in Figure 16. It can be seen that the network life cycle optimized by IMPA is the highest, reaching 4333 rounds, followed by the WOA algorithm, but the network coverage rate optimized by WOA is not ideal. Overall, the WSN optimized by the proposed IMPA algorithm can achieve higher network coverage and longer life cycle.Figure 15
Comparison of the life cycle of the five algorithms in undulating terrain.Figure 16
Comparison of the life cycle of the five algorithms in rugged terrain.


## 7. Conclusion and Future Work

This study proposed a multistrategy integrated improved marine predator algorithm (IMPA) for coverage optimization of WSN on the 3D surfaces. The fusion of the random opposition-based learning strategy and the differential evolution operator was employed to expand the search range and improve the diversity of the population to enhance the global search ability. Furthermore, the improved SFLA strategy, as a local search algorithm, was employed to replace the original FADs effect of MPA, reducing the blind search of individuals and improving the algorithm’s local exploitation ability. Moreover, a dynamic learning mechanism based on quasireflected opposition-based learning was proposed to improve the algorithm’s convergence speed and optimization search accuracy. The IMPA algorithm is compared with the other seven algorithms in 15 benchmark test functions. The results show that IMPA has significant advantages in solution accuracy and convergence speed in both 30 and 100 dimensions. Finally, IMPA is applied to the node deployment of WSN on two 3D surfaces. The simulation results show that the nodes optimized by IMPA are more uniform, have higher coverage, longer network life cycle, and have faster convergence speed. It also shows that the IMPA has strong adaptability to the two 3D simulated terrains and has certain practicability and effectiveness.Because of the complexity of the 3D surfaces and the existence of 3D perception blind area, it is difficult to ensure the connectivity of the deployed networks. In the future, we will study appropriate methods to ensure the connectivity of WSN deployed on 3D surfaces.



---

*Source: 9593103-2022-11-17.xml*