An Improved Strength Pareto Evolutionary Algorithm 2 with Adaptive Crossover Operator for Bi-Objective Distributed Unmanned Aerial Vehicle Delivery

Abstract

With the development of the e-commerce industry, using UAVs (unmanned aerial vehicles) to deliver goods has become more popular in transportation systems. This delivery method can reduce labor costs and improve the distribution efficiency, and UAVs can reach places that are difficult for humans to reach. Because some goods are perishable, the quality of the delivery will have an impact on the customer satisfaction. At the same time, the delivery time should also meet the needs of customers as much as possible. Therefore, this paper takes the distribution distance and customer satisfaction as the objective functions, establishes a bi-objective dynamic programming model, and proposes an improved SPEA2 (strength Pareto evolutionary algorithm 2). The improved algorithm introduces the local search strategy, on the basis of the original algorithm. It conducts a local search for the better non-dominated solutions obtained in each iteration. The new dominated solutions and non-dominated solutions are determined, and the crossover operator is improved, so that the local search ability is improved, on the basis of ensuring its global search ability. The numerical experiment results show that the improved algorithm achieves an excellent performance in three aspects: the Pareto front, generation distance, and spacing, and would have a high application value in UAV cargo delivery and other MOPs (multi-objective optimization problems). The average spacing value of the improved algorithm is more than 20% smaller than SPEA2 + SDE (strength Pareto evolution algorithm 2–shift-based density estimation), which is the second-best algorithm. In the comparison of the average generation distance value, this number reaches 30%.

1. Introduction

With the development of the e-commerce industry, great changes have taken place in the field of logistics and delivery. In order to reduce labor costs and improve the delivery efficiency, more and more logistics companies are beginning to use UAVs (Unmanned Aerial Vehicles) to deliver goods, although the UAV logistics technology is not mature. At the same time, compared with the traditional logistics industry, it faces a series of constraints, such as air control. However, the advantages of UAV logistics systems are very significant, such as flexibility, ease of use, and operability. These qualities ensure they have a place in the future development of the logistics industry [1].

In the existing UAV delivery research, the main focus is on the shortest total path of UAV flight [2]. Although this is an important factor affecting revenue, the impact of time on UAV delivery has not been given enough attention. Generally, customers have certain requirements for the delivery time of goods. If we can deliver the goods within the time required by customers, we can obtain a higher level of satisfaction. On the other hand, more and more people are beginning to buy perishable goods, such as fruits and vegetables, through online shopping [3]. The quality of such goods at the time of delivery will greatly affect the customer satisfaction. These goods can easily decay during transportation. In some countries, the loss rate during transportation is as high as 30%, which is even more serious in underdeveloped regions [4]. This will also have a high negative impact on the local economy.

In addition, with improvements in shopping platform functions, consumers can evaluate the products they buy on a platform. If the purchased items do not meet expectations in terms of the freshness, quality, integrity, etc., the customer will give poor comments to the seller. This will greatly reduce the desire of other consumers to shop at that business. Through observing the evaluation areas of several e-commerce companies, Zhang et al. [5] found that customers are usually most concerned with the freshness, quality, and delivery efficiency of goods.

To summarize, the object of this study is to identify not only the shortest delivery distance of a UAV, but also how customer satisfaction is affected by delayed delivery, and the deterioration of perishable products. In other words, this paper studies a two-objective optimization problem. Therefore, a MOEA (multi-objective evolutionary algorithm) will be used to solve the problem. In short, this paper has the following differences from the existing research of this type. First of all, when considering customer satisfaction, the existing research is usually based on the delivery time. However, in this paper, in addition to considering this factor, we acknowledge that the customer satisfaction is also affected by the quality of the goods delivered. Secondly, we propose an improved SPEA2 (strength Pareto evolutionary algorithm 2) to solve the problem. Through a local search and improved crossover operator, the improved algorithm avoids the defects of easy-to-fall-into local optimization and a slow convergence speed. The relationship between customer satisfaction and distance cost is expressed through the Pareto frontier, which is selected and judged by the sellers themselves. Finally, the concept of convergence distance is proposed, to judge the convergence speed of the algorithm.

Compared with existing similar studies, this paper makes the following contributions:

• This paper establishes a bi-objective model to solve the UAV delivery problem. The previous way of thinking, which only considered the total cost of transportation, is changed, to face this type of problem from a multi-objective perspective, and the effect of the freshness of perishable products is taken into account when calculating the customer satisfaction.
• An improved SPEA2 is proposed in this paper. Considering that SPEA2 itself performs well in global searching, this paper conducts local searching for the points with higher adaptation after each iteration, meaning that it takes both the global and local searching ability into account. In addition, the crossover operator is also improved, to increase the convergence speed of the algorithm, and avoid falling into local optimality.
• The concept of convergence distance is proposed in this paper. The convergence speed of the algorithm is judged by calculating the distance between the experimental results and the theoretical optimal results.

The other chapters of this article are as follows. In the second chapter, the relevant background of the MOEA, memetic algorithms, and UAV delivery is summarized, which provides a basis for the algorithm improvement and model construction. The third chapter constructs and analyzes the UAV delivery model and its constraint function. The fourth chapter improves the SPEA2. Based on the original algorithm, a local search strategy and new crossover operator are introduced, to improve the algorithm. The fifth chapter compares the improved algorithm with other MOEAs, under nine benchmarks [6], and verifies the superiority, stability, and convergence of the improved algorithm. Then, the parameters are determined through experiments, and the improved algorithm is used to solve the model. The improved algorithm can take into account both the distance cost and time cost, and is better than other MOEAs. The last chapter summarizes the work of the paper, points out the advantages and disadvantages, and proposes the focus of future work.

2. Related Work

2.1. Background of the MOEA

So far, there have been more than 20 methods proposed to solve MOPs (multi-objective optimization problems) in operational research. Most of them follow a fixed pattern. They use specific tech MOPs in multiple single-objective optimization problems. In this way, an optimal solution can be easily obtained. However, the solution of MOPs is not a single solution, but a solution set, and all the points in the set must meet the conditions [7]. In addition, results obtained by same method may be different. This is because the solution is closely related to “decision-makers”. Different decision-makers focus on different goals.

The evolutionary algorithm is a common method of solving optimization problems. Its advantage is that it can deal with a very large search space, which is in line with the characteristics of MOPs. Therefore, many scholars consider using an evolutionary algorithm to optimize multiple objectives. The MOEA is derived from this. Compared with traditional methods, it can obtain multiple optimal solutions in one run [8].

The research history of the MOEA is not long. Due to the different main characteristics of the algorithm [9], its evolution can be considered to comprise two stages [10]:

(1)

First-generation MOEA

In 1985, Schaffer [11] improved the single-objective genetic algorithm, and proposed Vega, which was the first real MOEA. Goldberg first introduced Pareto into the genetic algorithm in 1989, which had a significant and far-reaching impact on the development of multi-objective optimization. Subsequently, MOEAs sprung up: in 1993, the MOGA (multi-objective genetic algorithm) was proposed by Fonseca and Fleming [12]; in 1994, Srinivas et al. [13] proposed the NSGA (non-dominated sorting genetic algorithm) and improved the selection operator; in the same year, Horn et al. [14] proposed the NPGA (niched Pareto genetic algorithm). In brief, the first generation MOEA is characterized by the use of the selection strategy based on the Pareto level, and the diversity maintenance strategy based on the fitness sharing mechanism [15].

(2)

Second-generation MOEA

Since entering the 21st century, with the rapid development of the related technologies, many classical algorithms have been proposed. Zitzler et al. proposed the SPEA (strength Pareto evolution algorithm) [16] and the SPEA2 [17]. The SPEA2 is one of the representative classical algorithms of the second generation. The algorithm adds the external archive set, and maintains the scale of the external archive set through the clustering and pruning methods. In addition, the concept of Pareto dominance is also introduced into the assignment calculation of the population individual fitness. Deb et al. proposed the classical algorithm NSGA-II (non-dominated sorting genetic algorithm-II) [18] in the second generation MOEA. The algorithm uses better accounting strategies, according to various Pareto-dominated frontiers, which reduces the running time of the algorithm. The other major MOEAs also include the PAES (Pareto archived evolution strategy) [19] and the NPGA2 (niched Pareto genetic algorithm 2) [20], the PESA (Pareto envelope-based selection algorithm) [21] and the PESA-II (Pareto envelope-based selection algorithm-II) [22]. The remarkable feature of the second-generation MOEA is the elite individual preservation mechanism.

(3)

Classification of MOEA

According to different environmental selection strategies, there are four types of current MOEAs [23,24,25]: (a) dominance-based; (b) indicator-based; (c) decomposition-based; and (d) others.

The SPEA2 is a force-based approach among many MOEAs. After each iteration, a series of dominant and non-dominant solutions will be obtained, and then the applicability of the solution will be judged according to the quantity of its dominating non-dominant solutions. In addition, an external file will be set to store the non-dominant solution. In each iteration, these solutions are reevaluated, and the external archive is updated. Whether the solution will be used depends on its strength after each update [26]. Finally, the solutions in the external file, and the non-dominated solutions in the population make up the optimal solution set.

2.2. Background of Memetic Algorithms

Memetic algorithms (MAs) were first introduced by Moscato et al. [27,28]. Memetics represents the unity of cultural evolution, which can spread, and show local enhancement. Therefore, these algorithms are mostly used to improve the local search ability of genetic algorithms [29].

Initially, MAs focused on genetic algorithms. Later, they were extended to improve the local search capability of evolutionary strategies. They can also be used to improve swarm-based and other population-based metaheuristics [30]. In addition, several types of local search strategies have been adopted, such as simulated annealing [31], tabu search [32], and pattern search [33], to improve the solution refinement capabilities of the swarm-based metaheuristics.

One of the reasons for the increasing popularity of MAs is their performance on high-dimensional and non-differentiable search spaces, and their enhanced exploitation capabilities [34]. Therefore, many papers have proposed, improved, or used MAs for discrete, continuous, constrained, multi-objective, and other classes of optimization problems [34].

2.3. Background of Goods Delivery

As people’s consumption of perishable products increases, the research on the delivery of perishable goods is also increasing steadily. First of all, we will review the literature on single-objective problems in this field.

Tarantilis and Kiranoudis [35] solved the problem of fresh milk distribution with heterogeneous fleets in 2001. They proposed an algorithm based on the threshold acceptance, but the goal was to minimize the distance. Rabbani et al. [36] considered the loss of freshness in their research, and included it in the total cost. However, it was still a single-target problem, which could be solved by a genetic algorithm. Alvarez and Cordeau et al. [37] regarded the corrosion of goods as the inventory cost, which was one of the factors affecting the total profit, and then used the branch-cut algorithm to maximize the total profit. Liang et al. [38] wanted to use multiple distribution tools to reduce the transportation distance, so they considered the delivery of different perishable goods. When considering the decay of goods, they also considered the impact of the temperature and humidity on the goods. However, they did not study the impact of goods decay on customer satisfaction. Compared with the usual targets in the past, Meneghetti and Ceschia [39], when studying such problems, instead considered the fuel consumption of the transportation vehicles. Because the goods needed to be refrigerated during the distribution process, the fuel consumption would be increased. Based on this point, they conducted research to increase profits.

As mentioned above, these are all goods distribution and delivery problems with only a single goal, and the goal is usually an economic one; for example, the shortest distance, the largest total revenue, and the lowest fuel consumption. Although researchers have taken into account the impact of goods decay on earnings, this is only one factor. When it is considered as a single goal, a multi-objective optimization framework can provide new ideas and better solutions. Therefore, we will review the multi-objective research in this area.

Wang et al. [40] took the total cost and the number of vehicles as the two objectives of their study, and established a multi-database model. In addition to solving the model, they also added the service time, as a third goal, to the model, and built a three-goal model. When setting targets, Tirkolaee and Aydin [41] focused on gas emissions and job opportunities, in addition to costs. In their research, whether perishable goods deteriorated was an uncertain parameter, and they used the fuzzy method to solve it. In addition, in the delivery of perishables, the customer satisfaction will also have a great impact on revenue. This concept was first proposed by Cardozo [42]. He believed that customers would buy the products they were satisfied with again, and consider other products from the business. At the same time, customers would also recommend the products they were satisfied with to other customers, affecting further purchases. Rahami and Armand [43] took customer satisfaction into account when establishing a dual-objective model for perishable products. In their model, a delayed delivery would reduce the customer satisfaction, and the customer satisfaction could be expressed in this way. In the research of Wang et al. [44], time and freshness were two factors that affected the customer satisfaction. The researchers set priorities for customers’ needs, and solved this problem by modifying the priority function.

As far as we know, few studies have taken the transportation cost and customer satisfaction as the two objectives for goods delivery. In previous studies, the customer satisfaction was usually determined by the customer waiting time. However, when perishable goods are delivered, their quality has a great impact on the customer satisfaction [45]. This study establishes a model of the impact of the perishable product quality and the time window on the customer satisfaction. The details will be described in the next chapter.

3. Model Building

3.1. Problem Statement and Assumptions

In the problems studied in this paper, the supplier must achieve two goals; namely, reducing the distance cost of delivering the products, and improving the customer satisfaction. The customer satisfaction is affected by the quality of the product, and by the waiting time of the customer at the time of delivery. Therefore, each customer’s needs, and the expected delivery time, must be known. On this premise, we must study how to reduce costs and improve customer satisfaction at the same time. To solve this problem, the following assumptions must be made first [46]:

(1)

The number of UAVs in the distribution center is enough to complete all tasks.

(2)

Each customer can only submit one order, and can only be served once by a UAV.

(3)

The same order cannot be completed by multiple UAVs; that is, the cargo demand of any order cannot exceed the maximum load capacity of a single UAV.

(4)

All the products are fresh before delivery, and the quality decreases during transportation.

The goods in this paper are perishable products. Therefore, the relationship between the degree of product decay and the time must be known. These products can be divided into two categories: products with a fixed period of validity, such as beverages and biscuits, and products with an indefinite period of validity, such as fruits and vegetables. For the first type of product, the deterioration degree is easy to calculate, but the deterioration degree of the other type of product is affected by many factors. Labuza [47] has discovered that the deterioration of the second type of product conforms to the zero-order reaction function.

This paper takes fresh vegetables as the products to be delivered. Their deterioration follows the zero-order reaction function $q\left(t\right)=q_0\phi t$, where $q\left(t\right)$ is the quality at time $t$ and $q_0$ is the initial quality. $\phi$ is the deterioration rate of this type of vegetable, obtained by Wang and Li as 0.0216/h [48].

3.2. Notations

The mathematical symbols used in the model are as follows:

Sets:

$O$: Set of distribution centers $P$: Set of UAVs $C$: Set of delivery points $AS$: Set of all points

Parameters:

$V$: The average velocity of a UAV $V_k$: Unit distance cost of a UAV $k$, $k\in P$ $Q_k$: Maximum load of a UAV $k$, $k\in P$ $D_k$: Maximum mileage of a UAV $k$, $k\in P$ $t_{ij}$: Time from point $i$ to point $j$, $i, j\in AS,i\ne j$ $c_{ij}$: Distance from point $i$ to point $j$, $i, j\in AS,i\ne j$ $d_i$: Demand of customer $i$, $i\in C$ $b_j$: Pickup volume from distribution center $j$, $j\in O$ $q_{kj}^L$: Load of UAV $k$ leaving point $j$, $k\in P$, $j\in AS$ $\phi$: The deteriorate rate, 0.0216/h ${\phi }_i^k$: The deterioration degree of the goods transported by UAV $k$, received by customer $i$, $i\in C$, $k\in P$ $L{T}_i$: The latest delivery time that customer $i$ can accept, $i\in C$

Decision variables:

$X_{ijk}=\left\{\begin{array}{c}1 \mathrm{if} \mathrm{UAV} k \mathrm{flies} \mathrm{from} \mathrm{point} i \mathrm{to} \mathrm{point} j\hfill \\ 0 \mathrm{otherwise}\hfill \end{array}\right.$ $Y_{ik}=\left\{\begin{array}{c}1 \mathrm{if} \mathrm{UAV} k \mathrm{delivers} \mathrm{goods} \mathrm{to} \mathrm{customer} i\hfill \\ 0 \mathrm{otherwise}\hfill \end{array}\right.$ $T_{ki}^A$: The time of UAV $k$ arriving at point $i$ $Z_{ik}=\left\{\begin{array}{c}1 T_{ki}^A>L{T}_i\hfill \\ 0 T_{ki}^A\le L{T}_i\hfill \end{array}\right.$

3.3. The Mathematic Model

3.3.1. Objective Function

(1)

Minimize distance cost:

$$f_1=min{\displaystyle \sum }_{i\in AS}{\displaystyle \sum }_{j\in S}{\displaystyle \sum }_{k\in P}{X}_{ijk}{V}_k{c}_{ij}$$

(1)

(2)

Maximize customer satisfaction:

$$f_2=min{\displaystyle \sum }_{k\in P}{\displaystyle \sum }_{i\in C}{\phi }_i^k{d}_i/{\displaystyle \sum }_{i\in C}{d}_i+{\displaystyle \sum }_{k\in P}{\displaystyle \sum }_{i\in C}{Z}_{ik}\left(T_{ki}^A-L{T}_i\right)/{\displaystyle \sum }_{i\in C}L{T}_i$$

(2)

The objective function (1) minimizes the distance cost, which is a common objective function in such problems, and is still one of the objectives in this paper. The ultimate goal of the objective function (2) is to maximize the customer satisfaction, which is affected by the degree of goods deterioration and the delivery time. The deterioration degree of goods follows the zero-order reaction function. The more serious the deterioration degree, the lower the customer satisfaction. At the same time, the goods delivered within the time required by customers will not affect the customer satisfaction but, once that time is exceeded, the time penalty will be triggered, and the customer satisfaction will be reduced. Therefore, the objective function (2) minimizes these two variables, to maximize the customer satisfaction.

3.3.2. Constraint Function

Formula (3) means that the UAV must leave after completing the order of customer $j$. Equation (4) restricts each customer to being served only once. Constraint (5) requires that the goods transported by the UAV each time cannot exceed its maximum load. Equations (6) and (7) are the calculation methods for the transportation time and the deterioration of goods, respectively. Formula (8) means that the pickup volume of the UAV must be equal to the customer’s demand. Constraint (9) is the load constraint of the UAV [49,50].

$${\displaystyle \sum }_{i\in AS}{X}_{ijk}-{\displaystyle \sum }_{l\in AS}{X}_{jlk}=0$$

(3)

$${\displaystyle \sum }_{k\in P}{Y}_{ik}=1$$

(4)

$${\displaystyle \sum }_{i\in C}{Y}_{ik}{d}_i\le Q_k$$

(5)

$$t_{ij}=c_{ij}/V$$

(6)

$${\phi }_i^k=q_0\phi T_{ki}^A{Y}_{ik}$$

(7)

$${\displaystyle \sum }_{i\in C}{\displaystyle \sum }_{j\in O}{X}_{ijk}{b}_j={\displaystyle \sum }_{j\in O}{\displaystyle \sum }_{l\in C}{X}_{jlk}{d}_i$$

(8)

$${\displaystyle \sum }_{i\in AS}{\displaystyle \sum }_{j\in AS}{c}_{ij}{X}_{ijk}\le D_k$$

(9)

3.4. Model Analysis

When modeling the UAV delivery problem, scholars usually only consider the distance cost. In general, the shorter the distance, the lower the delivery cost. Therefore, the method that can obtain the shortest distance is often adopted. However, many goods need to be delivered at the specified time, such as fresh food and medicine. The delivery being too late may lead to food decay, or the best time to treat patients being missed. Fresh vegetables and fruits, especially, will inevitably decay during transportation. These consequences will lead to a reduction in customer satisfaction, thus affecting the long-term profits of businesses. Therefore, while focusing on the distance cost, this paper establishes a bi-objective model, with customer satisfaction as the object.

Through the study of the two objective functions, it is found that there is a negative correlation between them. To maximize the customer satisfaction, it is necessary to avoid the deterioration of goods as much as possible. This will make the transportation distance longer, and the distance cost higher. On the contrary, when the shortest route is used for transportation, some customers will wait for a long time. This will not only trigger the time penalty, but also lead to the serious deterioration of the goods. Therefore, we considered using the MOEA, which is more conducive to solving bi-objective problems with a negative correlation.

4. SPEA2 and Its Improvement

4.1. Overview of SPEA2

4.1.1. Pareto Optimal Solution Theory

In MOPs, the multiple objectives to be optimized are usually in conflict. Therefore, such problems do not have a specific optimal solution, and cannot be solved directly. Only by optimizing as many objectives as possible, can a compromise solution set be obtained, which is the Pareto optimal solution set. Then, according to different requirements, people can choose the desired solution from this set. Broadly speaking, there is no better solution than the Pareto solution set, so it is the optimal solution set [51].

MOPs can usually be expressed as Equation (10) [52]:

$$\mathrm{Max} \mathrm{or} \min F\left(x\right)=\left[f_1\left(x\right),f_2\left(x\right),\dots ,f_m\left(x\right)\right] \mathrm{s}. \mathrm{t}. x\in \mathsf{\Omega }$$

(10)

The decision variables are represented by $x$. By choosing the correct value of $x$, the ideal solution of the optimization problem can be obtained. The solution in the Pareto solution set cannot result in all the objectives reaching the optimal value but, when one objective is determined, the Pareto solution can allow the other objectives to reach the optimal value. Therefore, when choosing the optimal solution, we can only consider the Pareto solution set. Pareto solutions will degrade other objective functions when improving one objective function, which is completely consistent with the characteristics of MOPs. Pareto solutions are solutions that do not dominate other solutions, so they are also called non-dominated solutions. If there is no artificial choice, they have the same position in MOPs, so they are usually regarded as a whole, which is called the Pareto frontier. If the solution satisfies Equations (11) and (12), it is a Pareto solution [53]:

$$G\left(x^{\prime }\right)\le G\left(x\right)$$

(11)

$$\mathrm{Other} \mathrm{goals}: G_i\left(x^{\prime }\right)\ge G_i\left(x\right)$$

(12)

The Pareto dominance relation is defined as follows: for the minimized multi-objective vector a with $m$ objective components, for any two decision variables $F\left(x\right)=\left[f_1\left(x\right),···,f_m\left(x\right)\right]$ and $x_1,x_2\in X$, if and only if $i\in \left\{1,2,···,m\right\}$, there is $f_i\left(x_1\right)<f_i\left(x_2\right)$, then $x_1$ dominates $x_2$, which is recorded as $x_1>x_2$. In the whole solution space $S$, the Pareto solution set is a set composed of all the non-dominated solutions [54].

The Pareto frontier is defined as follows: a series of non-dominated solutions form Pareto solution set $P$, and the surface composed of their corresponding vectors is Pareto frontier $P_f$ [55].

4.1.2. Strength Pareto Evolutionary Algorithm 2

On the basis of the SPEA (strength Pareto evolutionary algorithm), Zitzler and others proposed SPEA2. Its pseudo code is similar to Algorithm 1 [56,57]:

Algorithm 1: Strength Pareto evolutionary algorithm 2

Input: $T$ (Maximum number of iterations), $n$ (Population number) Output: $P_T^{\ast }$ (non-dominated set) Process: Initialize evolutionary population $P_t$, initialize the external population $P_t^{\ast }$ as an empty set; Assign all individuals in $P_t$ and $P_t^{\ast }$ according to their fitness; Transfer non-dominated individuals from $P_t$ and $P_t^{\ast }$ to external population $P_{t+1}^{\ast }$. If the number of $P_{t+1}^{\ast }$ exceeds the preset capacity of external population, delete these non-dominated individuals; if the number of $P_{t+1}^{\ast }$ does not reach the preset capacity of external population, it is necessary to select the dominant individuals from the evolutionary population and store them in $P_{t+1}^{\ast }$; Judge whether the population evolution conditions are met. If so, $P_{t+1}^{\ast }$ will be output as the result and the algorithm ends; Otherwise, proceed to the next step; Using tournament selection method to select excellent individuals from $P_{t+1}^{\ast }$ as the parent to enter the mating pool; Cross and mutate such individuals to obtain a new evolutionary population $P_{t+1}$, make $t=t+1$, and transfer to step 2.

The SPEA2 flow is shown in Figure 1:

Compared with the SPEA, the SPEA2 features improvements in the population individual evaluation mechanism, and the environment selection mechanism. In order to avoid the evolutionary population dominated by the same external population individual having the same fitness value, the SPEA2 considers the dominant information and non-dominated information of the individual when calculating the individual fitness value. The fitness function is defined as Equation (13) [26]:

$$\left\{\begin{array}{c}f\left(i\right)=R\left(i\right)+D\left(i\right)\\ R\left(i\right)={\displaystyle \sum _{j\in P_t+P_t^{\ast },i<j}}S\left(j\right)\\ S\left(i\right)=\left|\left\{j|j\in P_t\cup P_t^{\ast },j<i\right\}\right|\\ D\left(i\right)=\frac{1}{{\sigma }_i^k+2}\end{array}\right.$$

(13)

In Equation (13), the number of individuals dominated by individual $i$ is expressed by $S\left(i\right)$. $R\left(i\right)$ represents the sum of the number of individuals $j$ dominating $i$ dominating other individuals in the evolutionary population $P_t$ and the external population $P_t^{\ast }$. The dominant individuals are sorted by $R\left(i\right)$, while the non-dominant individuals cannot be sorted by $R\left(i\right)$. Therefore, the individual density information $D\left(i\right)$ is added in the SPEA2, to sort the non-dominant individuals. ${\sigma }_i^k$ expresses the Euclidean distance from individual $i$ to the $kth$ individual [58], where $k=\sqrt{N+M}$. For $k=\sqrt{N+M}$; $N$ is the capacity of the evolutionary population, and $M$ is the capacity of the external population [15]. Through the definition of the fitness function in the above formula, we can not only obtain the dominant information and non-dominant information of individuals, but also understand the delivery of individuals in the population. The fitness calculation is shown in Figure 2 [59]:

In the figure, the hollow point is the non-dominant point, and the solid point is the dominant solution. The red number indicates the number of solutions dominated by the point $i$, which is called $S\left(i\right)$. The black number represents the sum of the solutions dominated by the points that dominate the point $i$, which is called $R\left(i\right)$.

In addition, in terms of the environmental selection mechanism, the capacity of the external population set by the SPEA2 is unchanged. When the selected non-dominated individuals exceed the capacity of the external population, they need to be deleted. The deletion strategy adopted by the SPEA2 is a density-based deletion strategy. The basic principle is that non-dominant individuals in dense areas are easier to delete than those in sparse areas [60]. Through this pruning strategy, compared with the SPEA, the SPEA2 can obtain a better-distributed external population.

4.2. Strategies to Improve SPEA2

The convergence ability of the MOEA will be different with different application problems. For some optimization problems with continuous smooth Pareto boundaries, the MOEA has a strong convergence performance. When the Pareto front of some MOPs is discontinuous and fragmented, the convergence performance of the evolutionary algorithm will not be good enough, which will reduce the accuracy of the optimal solution. The most important factor that leads to this phenomenon lies in the characteristics of MOEAs. As a probabilistic algorithm, the ability of MOEAs to search the whole solution space is more prominent. On this basis, the SPEA2 also expands the search scope, so that the optimization process will not fall into local optimization. This further weakens its local search ability. In addition, this characteristic is especially obvious when the MOP needs to be completed in an irregular search space [61].

To summarize, if the SPEA2 can improve the local search on the basis of maintaining the global search ability, then it can solve this type of problem excellently. Therefore, an improved SPEA2 is proposed. This algorithm mainly solves the problem of insufficient local search ability, and effectively improves its convergence ability to the global optimal solution.

4.2.1. Add a Local Search Policy

After each iteration, the algorithm will generate the non-dominated solution set $Q$. Firstly, the additional population $L$ is set to replicate the elements in $Q$. Then, a neighborhood search is carried out for the individuals in $Q$, to see whether there is a better solution. If such a solution can be found, then $Q$ will be updated. With the continuous repetition of this process in each iteration, the set in $Q$ gradually approaches the Pareto optimal boundary of MOPs. The idea of the local search is shown in Figure 3:

As shown in Figure 3, points $A,B$, and $C$ represent a portion of the population individuals generated after genetic evolution. We assume that point $D$ is the optimal solution, in theory. Obviously, points $A,B$, and $C$ are very close to point $D$, but several iterations are still required to get to point $D$. At this time, if a local search can be carried out on point $B$, then the optimal point $D$ can be easily obtained. Therefore, it is very useful to locally search some points obtained by the algorithm.

If the local search is used for all individuals in th e external archive set $Q$ after each iteration, the amount of calculation will be large. Therefore, $Q$ is sorted according to fitness. Then, a local search is only carried out on excellent individuals. The higher the fitness of individuals, the closer they are to the Pareto optimal boundary, so a local search for them is more meaningful than for any other individual. For a selected non-dominated solution, we take each decision variable value of the individual as the center of the circle, to search the neighborhood near their distinct search radius. The specific operation of a neighborhood search is as follows: we suppose an MOP has two decision variables, which are $X$ and $Y$, respectively. The selected coordinate of an excellent individual is $\left(x_i,y_i\right)$. According to the preset search radius, the search range is shown in Equation (14):

$$\left\{\begin{array}{c}{R}_x=x_i\pm r\\ R_y=y_i\pm r\end{array}\right.$$

(14)

where $R$ stands for the range of neighborhood to be searched, and the value of the search radius $r$ can be determined according to the specific problem. We divide the interval $R$ around each selected individual into $n$ equal parts, and this individual represents this interval. The following Equation (15) is the new individual coordinate set generated after the local search of the individual:

$$\left\{\begin{array}{c}{x}_{ij}=x_i-r+\frac{2rj}{n}\\ y_{ij}=y_i-r+\frac{2rj}{n}\end{array}\right.$$

(15)

where $n$ represents the search density of the neighborhood interval. The higher the density, the stronger the local search capability. However, the corresponding calculation amount will increase exponentially. $x_{ij}$ represents the coordinate of the $j$th point in the neighborhood of the $i$th individual, and $y_{ij}$ has a similar meaning. Through taking the cartesian product of the coordinate sets of the above two variables, a new set of individual coordinate sets can be obtained. These coordinate sets constitute an individual set after a local search.

In most cases, one individual alone cannot achieve the optimal solution via a local search. However, if this operation is performed on multiple individuals in each iteration, and with continuous iterations, the probability of finding the global optimal solution set will increase exponentially. In this way, it can make up for the defects of the original algorithm, in theory.

After the local search, the population $L$ and the evolutionary population $P$ jointly constitute the parent population of the next iteration. Then, the new parent population participates in other genetic operations in the algorithm. The use of the external population exchanges the excellent genes generated by the evolutionary population and the local search population, meaning that the balance state within the population can be broken through. This will help to prevent falling into local optimization, enrich the diversity of the alternative populations, and greatly increase the chance of finding the optimal solution set.

4.2.2. Improved Crossover Operator

In the general genetic algorithm, the crossover operator has a crossover probability. During operation, it is only necessary to recombine the genes of two individuals, according to this probability. The consequence is that similar, or even identical, individuals may be crossed. Obviously, this is meaningless for the population evolution. It will only increase the amount of calculation in the algorithm, cause the population evolution to stagnate, and lead to local optimization. To avoid this phenomenon in the improved algorithm, the two individuals which will be crossed need to be judged according to their similarity before the operation. If their difference is not large enough, they will not be crossed, to avoid meaningless crossing between two identical or similar individuals.

In general, the SPEA2 takes the same gene number of two chromosomes as the similarity of individuals. Although this method is simple, it ignores a very important problem. Two similar chromosomes usually have the same genes at most corresponding positions. However, in binary coding, having the same genes in most positions does not mean that the actual values of the two chromosomes are similar [62]. Therefore, this paper adopts the Euclidean distance as the measure of individual similarity, which is more scientific. The specific operation is as in Equation (16):

$$Crossover=\left\{\begin{array}{c}1, D_{ij}\ge T\\ 0, D_{ij}<T\end{array}\right.$$

(16)

$D_{ij}$ is the Euclidean distance from $i$ to $j$, 1 represents that the cross-operation needs to be carried out, and 0 represents that the values do not need to cross. $T$ stands for the set similarity threshold.

4.3. Steps to Improve SPEA2

The improvement of the algorithm is mainly reflected in the setting of the archive set $Q$, and adding the local search archive set $L$. Because the main scheme of algorithm improvement is a local search for outstanding individuals after each archiving, a special set should be set for the local search. After each local search, a fitness allocation should be carried out for individuals in the three sets, and non-dominant individuals should be retained.

The pseudo code of the improved SPEA2 is shown in Algorithm 2:

Algorithm 2: Improved strength Pareto evolutionary algorithm 2

Input: $T$ (Maximum iteration), $N$ (population), $M$ (archive size) Output: $P_T^{\ast }$ (non-dominated set) Process: Initialize population $P_0$, size $N$, set archive set $Q_0$ and local search archive set $L_0$ as empty set, set $Q_0$ and $L_0$ size $M$, $t=0$ Fitness allocation for all individuals in $P_t$, $Q_t$ and $L_t$ Save all non-dominant individuals in $P_t$, $Q_t$, and $L_t$ into $Q_{t+1}$. If the size of $Q_{t+1}$ is not $M$, the size of $Q_{t+1}$ needs to be adjusted to meet the requirement of delivery Assign $Q_{t+1}$ to $L_{t+1}$ If the termination condition is met, all non-dominant individuals in $Q_{t+1}$ are returned Select, cross and mutate $Q_{t+1}$, save the results to $P_{t+1}$, update some individuals on $P_{t+1}$, perform local search on $L_t$, $t=t+1$, go to Step 2

5. Experimental Results and Analysis

5.1. Experimental Environment

In this chapter, we first use the existing test functions to prove that the improved algorithm has a better optimization ability and convergence. Then, the crossover operator and mutation operator of the improved algorithm are determined by comparison. Finally, several algorithms are used to solve the model, proving that the improved algorithm can solve this type of problem, and provide a better solution.

The experimental platform of this paper is AMD Ryzen 7 5800H, 16.0 GB RAM, and the experimental environment is MATLAB R2021a.

5.2. Algorithm Comparison

To discuss the characteristics of the improved SPEA2, nine common test functions commonly used in MOEAs were selected in this paper. In these nine cases, the improved SPEA2 was compared with four MOEAs: the SPEA2, NSGA2, SPEA2 + SDE (strength Pareto evolution algorithm 2–shift-based density estimation), and SPEA [63]. All the experiments were conducted on Tian’s PlatEMO platform [63]. The experimental parameters were as follows: the initial population size was 50, the number of iterations was 100, the external archive set size was 30, the crossover probability was 0.5, and the mutation probability was 0.5. All the experiments adopted the binary tournament selection method. See

Table 1. The benchmark functions.

	Function	Dim	Range
${\mathit{F}}_{\mathbf{1}}$	$Min F_1\left(x\right)=\left(f_1\left(x\right),f_2\left(x\right)\right)$ $\left\{\begin{array}{c}{f}_1\left(x\right)=x\\ f_2\left(x\right)=g\left(1-\sqrt{f_1/g}\right)\end{array}\right.$ $g\left(x\right)=1+9{\displaystyle \sum _{i=2}^m}{x}_i/\left(m-1\right)$	30	[0, 1]
${\mathit{F}}_{\mathbf{2}}$	$Min F_2\left(x\right)=\left(f_1\left(x\right),f_2\left(x\right)\right)$ $\left\{\begin{array}{c}{f}_1\left(x\right)=x\\ f_2\left(x\right)=g\left(1-{\left(f_1/g\right)}^2\right)\end{array}\right.$ $g\left(x\right)=1+9{\displaystyle \sum _{i=2}^m}{x}_i/\left(m-1\right)$	30	[0, 1]
${\mathit{F}}_{\mathbf{3}}$	$Min F_3\left(x\right)=\left(f_1\left(x\right),f_2\left(x\right)\right)$ $\left\{\begin{array}{c}{f}_1\left(x\right)=x\\ f_2\left(x\right)=g\left(1-\sqrt{\frac{f_1}{g}}-\left(f_1/g\right)sin\left(10\pi f_1\right)\right)\end{array}\right.$ $g\left(x\right)=1+9{\displaystyle \sum _{i=2}^m}{x}_i/\left(m-1\right)$	30	[0, 1]
${\mathit{F}}_{\mathbf{4}}$	$Min F_4\left(x\right)=\left(f_1\left(x\right),f_2\left(x\right)\right)$ $\left\{\begin{array}{c}{f}_1\left(x\right)=1-e^{-4x}si{n}^6\left(6\pi x\right)\\ f_2\left(x\right)=g\left(1-{\left(f_1/g\right)}^2\right)\end{array}\right.$ $g\left(x\right)=1+9\left({\left({\displaystyle \sum _{i=2}^m}\frac{x_i}{m-1}\right)}^{0.25}\right)$	30	[0, 1]
${\mathit{F}}_{\mathbf{5}}$	$Min F_5\left(x\right)=\left(f_1\left(x\right),f_2\left(x\right)\right)$ $\left\{\begin{array}{c}{f}_1\left(x\right)=x\\ f_2\left(x\right)=g\left(1-\sqrt{f_1/g}\right)\end{array}\right.$ $g\left(x\right)=1+10\left(m-1\right)+{\displaystyle \sum _{i=2}^m}\left(x_i^2-10\cos \left(4\pi x_i\right)\right)$	10	[−5, 5]
${\mathit{F}}_{\mathbf{6}}$	$Min F_6\left(x\right)=\left(f_1\left(x\right),f_2\left(x\right)\right)$ $\left\{\begin{array}{c}{f}_1\left(x\right)=x^2\\ f_2\left(x\right)={\left(x-2\right)}^2\end{array}\right.$	100	$[-{10}^5, {10}^5$]
${\mathit{F}}_{\mathbf{7}}$	$Min F_7\left(x\right)=\left(f_1\left(x\right),f_2\left(x\right)\right)$ $\left\{\begin{array}{c}{f}_1\left(x\right)=1-e^{\left(-{\displaystyle \sum _{i=1}^3}\left(x_i-\frac{1}{\sqrt{3}}\right)^2\right)}\\ f_2\left(x\right)=1-e^{\left(-{\displaystyle \sum _{i=1}^3}\left(x_i+\frac{1}{\sqrt{3}}\right)^2\right)}\end{array}\right.$	3	[−4, 4]
${\mathit{F}}_{\mathbf{8}}$	$Min F_8\left(x\right)=\left(f_1\left(x\right),f_2\left(x\right)\right)$ $\left\{\begin{array}{c}{f}_1\left(x\right)=x_1\\ f_2\left(x\right)=\left(1+10x_2\right)\left[1-{\left(\frac{x_1}{1+10x_2}\right)}^2-\frac{x_1}{1+10x_2}sin\left(8\pi x_1\right)\right]\end{array}\right.$	20	[0, 1]
${\mathit{F}}_{\mathbf{9}}$	$Min F_9\left(x\right)=\left(f_1\left(x\right),f_2\left(x\right)\right)$ $\left\{\begin{array}{c}{f}_1\left(x\right)={\displaystyle \sum _{i=1}^{m-1}}\left(-10e^{-0.2\sqrt{x_i^2+x_{i+1}^2}}\right)\\ f_2\left(x\right)={\displaystyle \sum _{i=1}^m}\left({\left|x_i\right|}^{0.8}+5sin{x}_i^3\right)\end{array}\right.$	3	[−5, 5]

for all test functions [6]:

The optimal solution sets of the nine Benchmark functions are shown in Figure 4:

Figure 4 shows the Pareto optimal solution set of five functions under nine test functions. Under the same conditions, the delivery of the Pareto optimal solution of the SPEA2 and SPEA is messy, which is far from the theoretical optimal solution. Although the results of the NSGA2 and SPEA2 + SDE are better than the above two algorithms, they are still worse than the improved SPEA2. It can be seen that the optimal solution set of the improved SPEA2 is an ideal set, which is evenly distributed on the whole Pareto boundary. Compared with the other four algorithms, its results are closer to the Pareto front, and there is no phenomenon of aggregation and clustering in the optimal set points. This is mainly because the local search can search strongly near the solution set, and can achieve a more uniform solution set than the traditional algorithm.

In previous studies, scholars have proposed some indicators to measure the performance of an MOEA. This paper selects three representative ones [52]:

(1)

VN: The numbers of times that every algorithm can run successfully and obtain results [64]. According to these data, the ability of algorithms to find results can be understood. Generally speaking, algorithms that cannot find useful solutions stably and efficiently will be abandoned.

(2)

TT100%: Sayyad et al. proposed this index in 2013 [65]; it represents the time at which the algorithm reaches TT100%. This can be used to express the speed at which the algorithm converges with the Pareto front.

(3)

Hypervolume (HV): A common performance indicator in the MOEA field. The Pareto front dominates a solution space $P_f$, and $R={\left(r_1,r_2, . . .,r_n\right)}^T$ is a point in $P_f$. According to the definition of Zitzler [16], the HV represents the volume of $P_f$. Taking $R$ as the boundary, $HV\left(A\right)$ can be obtained using Formula (17):

$$HV\left(A\right)=V\left(U_{X\in A}\left[f_1\left(x\right),r_1\left]\times \right[f_2\left(x\right),r_2\left]\times ···\times \right[f_n\left(x\right),r_n\right]\right)$$

(17)

$n$ is the objective number, $F=\left(f_1\left(x\right), f_2\left(x\right), \dots , f_n\left(x\right)\right)$ is the objective function in $P_f$, and $V\left(S\right)$ represents the Lebesgue measure [66] of $P_f$. A large value of HV represents that the algorithm performs better in convergence, diversity, and uniformity.

The VN, TT100%, and HV of five algorithms under nine test functions are shown in

Table 2. The VN, TT100% and HV of five algorithms under nine test functions.

Functions	Index	Improved SPEA2	NSGA2	SPEA2	SPEA2 + SDE	SPEA
${\mathit{F}}_{\mathbf{1}}$	VN (20)	20	20	20	20	20
	TT100%	1.0	1.1	7.3	0.8	10.5
	HV	0.8934	0.8718	0.8439	0.8792	0.7957
${\mathit{F}}_{\mathbf{2}}$	VN (20)	20	20	20	20	20
	TT100%	1.5	1.6	9.7	1.2	14.4
	HV	0.8859	0.8731	0.8392	0.8917	0.8038
${\mathit{F}}_{\mathbf{3}}$	VN (20)	20	20	20	20	20
	TT100%	1.9	1.8	15.6	1.9	21.4
	HV	0.8126	0.7933	0.7279	0.8018	0.6817
${\mathit{F}}_{\mathbf{4}}$	VN (20)	20	20	20	20	20
	TT100%	2.6	2.7	24.8	2.4	35.7
	HV	0.8596	0.7151	0.6139	0.8638	0.5884
${\mathit{F}}_{\mathbf{5}}$	VN (20)	20	20	18	20	12
	TT100%	21.9	26.4	96.1	20.3	183.7
	HV	0.8394	0.7062	0.6028	0.8045	0.4119
${\mathit{F}}_{\mathbf{6}}$	VN (20)	20	20	20	20	20
	TT100%	7.9	8.3	41.9	8.2	67.5
	HV	0.8783	0.8526	0.8313	0.8617	0.8106
${\mathit{F}}_{\mathbf{7}}$	VN (20)	20	20	20	20	20
	TT100%	25.7	28.5	136.2	24.3	207.3
	HV	0.7908	0.7347	0.7214	0.7705	0.7008
${\mathit{F}}_{\mathbf{8}}$	VN (20)	20	20	20	20	20
	TT100%	3.0	2.8	17.1	2.4	37.6
	HV	0.8769	0.8573	0.8491	0.8707	0.8327
${\mathit{F}}_{\mathbf{9}}$	VN (20)	20	20	20	20	20
	TT100%	11.9	13.2	50.4	10.6	118.1
	HV	0.7839	0.7654	0.7486	0.7792	0.7218

:

In Table 2, in the 20 experiments for each test function, the VN value of almost all the algorithms is 20. This shows that nearly all the algorithms can obtain the answer every time. From the perspective of the TT100%, SPEA2 + SDE usually has the smallest value, indicating that it is the fastest to obtain the optimal solution. However, the value of the improved SPEA2 is only a little larger than the minimum value, which is greatly improved compared with the SPEA2 and SPEA. Finally, in eight of the nine test functions, the improved SPEA2 has the largest HV value. This shows that it is superior to the other four algorithms in many aspects.

In addition, this paper also selects the test functions $F_1$, $F_2$, $F_3$, and $F_4$ to compare the convergence of five algorithms. The convergence distance is the Euclidean distance between a point in the population after each iteration, and the point before the iteration.

In Figure 5, the convergence speed of the NSGA2 is very fast in four cases, but it is extremely unstable, and fluctuates greatly. The SPEA2 and SPEA usually show no obvious convergence trend at 100 iterations. Although the SPEA2 + SDE can converge faster, the convergence speed is usually slower than that of the improved SPEA2. Therefore, the improved SPEA2 is also better than the other four algorithms in the convergence speed and stability, which reflects the advantages of the improved algorithm.

5.3. Parameter Selection

The parameters that affect the improved SPEA2 are mainly the crossover operator and mutation operator. Next, the classical benchmark functions $F_1$, $F_2$, and $F_3$ will be taken as examples, to determine the values of these two parameters.

In order to compare these two parameters as widely as possible between 0 and 1, five cases were selected: 0.1, 0.3, 0.5, 0.7, and 0.9. We conducted orthogonal experiments on the crossover operator and mutation operator. Firstly, when the mutation operator was 0.1, 0.3, 0.5, 0.7, or 0.9, experiments were conducted on the crossover operators of 0.1, 0.3, 0.5, 0.7, and 0.9, respectively. Then, we compared the best crossover operator under the five mutation operators, to obtain the best combination. The following figures show the average results after 20 experiments.

Figure 6a–e show the comparison of the crossover operators, when the mutation operator was 0.1, 0.3, 0.5, 0.7, and 0.9, respectively. It can be seen that, regardless of the mutation operator, when the crossover operator was 0.5, the result was the best. So, a comparative experiment was conducted on the best results from the first five scenarios (Figure 6f). It can be seen that when the mutation operator and crossover operator were both 0.5, the results were the best.

Next, we performed the same experiment on $F_2$ and $F_3$.

In Figure 7 and Figure 8, we can see that the improved SPEA2 can achieve the best results when the mutation operator is 0.5 and the crossover operator is 0.5. Therefore, the mutation operator and crossover operator are both selected as 0.5.

5.4. Simulation Experiment

In Section 5.2, the improved algorithm was tested on nine benchmark functions. From those experiments, we can see that the improved algorithm has certain advantages. In addition, we also wanted to verify whether the improved algorithm still achieved better results in the established UAV delivery model. Therefore, a certain number of distribution centers and delivery points were randomly generated in the city. We then compared the improved algorithm with the SPEA2, NSGA2, SPEA2 + SDE, and SPEA in the scene.

We set the crossover probability as 0.5, mutation probability as 0.5, population as 100, iteration times as 100, and test times as 20. Figure 8, below, shows the average results of the 20 experiments on these five algorithms, under the UAV delivery model in this paper.

Figure 9 shows the experimental results when the number of tasks is 30 and 100, respectively. The two factors that affect customer satisfaction are actually related to time. The goods decay over time, and late delivery will also have an impact. Therefore, the ‘Time cost’ on the vertical axis does not mean that time is a goal, but refers to the customer satisfaction. The abscissa is the other object, namely the ‘Distance cost’. From the perspective of the algorithm, the improved algorithm has better Pareto frontiers. It has better continuity, and is closer to the origin. This indicates that the improved algorithm is stable, and gives better results. On the other hand, from the perspective of the simulation, when the total distance of the UAVs is the same, the smaller the value of the improved SPEA2 in the vertical coordinate, the higher the customer satisfaction. When the customer satisfaction is the same, the path planned by the improved SPEA2 is shorter. This comparison is especially evident when the number of tasks is 100; for example, when the distance cost is 1, the customer satisfaction of the improved algorithm is more than 20% higher than that of other algorithms. That is to say, the improved algorithm always yields more benefits, regardless of the options chosen by the decision-maker.

To illustrate the comparison of the five algorithms, the generation distance proposed by Deb [18], and the spacing index proposed by Schott [67] are used in the quantitative analysis of the algorithms. The comparison results are shown in Figure 10 and Figure 11:

Their 95% confidence intervals are shown in

Table 3. The 95% confidence intervals for the generational distance values and spacing values of the five algorithms.

Index	Improved SPEA2	NSGA2	SPEA2	SPEA2 + SDE	SPEA
Generational distance	[1.34 × 10−3, 2.60 × 10−3]	[2.34 × 10−3, 5.43 × 10−3]	[1.09 × 10−2, 1.56 × 10−2]	[1.81 × 10−3, 3.48 × 10−3]	[1.43 × 10−2, 1.92 × 10−2]
Spacing	[8.55 × 10−4, 1.93 × 10−3]	[2.00 × 10−3, 3.97 × 10−3]	[3.27 × 10−3, 5.10 × 10−3]	[1.34 × 10−3, 2.60 × 10−3]	[4.40 × 10−3, 6.31 × 10−3]

:

Figure 10 and Figure 11 are box diagrams, drawn according to the values of the generational distance and spacing that were obtained by testing the objective function 20 times, and with the number of tasks as 30. All the box diagrams were drawn using Origin software. It can be seen that the values of the generation distance and spacing obtained by the improved SPEA2 are always less than those obtained by the SPEA2, NSGA2, SPEA2 + SDE, and SPEA. The average generational distance value of the improved algorithm is about 30% smaller than the second one, and the average spacing value is about 20% smaller. At the same time, the volatility of the improved SPEA2 is relatively low, which shows that it has advantages and stability in solving the problems of UAV delivery.

6. Conclusions

With the rise of the e-commerce industry, the UAV delivery of goods has gradually become more and more a delivery method that is discussed in daily life. It can reduce the dependence on manpower, reduce the transportation distance, and plan flexible routes. In addition, with the increasing demand for perishables, it has become an important issue to provide customers with goods as fresh as possible. There is no doubt that vehicles with a better refrigeration effect can be used for transportation; that is, keeping perishables fresh as much as possible during transportation. However, this paper aimed to improve the customer satisfaction without changing the transport vehicle, which required planning a more reasonable path. Therefore, a bi-objective model was established, with the distance cost and customer satisfaction as the objectives. Based on analysis of the model, the SPEA2 was adopted and improved. When improving the SPEA2, the idea of a local search was introduced, to determine the dominant and non-dominant solutions. In this way, on the basis of ensuring the excellent global search ability of the SPEA2, its poor local search ability was also strengthened. The improved algorithm not only maintains the diversity of the solution set, but also obtains results closer to the Pareto front. It shows a good performance, in theory and in experiments. Therefore, it has a high application value in UAV cargo delivery, and other MOPs.

(1)

The improved SPEA2, SPEA2, NSGA2, SPEA2 + SDE, and SPEA were tested through nine classical test functions. Compared with the other four algorithms, the convergence of the improved SPEA2, and the stability of the Pareto solution set are significantly improved.

(2)

In the empirical analysis, the improved SPEA2 is nearer to the theoretical Pareto front than the SPEA2, NSGA2, SPEA2 + SDE, and SPEA, and the result delivery is more uniform and stable. The convergence and stability of the algorithm were verified according to the aspects of the Pareto front, generational distance, and spacing. Through comprehensive analysis, the improved SPEA2 has been made more effective in dealing with the delivery of UAVs with hard time windows, and can achieve better results.

Although the improved SPEA2 can effectively solve the bi-objective UAV delivery problem, it failed to perform well on models with more objectives. We will further improve the multi-objective algorithm so that it can perform well when facing problems with more than five objectives. On the other hand, the effect of the traffic conditions on the problem is not taken into account in the model of this paper, because it involves a dynamic environment. In our follow-up work, we will try to incorporate the dynamic environment, to solve more complex problems.