# Probabilistic Analysis of Search Performance of Differential Evolution Algorithm in Low-Dimensional Case

Authors: Congjiao Wang; Yuhan Wei; Yanlin Hou; Wenjia Wu; Cong Huang; Jing Chen
Journal: Wireless Communications and Mobile Computing
 (2022)

---

## Abstract

In this paper, a method of establishing and analyzing the probability model of its mutation and crossoperation process for differential evolution (DE) algorithm is proposed; especially, the escape ability and further optimization ability of individuals trapped in the local optimal neighborhood are deduced in detail, and the characteristic curves of the influence of key parameters such as population size, scaling factor, and crossfactor on the search performance of the algorithm are obtained. It provides a theoretical reference for the application of the algorithm.


---

## Body

## 1. Introduction

As a rising star in intelligent optimization algorithms, differential evolution (DE) algorithm [1] has the advantages of simple structure, few controlled parameters, and high parallelism. It has been successfully applied to various scientific research and practical engineering fields such as data mining, pattern recognition, signal processing, control, and scheduling. The research on its theory and application has important academic significance and engineering value [2].The standard version of DE algorithm was first proposed by Rainer Storn and Kenneth Price of Berkeley University in 1995. Because of its unique differential mutation operator, it was originally used to solve the Chebyshev polynomial fitting problem, but it was later found that it is also very effective in solving complex optimization problems. In 2004, Vesterstrom and Thomsen conducted an in-depth comparative study of DE algorithm and other intelligent optimization methods by using 34 classical test functions in the literature [3]. The results show that the performance of DE algorithm is the best. Standard DE also adopts basic genetic operations such as mutation, hybridization, and selection, but it reduces the complexity of operation and absorbs the idea of swarm intelligence algorithm using swarm information sharing, so as to further improve the global optimization efficiency of the algorithm.In order to further optimize its performance and expand its application, researchers in various countries have carried out a series of improvement research, but only a few scholars have carried out relevant theoretical analysis on the optimization mechanism of DE algorithm. Based on the knowledge of probability distribution, Zaharie analyzed the influence of the CR value on different crossover operators. Dasgupta established the dynamic system model of standard DE and tried to analyze its stability and convergence speed theoretically. It shows that the DE algorithm is similar to the classical gradient descent method, so it can quickly approach the global optimal solution. Generally speaking, there are few theoretical analysis attempts on DE algorithm in academic circles, and the results are very limited.In this paper, a method of establishing and analyzing the probability model of its mutation and crossoperation transient process is proposed; especially, the escape ability and further optimization ability of individuals trapped in the local optimal neighborhood are deduced in detail, and the characteristic curves of the influence of key parameters such as population size, scaling factor, and crossfactor on the search performance of the algorithm are obtained, which provides a theoretical reference for the application of the algorithm.


### 1.1. Research Status

At present, the researchers of DE algorithm have carried out further research on it mainly from the selection of control parameters, the design of difference strategy, parallel distributed computing, and the combination with other algorithms and promoted the development of DE algorithm in complex environment.


#### 1.1.1. Research on Control Parameters

One of the advantages of DE algorithm is less control parameters. In the standard DE algorithm, there are only three control parameters: population size NP, scaling factorF, and crossover probability factor CR. The traditional parameter setting is mostly based on experience and adopts the trial and error method to select the appropriate value, but this method is time-consuming and labor-consuming and has the problems of difficult decision-making and strain. Literatures [4, 5] analyze the influence of various parameters on the optimization performance of the algorithm through a large number of experimental studies and give useful setting ranges or recommended values. However, there are differences in the conclusions of various literatures, which cannot well deal with specific practical problems. Literature [6] pointed out that the algorithm parameters are sensitive to problems, and different settings are required for problems with different properties and characteristics. Unreasonable parameters will lead to premature or stagnation.Therefore, many scholars have proposed dynamic or adaptive parameter control strategies. Reference [7] proposed an adaptive DE algorithm based on fuzzy principle, which uses fuzzy logic controller to adaptively control F and CR. Reference [8] proposed a control strategy based on the change of fitness value, which dynamically adjusts F by using the ratio of individual maximum fitness value to minimum fitness value in the optimization process, so as to adjust the step size according to the search information. Literature [9] uses mutation operator similar to DE algorithm to adaptively control F. Document [10] proposed a DE algorithm with dynamic adaptive population size, which uses two different coding mechanisms: relative and absolute to adaptively control the population size. In reference [11], it is set to be fixed, while F and CR are directly encoded into the solution individual to participate in the evolutionary operation. By testing some standard functions and comparing with other algorithms, the author verifies the efficiency of this method.Literature [12] proposed JADE algorithm. Each individual corresponds to an F subject to Cauchy distribution and a CR subject to normal distribution, records the success value of making the subindividual enter the next generation smoothly, and then adaptively adjusts the parameters of the two distributions by means of connection update, so as to guide the direction of parameter adaptive adjustment. Experiments show that this method effectively improves the robustness and convergence of the algorithm. According to the information in the evolution process, SaDE algorithm in document [13] adopts different adaptive mechanisms for parameters at different evolution times and different differential strategies, so that the algorithm has achieved good optimization results.


#### 1.1.2. Research on Differential Strategy

Differential mutation operator is the core operator and main feature of DE algorithm. Price and Storn have proposed more than ten different difference operators to realize mutation operation. Then, researchers have carried out rich research on operator design and improvement. Reference [14] gives a trigonometric mutation operator for DE algorithm, which regards the individual as the center point of a supertriangle and moves along the three sides of the supertriangle composed of three groups of weighted difference vectors in different steps to produce new mutated individuals. Literature [15] proposed a mutation operator based on polynomials for multiobjective optimization problems. Reference [16] designed a new DE/target-to-best operator based on domain search in DEGL algorithm, which is weighted combined with traditional mutation operator to strengthen the local optimization ability of the algorithm. The JADE algorithm mentioned above makes full use of the information of the local optimal solution to construct the mutation operator. At the same time, it introduces an archive set for saving the inferior solution to realize the two mutation operations of archiving and nonarchiving. The results show that the algorithm is effective in solving high-dimensional and multimode problems. In recent years, inspired by parameter adaptation, literatures [17–19] have proposed the adaptive mechanism of differential mutation strategy. Its main purpose is to adaptively select the most appropriate mutation strategy from the pool of mutation strategies at different times.


### 1.2. Current Problems

The DE algorithm still needs further discussion and deepening in the following aspects:(1)
Theoretical researchThe theoretical research on convergence, convergence speed, parameter robustness, and iteration termination criterion of the DE algorithm is quite limited. The main research method is verified by experiments, and there is a lack of generally applicable algorithm convergence theory. In order to deeply explore the internal optimization mechanism of DE and improve the optimization performance of DE, it is necessary to speed up its theoretical research.(2)
Improvement researchLike other evolutionary algorithms, the DE algorithm is prone to premature convergence or search stagnation. The main reasons are as follows:(1)
The DE algorithm is sensitive to parameter setting. Whether the parameters are appropriate or not will affect the optimization efficiency of the algorithm, and the selection of parameters is closely related to the nature of specific problems and different search stages(2)
Selection of mutation strategy is difficult. Different mutation strategies will have different guiding effects on the search process of the algorithm. It is difficult to select the best mutation strategy according to the actual problems(3)
Insufficient local search capability. The DE algorithm is based on population and random search mode, which makes it have strong global search ability, but the algorithm itself lacks local search power, which leads to slow convergence speed in the later stage of evolution, and it is difficult to quickly converge to the optimal solution of the problem with less evaluation times of fitness function. This is due to the simple structure, insufficient use of information, and lack of consideration of the analytical properties of the objective function of evolutionary algorithms, including the DE algorithm. These defects are also the unavoidable cost of evolutionary algorithmsThe existing improved DE has improved the overall performance of the algorithm to a certain extent, but these improvements have different characteristics. In the face of some practical application problems, the optimization results are still not satisfactory. Therefore, more efficient and systematic research is needed on how to design parameter strategies and construct improved operators that balance global search and local search and various hybrid forms of DE. At the same time, in view of the complexity of practical problems, paying attention to the algorithm research in complex environment is of great significance to engineering practice.


## 1.1. Research Status

At present, the researchers of DE algorithm have carried out further research on it mainly from the selection of control parameters, the design of difference strategy, parallel distributed computing, and the combination with other algorithms and promoted the development of DE algorithm in complex environment.


### 1.1.1. Research on Control Parameters

One of the advantages of DE algorithm is less control parameters. In the standard DE algorithm, there are only three control parameters: population size NP, scaling factorF, and crossover probability factor CR. The traditional parameter setting is mostly based on experience and adopts the trial and error method to select the appropriate value, but this method is time-consuming and labor-consuming and has the problems of difficult decision-making and strain. Literatures [4, 5] analyze the influence of various parameters on the optimization performance of the algorithm through a large number of experimental studies and give useful setting ranges or recommended values. However, there are differences in the conclusions of various literatures, which cannot well deal with specific practical problems. Literature [6] pointed out that the algorithm parameters are sensitive to problems, and different settings are required for problems with different properties and characteristics. Unreasonable parameters will lead to premature or stagnation.Therefore, many scholars have proposed dynamic or adaptive parameter control strategies. Reference [7] proposed an adaptive DE algorithm based on fuzzy principle, which uses fuzzy logic controller to adaptively control F and CR. Reference [8] proposed a control strategy based on the change of fitness value, which dynamically adjusts F by using the ratio of individual maximum fitness value to minimum fitness value in the optimization process, so as to adjust the step size according to the search information. Literature [9] uses mutation operator similar to DE algorithm to adaptively control F. Document [10] proposed a DE algorithm with dynamic adaptive population size, which uses two different coding mechanisms: relative and absolute to adaptively control the population size. In reference [11], it is set to be fixed, while F and CR are directly encoded into the solution individual to participate in the evolutionary operation. By testing some standard functions and comparing with other algorithms, the author verifies the efficiency of this method.Literature [12] proposed JADE algorithm. Each individual corresponds to an F subject to Cauchy distribution and a CR subject to normal distribution, records the success value of making the subindividual enter the next generation smoothly, and then adaptively adjusts the parameters of the two distributions by means of connection update, so as to guide the direction of parameter adaptive adjustment. Experiments show that this method effectively improves the robustness and convergence of the algorithm. According to the information in the evolution process, SaDE algorithm in document [13] adopts different adaptive mechanisms for parameters at different evolution times and different differential strategies, so that the algorithm has achieved good optimization results.


### 1.1.2. Research on Differential Strategy

Differential mutation operator is the core operator and main feature of DE algorithm. Price and Storn have proposed more than ten different difference operators to realize mutation operation. Then, researchers have carried out rich research on operator design and improvement. Reference [14] gives a trigonometric mutation operator for DE algorithm, which regards the individual as the center point of a supertriangle and moves along the three sides of the supertriangle composed of three groups of weighted difference vectors in different steps to produce new mutated individuals. Literature [15] proposed a mutation operator based on polynomials for multiobjective optimization problems. Reference [16] designed a new DE/target-to-best operator based on domain search in DEGL algorithm, which is weighted combined with traditional mutation operator to strengthen the local optimization ability of the algorithm. The JADE algorithm mentioned above makes full use of the information of the local optimal solution to construct the mutation operator. At the same time, it introduces an archive set for saving the inferior solution to realize the two mutation operations of archiving and nonarchiving. The results show that the algorithm is effective in solving high-dimensional and multimode problems. In recent years, inspired by parameter adaptation, literatures [17–19] have proposed the adaptive mechanism of differential mutation strategy. Its main purpose is to adaptively select the most appropriate mutation strategy from the pool of mutation strategies at different times.


## 1.1.1. Research on Control Parameters

One of the advantages of DE algorithm is less control parameters. In the standard DE algorithm, there are only three control parameters: population size NP, scaling factorF, and crossover probability factor CR. The traditional parameter setting is mostly based on experience and adopts the trial and error method to select the appropriate value, but this method is time-consuming and labor-consuming and has the problems of difficult decision-making and strain. Literatures [4, 5] analyze the influence of various parameters on the optimization performance of the algorithm through a large number of experimental studies and give useful setting ranges or recommended values. However, there are differences in the conclusions of various literatures, which cannot well deal with specific practical problems. Literature [6] pointed out that the algorithm parameters are sensitive to problems, and different settings are required for problems with different properties and characteristics. Unreasonable parameters will lead to premature or stagnation.Therefore, many scholars have proposed dynamic or adaptive parameter control strategies. Reference [7] proposed an adaptive DE algorithm based on fuzzy principle, which uses fuzzy logic controller to adaptively control F and CR. Reference [8] proposed a control strategy based on the change of fitness value, which dynamically adjusts F by using the ratio of individual maximum fitness value to minimum fitness value in the optimization process, so as to adjust the step size according to the search information. Literature [9] uses mutation operator similar to DE algorithm to adaptively control F. Document [10] proposed a DE algorithm with dynamic adaptive population size, which uses two different coding mechanisms: relative and absolute to adaptively control the population size. In reference [11], it is set to be fixed, while F and CR are directly encoded into the solution individual to participate in the evolutionary operation. By testing some standard functions and comparing with other algorithms, the author verifies the efficiency of this method.Literature [12] proposed JADE algorithm. Each individual corresponds to an F subject to Cauchy distribution and a CR subject to normal distribution, records the success value of making the subindividual enter the next generation smoothly, and then adaptively adjusts the parameters of the two distributions by means of connection update, so as to guide the direction of parameter adaptive adjustment. Experiments show that this method effectively improves the robustness and convergence of the algorithm. According to the information in the evolution process, SaDE algorithm in document [13] adopts different adaptive mechanisms for parameters at different evolution times and different differential strategies, so that the algorithm has achieved good optimization results.


## 1.1.2. Research on Differential Strategy

Differential mutation operator is the core operator and main feature of DE algorithm. Price and Storn have proposed more than ten different difference operators to realize mutation operation. Then, researchers have carried out rich research on operator design and improvement. Reference [14] gives a trigonometric mutation operator for DE algorithm, which regards the individual as the center point of a supertriangle and moves along the three sides of the supertriangle composed of three groups of weighted difference vectors in different steps to produce new mutated individuals. Literature [15] proposed a mutation operator based on polynomials for multiobjective optimization problems. Reference [16] designed a new DE/target-to-best operator based on domain search in DEGL algorithm, which is weighted combined with traditional mutation operator to strengthen the local optimization ability of the algorithm. The JADE algorithm mentioned above makes full use of the information of the local optimal solution to construct the mutation operator. At the same time, it introduces an archive set for saving the inferior solution to realize the two mutation operations of archiving and nonarchiving. The results show that the algorithm is effective in solving high-dimensional and multimode problems. In recent years, inspired by parameter adaptation, literatures [17–19] have proposed the adaptive mechanism of differential mutation strategy. Its main purpose is to adaptively select the most appropriate mutation strategy from the pool of mutation strategies at different times.


## 1.2. Current Problems

The DE algorithm still needs further discussion and deepening in the following aspects:(1)
Theoretical researchThe theoretical research on convergence, convergence speed, parameter robustness, and iteration termination criterion of the DE algorithm is quite limited. The main research method is verified by experiments, and there is a lack of generally applicable algorithm convergence theory. In order to deeply explore the internal optimization mechanism of DE and improve the optimization performance of DE, it is necessary to speed up its theoretical research.(2)
Improvement researchLike other evolutionary algorithms, the DE algorithm is prone to premature convergence or search stagnation. The main reasons are as follows:(1)
The DE algorithm is sensitive to parameter setting. Whether the parameters are appropriate or not will affect the optimization efficiency of the algorithm, and the selection of parameters is closely related to the nature of specific problems and different search stages(2)
Selection of mutation strategy is difficult. Different mutation strategies will have different guiding effects on the search process of the algorithm. It is difficult to select the best mutation strategy according to the actual problems(3)
Insufficient local search capability. The DE algorithm is based on population and random search mode, which makes it have strong global search ability, but the algorithm itself lacks local search power, which leads to slow convergence speed in the later stage of evolution, and it is difficult to quickly converge to the optimal solution of the problem with less evaluation times of fitness function. This is due to the simple structure, insufficient use of information, and lack of consideration of the analytical properties of the objective function of evolutionary algorithms, including the DE algorithm. These defects are also the unavoidable cost of evolutionary algorithmsThe existing improved DE has improved the overall performance of the algorithm to a certain extent, but these improvements have different characteristics. In the face of some practical application problems, the optimization results are still not satisfactory. Therefore, more efficient and systematic research is needed on how to design parameter strategies and construct improved operators that balance global search and local search and various hybrid forms of DE. At the same time, in view of the complexity of practical problems, paying attention to the algorithm research in complex environment is of great significance to engineering practice.


## 2. Introduction to Standard Differential Evolution Algorithm

### 2.1. Standard Process of Differential Evolution Algorithm

Standard DE is essentially a special genetic algorithm based on real number coding and optimization preserving greedy strategy. Its principle and structure are similar to genetic algorithm. Similarly, starting from the randomly initialized population, genetic operations such as mutation, crossover, and selection are adopted, so as to solve the optimization problem iteratively.Constrained optimization problems are usually expressed as(1)minfxs.t.cx=0,i∈Ecx≥0,i∈I,where E and I are the index sets of equality and inequality constraints, respectively, and cix are constraint functions. The standard process of solving the above general optimization problems by DE is described as below.
(1)
Population initializationAssuming that the population size is NP (i.e., there are NP individuals participating in evolution) and the evolutionary algebraG=0,1,⋯,Gmax, the ith individual in the population is expressed as
(2)Xi,G=x1,i,G,x2,i,G,x3,i,G,⋯,xD,i,G,where D represents the dimension of individual variables. In the initialization stage, in order to make individuals cover the whole search space as much as possible, each individual in the initial population (i.e., G=0) is randomly and evenly generated within the value range Xmin,Xmax specified by the objective function:
(3)xj,i,0=xj,min+randi,j0,1⋅xj,max−xj,min,where xj,max and xj,min are the upper and lower bounds of the jth dimension of individual variable Xi,G, respectively. randi,j0,1 is a random number uniformly distributed in the interval [0,1].
(2)
Mutation operationMutation is the core operation of the DE algorithm. The purpose is to generate new candidate individuals, enhance the diversity of the population, and guide individuals to a favorable search direction. Compared with other evolutionary algorithms, the mutation operation of DE makes full use of population information, which generates mutation vectora based on the difference vector between different individuals. DE/a/b is usually used to represent different mutation strategies, in which a represents the selection mode of basis vector, often including “rand” and “best,” representing the basis vector selected randomly and optimally, respectively; b represents the number of difference vectors. At present, the commonly used mutation strategies areDE/rand/1:
(4)Vi,G=Xr1,G+F⋅Xr2,G−Xr3,G,DE/best/1:
(5)Vi,G=Xbest,G+F⋅Xr2,G−Xr3,G,DE/best/2:
(6)Vi,G=Xbest,G+F⋅Xr1,G−Xr2,G+F⋅Xr3,G−Xr4,G,DE/current−to−best/1:
(7)Vi,G=Xi,G+F⋅Xbest,G−Xi,G+F⋅Xr1,G−Xr2,G,where r1,r2,r3,r4∈1,2,⋯,NP is a random integer that is not equal to i and different from each other, representing the index number of different individuals in the population; Xbest,G is the individual with the best fitness in generation G population; F∈0,2 is the scaling factor, which is used to control the search step size, and its value directly affects the convergence performance of the algorithm.
(3)
CrossoperationThe purpose of crossover is to disturb the individuals in the original population, so as to enhance the exploration of local areas. Binomial crossover is often used in standard DE, and the target vectorXi,G and variation vector Vi,G are discretely crossed in each dimension, resulting in test vectorUi,G=u1,i,G,u2,i,G,⋯,uD,i,G:
(8)uj,i,G=vj,i,G,randj≤CR⋯or⋯j=randni,xj,i,G,randj>CR⋯and⋯j≠randni,where randj∈0,1 is a uniformly distributed random number; randni∈1,2,⋯,D is the randomly selected dimension variable index to ensure that at least one bit of Ui,G is contributed by Vi,G, while for other bits, it is determined by the crossover probability factor CR; CR∈0,1; its value is often related to the specific nature of the problem (such as multimode problem or single-mode problem and independent variable problem or independent variable interdependence problem).
(4)
Select operationBased on the natural selection mechanism of fittest survival, standard DE selects individuals with better fitness to reproduce under the condition of keeping NP unchanged, so as to realize the inheritance of excellent characteristics, and then guide the population to the optimal region and gradually approach the optimal solution. The selection policy can be expressed as(9)Xi,G+1=Ui,G,fUi,G<fXi,G,Xi,G,fUi,G≥fXi,G,where Xi,G+1 is the ith individual of the next-generation population and f is the fitness function. Obviously, the selection operation ensures that the new individual is at least no worse than the original individual and has the function of memorizing the individual optimal solution.
(5)
Algorithm terminationAfter mutation, crossover, and selection, DE completes an update of the population. At this time, it needs to judge whether it meets the operation termination conditions. If satisfied, the current optimal solution is output; otherwise, continue the search. The conditions for determining the termination of the algorithm can be as follows:(a)
The fitness function value of the optimal individual in the two generations of the population is within a set range in successive iterations(b)
Reach the maximum evolutionary algebra set by the algorithmTo sum up, the overall flow of the standard DE algorithm is shown in Figure1.Figure 1
The flowchart of the standard DE algorithm.


### 2.2. Basic Characteristics of Differential Evolution Algorithm

Compared with other evolutionary algorithms, DE has its own uniqueness:(1)
The structure is simple and easy to implement: the control parameters of DE are relatively few, which reduces the difficulty of parameter selection and adjustment in practical application. At the same time, its core difference mutation operator only performs addition and subtraction operations between individual vectors, so it is easy to use(2)
Using real number encoding: compared with binary encoding of GA algorithm, real number encoding is very suitable for dealing with optimization problems in continuous space(3)
Selection operation based on optimization: the selection operation of DE is greedy, which can ensure that the excellent solution will not be lost in the evolution process. At the same time, the one-to-one selection competition mechanism can maintain the population diversity better than the sorting selection and bidding selection of the GA algorithm(4)
Variation operation based on difference: this operation enables the variation step size and search direction to be adjusted adaptively according to different objective functions, which is more conducive to global search


## 2.1. Standard Process of Differential Evolution Algorithm

Standard DE is essentially a special genetic algorithm based on real number coding and optimization preserving greedy strategy. Its principle and structure are similar to genetic algorithm. Similarly, starting from the randomly initialized population, genetic operations such as mutation, crossover, and selection are adopted, so as to solve the optimization problem iteratively.Constrained optimization problems are usually expressed as(1)minfxs.t.cx=0,i∈Ecx≥0,i∈I,where E and I are the index sets of equality and inequality constraints, respectively, and cix are constraint functions. The standard process of solving the above general optimization problems by DE is described as below.
(1)
Population initializationAssuming that the population size is NP (i.e., there are NP individuals participating in evolution) and the evolutionary algebraG=0,1,⋯,Gmax, the ith individual in the population is expressed as
(2)Xi,G=x1,i,G,x2,i,G,x3,i,G,⋯,xD,i,G,where D represents the dimension of individual variables. In the initialization stage, in order to make individuals cover the whole search space as much as possible, each individual in the initial population (i.e., G=0) is randomly and evenly generated within the value range Xmin,Xmax specified by the objective function:
(3)xj,i,0=xj,min+randi,j0,1⋅xj,max−xj,min,where xj,max and xj,min are the upper and lower bounds of the jth dimension of individual variable Xi,G, respectively. randi,j0,1 is a random number uniformly distributed in the interval [0,1].
(2)
Mutation operationMutation is the core operation of the DE algorithm. The purpose is to generate new candidate individuals, enhance the diversity of the population, and guide individuals to a favorable search direction. Compared with other evolutionary algorithms, the mutation operation of DE makes full use of population information, which generates mutation vectora based on the difference vector between different individuals. DE/a/b is usually used to represent different mutation strategies, in which a represents the selection mode of basis vector, often including “rand” and “best,” representing the basis vector selected randomly and optimally, respectively; b represents the number of difference vectors. At present, the commonly used mutation strategies areDE/rand/1:
(4)Vi,G=Xr1,G+F⋅Xr2,G−Xr3,G,DE/best/1:
(5)Vi,G=Xbest,G+F⋅Xr2,G−Xr3,G,DE/best/2:
(6)Vi,G=Xbest,G+F⋅Xr1,G−Xr2,G+F⋅Xr3,G−Xr4,G,DE/current−to−best/1:
(7)Vi,G=Xi,G+F⋅Xbest,G−Xi,G+F⋅Xr1,G−Xr2,G,where r1,r2,r3,r4∈1,2,⋯,NP is a random integer that is not equal to i and different from each other, representing the index number of different individuals in the population; Xbest,G is the individual with the best fitness in generation G population; F∈0,2 is the scaling factor, which is used to control the search step size, and its value directly affects the convergence performance of the algorithm.
(3)
CrossoperationThe purpose of crossover is to disturb the individuals in the original population, so as to enhance the exploration of local areas. Binomial crossover is often used in standard DE, and the target vectorXi,G and variation vector Vi,G are discretely crossed in each dimension, resulting in test vectorUi,G=u1,i,G,u2,i,G,⋯,uD,i,G:
(8)uj,i,G=vj,i,G,randj≤CR⋯or⋯j=randni,xj,i,G,randj>CR⋯and⋯j≠randni,where randj∈0,1 is a uniformly distributed random number; randni∈1,2,⋯,D is the randomly selected dimension variable index to ensure that at least one bit of Ui,G is contributed by Vi,G, while for other bits, it is determined by the crossover probability factor CR; CR∈0,1; its value is often related to the specific nature of the problem (such as multimode problem or single-mode problem and independent variable problem or independent variable interdependence problem).
(4)
Select operationBased on the natural selection mechanism of fittest survival, standard DE selects individuals with better fitness to reproduce under the condition of keeping NP unchanged, so as to realize the inheritance of excellent characteristics, and then guide the population to the optimal region and gradually approach the optimal solution. The selection policy can be expressed as(9)Xi,G+1=Ui,G,fUi,G<fXi,G,Xi,G,fUi,G≥fXi,G,where Xi,G+1 is the ith individual of the next-generation population and f is the fitness function. Obviously, the selection operation ensures that the new individual is at least no worse than the original individual and has the function of memorizing the individual optimal solution.
(5)
Algorithm terminationAfter mutation, crossover, and selection, DE completes an update of the population. At this time, it needs to judge whether it meets the operation termination conditions. If satisfied, the current optimal solution is output; otherwise, continue the search. The conditions for determining the termination of the algorithm can be as follows:(a)
The fitness function value of the optimal individual in the two generations of the population is within a set range in successive iterations(b)
Reach the maximum evolutionary algebra set by the algorithmTo sum up, the overall flow of the standard DE algorithm is shown in Figure1.Figure 1
The flowchart of the standard DE algorithm.


## 2.2. Basic Characteristics of Differential Evolution Algorithm

Compared with other evolutionary algorithms, DE has its own uniqueness:(1)
The structure is simple and easy to implement: the control parameters of DE are relatively few, which reduces the difficulty of parameter selection and adjustment in practical application. At the same time, its core difference mutation operator only performs addition and subtraction operations between individual vectors, so it is easy to use(2)
Using real number encoding: compared with binary encoding of GA algorithm, real number encoding is very suitable for dealing with optimization problems in continuous space(3)
Selection operation based on optimization: the selection operation of DE is greedy, which can ensure that the excellent solution will not be lost in the evolution process. At the same time, the one-to-one selection competition mechanism can maintain the population diversity better than the sorting selection and bidding selection of the GA algorithm(4)
Variation operation based on difference: this operation enables the variation step size and search direction to be adjusted adaptively according to different objective functions, which is more conducive to global search


## 3. Probabilistic Analysis of Algorithm Performance in Low-Dimensional Case

### 3.1. Probability Analysis of Algorithm Variation Performance

The feasibility analysis of the optimization algorithm is mainly divided into the analysis of the convergence of the algorithm and the analysis of its search ability. According to the detailed theoretical derivation of DE convergence in literature, the algorithm converges asymptotically with probability 1. Therefore, this paper focuses on the search performance of the algorithm and analyzes two performance indicators: one is the ability to search for the optimal solution in the current local region, and the other is the ability to search other regions out of the current local region. When the population falls into a local region, if the optimal position within the local range can be quickly searched, it means that the algorithm has good local optimization force. If the algorithm can jump out of the local optimal region quickly, it means that the algorithm has strong global optimization force. Clearly, an algorithm that balances the two can be called excellent performance. At present, there are few quantitative theoretical analyses about these two performance indexes in algorithm field. Therefore, based on the principle of probability and statistics, this chapter proposes a visual probability model to quantitatively analyze the behavior ability of the DE algorithm in local areas and discusses the specific influence of algorithm parameters on algorithm performance.The core operation of DE algorithm evolution is the driving differential mutation operation. This chapter first analyzes mutation operationDE/best/1: ViG=XbG+F⋅XjG−XkG.


#### 3.1.1. Analysis Based on Regular Circular Local Area Model

(1) Model Establishment of Local Optimal Neighborhood. Figure 2 shows three kinds of classical test function surfaces and contour lines. It can be seen that local single valley (peak) areas are mostly close to circle or ellipse, which can be simplified as “bowl”-shaped surfaces in modeling.Figure 2
Surfaces and contours of three classical test functions.
(a)
Rastrigin function(b)
Ackley function(c)
Egg holder functionIn this paper, binary space is taken as an example, and similarly, high-dimensional space can be speculated. Figures3(a) and 3 (c) indicate that when the algorithm progresses to a certain stage, all the individuals in the population fall into a local single valley (peak) region in the search space. In order to facilitate modeling and analysis, the neighborhood was idealized into a circle centered on the local extreme point and made with a radius that could include the single valley (peak) region, as shown in Figures 3(b) and 3 (d), where the gray area was a ring domain composed of fitness value contour lines where the contemporary optimal solution and the worst solution were located, so all the solutions in the population were in this ring domain.Figure 3
Idealized model of test function surface within local neighborhood.
(a)(b)(c)(d)The relative ratios ofrb, rw, and rlim in the figure above have a wide variety of possibilities and can be extrapolated from each other. It is neither feasible nor necessary to parse all cases. For the proportion, therefore, this article selects the typical case, evenrb:rw:rlim=0.5:0.8:1, at this point A, B, and C; three areas for the area ratio ofsbetter:swaste:sworse=0.25:0.4:0.35, population’s band from distance, or the area proportion of moderate, away from the local optimum location, and local area border have similar distance to their own band width, facilitate local optimization algorithm, had been studied, and jump out of the performance of the local neighborhood. Thus, the local neighborhood was divided into four parts: region A was the region to be searched for A better solution in the local neighborhood, region B was the ring region where the current population was located, region C was the region with worse fitness than the current population in the neighborhood, and region D was other regions outside the neighborhood.(2) Analysis of Probability Density Distribution of Difference VectorsXjG−XkG. Convention of concept name is shown in Figure 4.Figure 4
Annulus probability distribution of differential vector.The optimal sitting of the local neighborhood is labeled asXlocalmin.The coordinate system is established by taking contemporary random solution vectorXkG as the center coordinate and the line direction between this point and center Xlocalmin as the axis.The distance between the contemporary random solution vectorXkG and the optimal coordinate vector Xlocalmin in the local domain is denoted as rkG.The distance between contemporary optimal solution vectorXbG and optimal coordinate vector Xlocalmin in the local domain is denoted as rbG.The distance between the contemporary worst solution vectorXwG and the optimal coordinate vector Xlocalmin in the local domain is denoted as rwG.The difference vector modulus of contemporary random solution vectorsXjG and XkG is rG,r∈0,rwG+rkG.
(a)
Calculate the probability density of modulusrG of difference vector XjG−XkGThe probability density ofrG is the ratio of the line area srdr (line width micro element is dr) across the belt region to the band area So=π⋅rwG2−π⋅rbG2 with the subtraction vector XjG as the center of the circle and the modulus value of vector XjG−XkG as the radius rG, namely,
(10)Pr=srdrπ⋅rwG2−π⋅rbG2,where the expression of line length sr is
(11)sr=∫αβφ′2t+ψ′2tdt=∫αβdxdt2+dydt2dt.Based on the geometric properties of hollow bands during integration, the analytic analysis of probability density function is discussed in the following four stages:(1)
WhenrG∈0,rkG−rbG, as shown in Figure 5(a), the curve traversed by the difference vector at this stage is a closed circle:(12)PrG=srSo=2⋅∫0πdx/dθ2+dy/dθ2dθSo=2rGSo(2)
WhenrG∈rkG−rbG,rwG−rkG, as shown in Figure 5(b), part of the curve traversed by the difference vector at this stage is subtracted from the inner circle:(13)PrG=srSo=2⋅π−arccosrG2+rkG2−rbG2/2rGrkG⋅rGSo(3)
WhenrG∈rwG−rkG,rkG+rbG, as shown in Figure 6(a), part of the curve traversed by the difference vector at this stage is subtracted by both the inner circle and the outer circle:(14)PrG=2⋅π−arccosrG2+rkG2−rbG2/2rGrkG−arccosrG2+rkG2−rwG2/2rGrkG⋅rGSo(4)
WhenrG∈rkG+rbG,rkG+rwG, as shown in Figure 6(b), part of the curve traversed by the difference vector at this stage is subtracted from the outer circle:(15)PrG=srSo=2⋅arccosrkG2+rG2−rwG2/2rkGrG⋅rGSoFigure 5
Calculation schematic (1) ofrG probability density.
(a)(b)Figure 6
Calculation schematic (2) ofrG probability density.
(a)(b)According to the above analytical formula, the probability density distribution curve ofW can be obtained by numerical integration method, as shown in Figure 7. In the figure, the significant depression in the middle region of the curve is caused by the distribution characteristics of the ring region of the difference vector (the dotted line is the curve characteristic when the inner ring region is not affected). When the outer ring of the population is relatively unchanged, the size of the inner ring (the population is relatively far from the local optimum in the neighborhood) is positively correlated with the width of the middle depression interval.
(b)
Calculate the probability density of angleαG of difference vector XjG−XkG(1)
Whenα∈0,arccosrbG/rkG, as shown in Figure 8(a), the vector does not pass through the inner region of the ring; then,
(16)PαG=sαSo=rkGsinαG+rwG2−rkGcosαG2dsSo(2)
Whenα∈arccosrbG/rkG,π−arccosrbG/rkG, as shown in Figure 8(b), the vector passes through the inner region of the ring; then,
(17)PαG=sαSo=rkGsinαG+rwG2−rkGcosαG2dsSo(3)
Whenα∈π−arccosrbG/rkG,π, symmetric with case (1), the vector does not pass through the inner region of the ring: then,
(18)PαG=rkGsinαG+rwG2−rkGcosαG2−2⋅rbG2−rkGcosαG2dsSo(4)
Whenα∈π,2π, as shown in Figure 9, the vector direction deviates from the inner region of the ring; then,
(19)PαG=sαSo=rwG2−rkGcosαG2−rkGsinαGdsSoFigure 7
Probability density distribution curve of differential vector modulus.Figure 8
Calculation schematic (1) ofαG probability density.
(a)(b)Figure 9
Calculation schematic (2) ofαG probability density.According to the above analytical formula, the probability density distribution curve ofαG is obtained by numerical integration method, as shown in Figure 10.Figure 10
Probability density distribution curve of differential vector angle.It can be seen that, similar to the probability density distribution ofrG, αG is also affected by the size of the inner ring and is reduced in some regions. It can be seen that the modulus value and direction of difference vectorXjG−XkGhave no equal probability distribution, and they are highly coupled with each other. The probability distribution of the vector synthesized with the basis vector will be quite complex, which needs to be discussed in detail.(3) Synthetic Analysis of Difference VectorXjG−XkG and Basis Vector XbG. The combination of basis vector XbG and difference vector XjG−XkG belongs to vector superposition. The process of vector superposition is shown in Figure 11. XbG is randomly distributed on the inner circle with uniform probability, and XkG is randomly distributed on the circle with radius rkG. In the mutation operation, vector XbG and XjG−XkG are superimposed at any included angle β, which is evenly distributed in the interval of 0°,360°. Therefore, as shown in Figure 12, coordinate XkG in the probability distribution area of difference vector XjG−XkG and coordinate XbG in the coordinate system where basis vector XbG reside should be overlapped and rotated once along this point to calculate the probability distribution of the variation vector. The probability sum of the four regions A, B, C, and D is the probability of the variation vector moving into each region.Figure 11
Geometric process of mutation operation (withF=1 as example).Figure 12
All possibilities of superposition of base vector and differential vector.In view of the complexity of the area integral of the overall rotation process in the distribution region of difference vectors, it is neither practical nor necessary to integrate all the distribution probabilities of variation vectors at all included angles. Then, 8 typical samples can be taken at equal intervals (every 45°) in the rotation process, as shown in Figure13. Integrate the probability area occupied by each region of A, B, C, and D in each sample, and then, perform equal-weight averaging to obtain the required results. In the figure, green is area A (better area), and its relative area is ΔS↑. Yellow is zone B (wandering zone) with a relative area of ΔS↕. Red is region C (degraded region), with a relative area of ΔS↓. Blue is area D (jump out area); its relative area is ΔS↗.Figure 13
Equal interval sampling of superposition of base vector and differential vector.Figure13 takes F=1.0 as an example to show the variation of probability distribution in each region when β is from 0 to 360°. In this chapter, the influence of scaling factor F on algorithm search performance is studied. The probability distributions of F=0.25, 0.5, 0.625, 0.75, 1.0, 1.5, and 2.0 are shown in Figure 14 (β=0° and β=45° are taken as examples). It can be seen from several literatures that the classical value of F is between 0.5 and 1.0 and is mostly in the range of 0.5 to 0.7. Therefore, 0.625 is added between 0.5 and 0.75 in this chapter to better study the influence of F on variation performance.Figure 14
Four areal situation for differentβ from F=0.25 to 2.0.
(a)
β=0°(b)
β=45°It can be seen from Figure13 that the areas of 180°~360° and 0°~180° are symmetric, so only the unilateral case of 0°~180° needs to be calculated. Figure 14 shows the general situation of any position in rkG∈rbG,rwG, while different values of rkG may cause significant differences in the probability distribution of variation vectors falling into different regions, as shown in Figure 15.Figure 15
Distribution changes of each area with differentrkG (with F=1 as example)
(a)
rkG=rbG(b)
rkG=rwGTherefore, two extreme cases need to be calculated for each combination ofF and β, as shown in Figure 16.Figure 16
Distribution of two extremes with all kinds ofF andβ.By integrating the probability distribution of each region for the two extreme cases under eachF, six probability distribution curves can be obtained as shown in Figure 17, where max is rkG=rwG and min is rkG=rbG. Both zone B and zone C belong to invalid operation zone and are combined into one consideration in probability distribution.Figure 17
Probability density distribution of different areas which the mutation vectors fall into.By taking the intermediate value ofmax and min, the characteristic curve of the probability distribution of variation vectors falling into each region as F changes as shown in Figure 18 can be obtained.Figure 18
Synthesis of probability density distribution.As can be seen from the figure,F is the best at the junction of blue and green lines (0.7~0.75) from a single solution, and both local optimization and local escape can be considered at this time. Since the DE algorithm is based on population mutation operation, the probability distribution of population operation results needs to be calculated. The probability of a population searching for a local better region after mutation operation and the probability of mutation jumping out of a local region are shown in equations (20) and (21). According to formula (20), the probability P↑NP of the population searching for a locally better region after mutation operation can be calculated when NP and F are combined, as shown in Table 1.
(20)P↑NP=1−1−P↑1NP,(21)P↗NP=1−1−P↗1NP.Table 1
ProbabilityP↑ of searching to better local area with different combinations of NP and F.
FNP=10NP=20NP=30NP=40NP=500.250.99800.9999961110.50.98780.9998510.9999980.99999810.6250.97970.9995880.9999920.99999210.750.94990.9974860.9998740.99987411.00.83400.9724470.9954260.9954260.9998741.50.64280.8723830.9544110.9544110.9941822.00.50710.7570130.8802230.8802230.970896FNP=60NP=70NP=80NP=90NP=1000.25111110.5111110.625111110.75111111.00.999870.9999970.999999111.50.994180.9992580.9997350.9999112.00.970890.9929280.9965140.998281According to the data in Table1, the probability surface composed of P↑ can be drawn, as shown in Figure 19.Figure 19
Probability distribution of combinations of NP andF effect on local optimizing performance (top view and side view).The parameter combination of NP andF can better guarantee the local optimization performance of algorithm mutation operation in this area. Similarly, according to formula (21), the probability P↗ of the population jumping out of local neighborhood after mutation operation can be calculated for various combinations of NP and F, as shown in Table 2.Table 2
ProbabilityP↗ of escaping from local neighborhood with different combinations of NP and F.
FNP=10NP=20NP=30NP=40NP=500.25000000.50.668570.890150.963590.9635930.9960010.6250.906610.991280.999190.9991850.9999930.750.965810.998830.999960.9999611.00.997190.999991111.50.9999911112.011111FNP=60NP=70NP=80NP=90NP=1000.25000000.50.996000.999560.9998540.9999510.6250.9999911110.75111111.0111111.5111112.011111According to the data in Table2, the probability surface composed of P↗ can be drawn, as shown in Figure 20. The combination of NP and F can better guarantee the local escape ability of the mutation operation of the algorithm in this region.Figure 20
Probability distribution of combinations of NP andF effect on local escaping performance.


#### 3.1.2. Analysis Based on Noncircular Local Region Model

In the first section, the probability analysis of variation performance of the algorithm in regular circular local area is carried out by numerical analysis. However, to simplify the local area into a regular circle, the prerequisite condition is too strong for the application background of the algorithm. It can be seen from Figure21 that most of the ellipse neighborhoods with nonzero eccentricity exist in the surface, so this section carries out probability analysis on the local regional variation performance of ellipses under typical parameters following the thought in the previous section.Figure 21
Elliptic valley model with eccentricity being 0.5.The elliptic valley region can be obtained by centrifugal deformation of the circular local single-valley region model in Figure3, as shown in Figure 21, that is, the elliptic model with an eccentricity of 0.5.When the valley region is ellipse (eccentricity is not 0), the vector representation of variation operation of the DE algorithm is shown in Figure22. It can be seen that the coordinate positions of contour lines in this time domain no longer have equal geometric properties, and the possibility of geometric expression process of differential variation operation is more complex than that of circular operation. The optimal solution can be divided into three different situations: the minimum curvature position, the maximum curvature position, and the middle curvature position of the ellipse. The basis points of difference vectors are also divided into the minimum curvature position, the maximum curvature position, and some curvature position in the middle of the ellipse. Figure 23 shows nine typical geometric collocations of mutation operations.Figure 22
Geometric process of mutation operation with eccentricity being 0.5 (withF=1 as example).Figure 23
Nine typical geometry matches.According to the geometric collocation scheme shown in Figure23, numerical integration and analysis are carried out on the important probabilities (the probability of seeking local better and jumping out of local region) of two representative algorithm performances under three typical cases F=0.5, F=1.0, and F=1.5, respectively. Take F=1.0 as an example. The green and blue areas in Figure 24 are shown.Figure 24
Integral areas withF=1.0.According to the integration region of differentF, the curves of two probability values changing with the change of F can be obtained by integration, as shown in Figure 25.Figure 25
Probability density distribution of different areas which the mutation vectors fall into.Comparing Figure25 with Figure 17, it can be seen that the probability distribution of variation vectors falling into “A” and “D” regions almost does not change with the change of eccentricity. In other words, local optimization performance and local jump performance are mainly related to scaling factor and population size when solvers perform differential variation operation in such single-valley/peak region.Mutation operation is the core of the whole differential evolution operation, which determines the main performance of the DE algorithm. The subsequent crossover operation is a further probability optimization operator consistent with the DE algorithm.


### 3.2. Probability Analysis of Algorithm Crossover Performance

After mutation operation, the DE algorithm performs crossover operation to get the test vector. The essence of crossover operation is to selectively merge the vectors before and after mutation in all dimensions, so that the test vector after operation not only has the jump of mutation but also partly retains the traces before mutation; that is, the mutation operation is carried out by multidimensional carrier processing. Obviously, whenCR=1, crossover operation actually does not work, and the crossover result is the variation vector before crossover. When CR=0, crossover operation maximizes to play a conservative role, and only one dimension in crossover result comes from the variation vector, while the other dimensions are the same as the original vector before mutation. When CR value is between 0 and 1, the crossover results are abundant, which requires quantitative discussion.As shown in Figure26, coordinate points UixG and UiyG corresponding to various possible test vectors can be found according to the coordinates of crossvector XiG and variation vector ViG=XbG+XjG−XkG. Figure 26 shows the Cartesian coordinate system of three different angles, and three groups of different test vectors are obtained. Since the direction of the mutation operation has the property of isoprobability in each direction, the included angle between the coordinate axis and the mutation direction of the rectangular coordinate system has equal probability distribution between 0° and 360°. At the same time, all test vectors are evenly distributed on the circular line with the difference vector as diameter because the triangle surrounded by the coordinate point UixG (or UiyG), XiG, and ViG of the test vector must be a right triangle. The actual variation step is modXbG+XjG−XkG−XbG⋅sinγ, where γ is the included angle between the variation vector and the Y-axis of the coordinate system.Figure 26
“Probability circle” which the solution individual falls into by crossover operation.The crossover operation determines the value of the test vector according to the crossover probability factorCR. The value probabilities are shown in equations (22) and (23), respectively:
(22)pXbG+XjG−XkG=50%⋅CR+50%⋅CR=CR,(23)pUbxG=pUbyG=1−pXbG+XjG−XkG2=1−CR2.In order to establish the probability model of crossover operation, after mutation operation, the local region model was rotated around the geometric center until the mutation vector fell on the horizontal dashed line at the right end of the geometric center of the model. As shown in Figure27, the red dot is the position of ViG. The white point is the representative coordinate of XiG, and the dotted blue line shows where UiG might fall after the crossover operation.Figure 27
Various possibilities of crossover operation.In this chapter, the positions of variation vectors that may be generated under different scaling factorsF and the representative coordinates of XiG are traversed numerically, so as to obtain the probabilities of test vector UiG falling into the local better region and jumping out of the local region, respectively. The results are shown in Figure 28.Figure 28
Probability distributions of falling into local better area (a) and escaping from the local area (b).
(a)(b)According to Figure28(a), when the variation vector is in the locally better region (1#~3#), the probability of the test vector falling into the locally better region is large and increases with the increase of CR. When the variation vector is outside the locally preferred region (4#~14#), the probability of the test vector falling into the locally preferred region is very small and decreases with the increase of CR. According to Figure 28(b), when the variation vector is in the locally better region (1#~3#), the probability of the test vector jumping out of the local region is 0. When the variation vector is located outside the local preferred region (4#~14#), the probability of the test vector jumping out of the local region increases rapidly with the variation vector deviating from the center of the local region and increases with the increase of CR.


### 3.3. Probability Analysis of Comprehensive Performance of the Algorithm

Combined with the quantitative results obtained in the previous two sections, this section carries out a comprehensive probability analysis of algorithm performance and jointly examines the impact of scaling factor and crossover probability factor on two important performance indicators of the algorithm. According to the quantitative analysis of individual variation behavior (see Figure18 for details), the probability of variation vectors falling into different regions can be obtained as shown in Table 3.Table 3
Probability of mutation vectors falling into different areas.
The zoom factorThe probability of a mutant falling into a different regionLocally superior regionLocally nonsuperior regionJump out of local areaF=“0.25”43.22%56.78%0.00%F=“0.50”37.50%56.50%6.00%F=“0.625”34.04%50.32%15.65%F=“0.75”29.92%45.96%24.12%F=“1.00”18.59%43.62%37.79%F=“1.50”10.97%25.80%63.24%F=“2.00”7.65%17.13%75.22%According to the quantitative analysis of individual crossover behavior (see Figure28 for details), the probabilities of test vectors falling into different regions can be obtained as shown in Table 4.Table 4
Situations of crossover operation after mutation into different areas.
The condition after crossover operation of variation vectors falling into different regionsCRThe probability of jumping out of a local area after a crossover operationThe probability of falling into a locally better region after crossover operationThe vector of variation is locally superiorThe vector of variation is in a locally nonsuperior regionThe vector of variation is out of placeThe vector of variation is locally superiorThe vector of variation is in a locally nonsuperior regionThe vector of variation is out of place00.00%10.39%63.57%56.32%20.25%8.59%0.10.00%9.35%67.21%60.69%18.22%7.73%0.20.00%8.31%70.85%65.06%16.20%6.87%0.30.00%7.27%74.50%69.42%14.17%6.01%0.40.00%6.23%78.14%73.79%12.15%5.15%0.50.00%5.20%81.78%78.16%10.12%4.29%0.60.00%4.16%85.43%82.53%8.10%3.43%0.70.00%3.12%89.07%86.90%6.07%2.58%0.80.00%2.08%92.71%91.26%4.05%1.72%0.90.00%1.04%96.36%95.63%2.02%0.86%1.00.00%0.00%100.00%100.00%0.00%0.00%Joint calculation of the probability distribution in Tables3 and 4 can obtain the probability of a single solution individual jumping out of the local area and falling into the local better area under different F and CR, as shown in Figure 29.Figure 29
Probability of single solution individual falling into local better area (a) and escaping from the local area (b).
(a)(b)It can be seen from Figure29(a) that the probability of the solution individual jumping out of the local area is highly positively correlated with the scaling factor F, highly negatively correlated with the crossover factor when F is small, and gradually highly positively correlated with the increase of F.It can be seen from Figure29(b) that the probability of solution individuals falling into local optimization is highly inversely correlated with the scaling factor F, highly positively correlated with the crossover factor when F is small, and gradually highly negatively correlated with the increase of F.It can be seen from the above that ifF and CR are both constant values, in order to give consideration to the two performance of the algorithm, it is appropriate to select the middle value of the value range of F, which can be F∈0.5,1. IfFand CR are set by algorithm for adaptive adjustment, it can be set up in the early stage of the evolution of the largerFand CR to emphasize particularly on the performance of the play to jump out of local enhancement algorithm, namely, initial exploration ability; thus improve the odds of finding the global optimal area; and take smaller F in the late evolution and CR, focusing on the local optimization ability, namely, the improved algorithm later production ability, in order to perform a fine search near the optimal solution. At the same time, it can be found that CR always tends to take a larger value, and CR∈0.7,0.9 can be selected.According to formulas (20) and (21), the population size and the comprehensive characteristics of the two search performance are considered in fusion, as shown in Figure30. The geometric average method is used to synthesize the local optimization and local escape, which is better than the arithmetic average method to show the equilibrium degree of the two performances. It can be seen from Figure 30 that when the population size is small (10 solution individuals), F at about 0.75 is obviously superior to other values. Comparatively speaking, the comprehensive search performance is not very sensitive to the value of CR. When the population size reaches 40, the search performance of the algorithm is significantly improved in both aspects. The data conclusions in this section are derived from the circular local region model. It has been verified that the elliptic local region model can obtain a similar result curve.Figure 30
Probability of solution individual escaping from local area and falling into local better area when the population size is 10 (a) and 40 (b).
(a)(b)


## 3.1. Probability Analysis of Algorithm Variation Performance

The feasibility analysis of the optimization algorithm is mainly divided into the analysis of the convergence of the algorithm and the analysis of its search ability. According to the detailed theoretical derivation of DE convergence in literature, the algorithm converges asymptotically with probability 1. Therefore, this paper focuses on the search performance of the algorithm and analyzes two performance indicators: one is the ability to search for the optimal solution in the current local region, and the other is the ability to search other regions out of the current local region. When the population falls into a local region, if the optimal position within the local range can be quickly searched, it means that the algorithm has good local optimization force. If the algorithm can jump out of the local optimal region quickly, it means that the algorithm has strong global optimization force. Clearly, an algorithm that balances the two can be called excellent performance. At present, there are few quantitative theoretical analyses about these two performance indexes in algorithm field. Therefore, based on the principle of probability and statistics, this chapter proposes a visual probability model to quantitatively analyze the behavior ability of the DE algorithm in local areas and discusses the specific influence of algorithm parameters on algorithm performance.The core operation of DE algorithm evolution is the driving differential mutation operation. This chapter first analyzes mutation operationDE/best/1: ViG=XbG+F⋅XjG−XkG.


### 3.1.1. Analysis Based on Regular Circular Local Area Model

(1) Model Establishment of Local Optimal Neighborhood. Figure 2 shows three kinds of classical test function surfaces and contour lines. It can be seen that local single valley (peak) areas are mostly close to circle or ellipse, which can be simplified as “bowl”-shaped surfaces in modeling.Figure 2
Surfaces and contours of three classical test functions.
(a)
Rastrigin function(b)
Ackley function(c)
Egg holder functionIn this paper, binary space is taken as an example, and similarly, high-dimensional space can be speculated. Figures3(a) and 3 (c) indicate that when the algorithm progresses to a certain stage, all the individuals in the population fall into a local single valley (peak) region in the search space. In order to facilitate modeling and analysis, the neighborhood was idealized into a circle centered on the local extreme point and made with a radius that could include the single valley (peak) region, as shown in Figures 3(b) and 3 (d), where the gray area was a ring domain composed of fitness value contour lines where the contemporary optimal solution and the worst solution were located, so all the solutions in the population were in this ring domain.Figure 3
Idealized model of test function surface within local neighborhood.
(a)(b)(c)(d)The relative ratios ofrb, rw, and rlim in the figure above have a wide variety of possibilities and can be extrapolated from each other. It is neither feasible nor necessary to parse all cases. For the proportion, therefore, this article selects the typical case, evenrb:rw:rlim=0.5:0.8:1, at this point A, B, and C; three areas for the area ratio ofsbetter:swaste:sworse=0.25:0.4:0.35, population’s band from distance, or the area proportion of moderate, away from the local optimum location, and local area border have similar distance to their own band width, facilitate local optimization algorithm, had been studied, and jump out of the performance of the local neighborhood. Thus, the local neighborhood was divided into four parts: region A was the region to be searched for A better solution in the local neighborhood, region B was the ring region where the current population was located, region C was the region with worse fitness than the current population in the neighborhood, and region D was other regions outside the neighborhood.(2) Analysis of Probability Density Distribution of Difference VectorsXjG−XkG. Convention of concept name is shown in Figure 4.Figure 4
Annulus probability distribution of differential vector.The optimal sitting of the local neighborhood is labeled asXlocalmin.The coordinate system is established by taking contemporary random solution vectorXkG as the center coordinate and the line direction between this point and center Xlocalmin as the axis.The distance between the contemporary random solution vectorXkG and the optimal coordinate vector Xlocalmin in the local domain is denoted as rkG.The distance between contemporary optimal solution vectorXbG and optimal coordinate vector Xlocalmin in the local domain is denoted as rbG.The distance between the contemporary worst solution vectorXwG and the optimal coordinate vector Xlocalmin in the local domain is denoted as rwG.The difference vector modulus of contemporary random solution vectorsXjG and XkG is rG,r∈0,rwG+rkG.
(a)
Calculate the probability density of modulusrG of difference vector XjG−XkGThe probability density ofrG is the ratio of the line area srdr (line width micro element is dr) across the belt region to the band area So=π⋅rwG2−π⋅rbG2 with the subtraction vector XjG as the center of the circle and the modulus value of vector XjG−XkG as the radius rG, namely,
(10)Pr=srdrπ⋅rwG2−π⋅rbG2,where the expression of line length sr is
(11)sr=∫αβφ′2t+ψ′2tdt=∫αβdxdt2+dydt2dt.Based on the geometric properties of hollow bands during integration, the analytic analysis of probability density function is discussed in the following four stages:(1)
WhenrG∈0,rkG−rbG, as shown in Figure 5(a), the curve traversed by the difference vector at this stage is a closed circle:(12)PrG=srSo=2⋅∫0πdx/dθ2+dy/dθ2dθSo=2rGSo(2)
WhenrG∈rkG−rbG,rwG−rkG, as shown in Figure 5(b), part of the curve traversed by the difference vector at this stage is subtracted from the inner circle:(13)PrG=srSo=2⋅π−arccosrG2+rkG2−rbG2/2rGrkG⋅rGSo(3)
WhenrG∈rwG−rkG,rkG+rbG, as shown in Figure 6(a), part of the curve traversed by the difference vector at this stage is subtracted by both the inner circle and the outer circle:(14)PrG=2⋅π−arccosrG2+rkG2−rbG2/2rGrkG−arccosrG2+rkG2−rwG2/2rGrkG⋅rGSo(4)
WhenrG∈rkG+rbG,rkG+rwG, as shown in Figure 6(b), part of the curve traversed by the difference vector at this stage is subtracted from the outer circle:(15)PrG=srSo=2⋅arccosrkG2+rG2−rwG2/2rkGrG⋅rGSoFigure 5
Calculation schematic (1) ofrG probability density.
(a)(b)Figure 6
Calculation schematic (2) ofrG probability density.
(a)(b)According to the above analytical formula, the probability density distribution curve ofW can be obtained by numerical integration method, as shown in Figure 7. In the figure, the significant depression in the middle region of the curve is caused by the distribution characteristics of the ring region of the difference vector (the dotted line is the curve characteristic when the inner ring region is not affected). When the outer ring of the population is relatively unchanged, the size of the inner ring (the population is relatively far from the local optimum in the neighborhood) is positively correlated with the width of the middle depression interval.
(b)
Calculate the probability density of angleαG of difference vector XjG−XkG(1)
Whenα∈0,arccosrbG/rkG, as shown in Figure 8(a), the vector does not pass through the inner region of the ring; then,
(16)PαG=sαSo=rkGsinαG+rwG2−rkGcosαG2dsSo(2)
Whenα∈arccosrbG/rkG,π−arccosrbG/rkG, as shown in Figure 8(b), the vector passes through the inner region of the ring; then,
(17)PαG=sαSo=rkGsinαG+rwG2−rkGcosαG2dsSo(3)
Whenα∈π−arccosrbG/rkG,π, symmetric with case (1), the vector does not pass through the inner region of the ring: then,
(18)PαG=rkGsinαG+rwG2−rkGcosαG2−2⋅rbG2−rkGcosαG2dsSo(4)
Whenα∈π,2π, as shown in Figure 9, the vector direction deviates from the inner region of the ring; then,
(19)PαG=sαSo=rwG2−rkGcosαG2−rkGsinαGdsSoFigure 7
Probability density distribution curve of differential vector modulus.Figure 8
Calculation schematic (1) ofαG probability density.
(a)(b)Figure 9
Calculation schematic (2) ofαG probability density.According to the above analytical formula, the probability density distribution curve ofαG is obtained by numerical integration method, as shown in Figure 10.Figure 10
Probability density distribution curve of differential vector angle.It can be seen that, similar to the probability density distribution ofrG, αG is also affected by the size of the inner ring and is reduced in some regions. It can be seen that the modulus value and direction of difference vectorXjG−XkGhave no equal probability distribution, and they are highly coupled with each other. The probability distribution of the vector synthesized with the basis vector will be quite complex, which needs to be discussed in detail.(3) Synthetic Analysis of Difference VectorXjG−XkG and Basis Vector XbG. The combination of basis vector XbG and difference vector XjG−XkG belongs to vector superposition. The process of vector superposition is shown in Figure 11. XbG is randomly distributed on the inner circle with uniform probability, and XkG is randomly distributed on the circle with radius rkG. In the mutation operation, vector XbG and XjG−XkG are superimposed at any included angle β, which is evenly distributed in the interval of 0°,360°. Therefore, as shown in Figure 12, coordinate XkG in the probability distribution area of difference vector XjG−XkG and coordinate XbG in the coordinate system where basis vector XbG reside should be overlapped and rotated once along this point to calculate the probability distribution of the variation vector. The probability sum of the four regions A, B, C, and D is the probability of the variation vector moving into each region.Figure 11
Geometric process of mutation operation (withF=1 as example).Figure 12
All possibilities of superposition of base vector and differential vector.In view of the complexity of the area integral of the overall rotation process in the distribution region of difference vectors, it is neither practical nor necessary to integrate all the distribution probabilities of variation vectors at all included angles. Then, 8 typical samples can be taken at equal intervals (every 45°) in the rotation process, as shown in Figure13. Integrate the probability area occupied by each region of A, B, C, and D in each sample, and then, perform equal-weight averaging to obtain the required results. In the figure, green is area A (better area), and its relative area is ΔS↑. Yellow is zone B (wandering zone) with a relative area of ΔS↕. Red is region C (degraded region), with a relative area of ΔS↓. Blue is area D (jump out area); its relative area is ΔS↗.Figure 13
Equal interval sampling of superposition of base vector and differential vector.Figure13 takes F=1.0 as an example to show the variation of probability distribution in each region when β is from 0 to 360°. In this chapter, the influence of scaling factor F on algorithm search performance is studied. The probability distributions of F=0.25, 0.5, 0.625, 0.75, 1.0, 1.5, and 2.0 are shown in Figure 14 (β=0° and β=45° are taken as examples). It can be seen from several literatures that the classical value of F is between 0.5 and 1.0 and is mostly in the range of 0.5 to 0.7. Therefore, 0.625 is added between 0.5 and 0.75 in this chapter to better study the influence of F on variation performance.Figure 14
Four areal situation for differentβ from F=0.25 to 2.0.
(a)
β=0°(b)
β=45°It can be seen from Figure13 that the areas of 180°~360° and 0°~180° are symmetric, so only the unilateral case of 0°~180° needs to be calculated. Figure 14 shows the general situation of any position in rkG∈rbG,rwG, while different values of rkG may cause significant differences in the probability distribution of variation vectors falling into different regions, as shown in Figure 15.Figure 15
Distribution changes of each area with differentrkG (with F=1 as example)
(a)
rkG=rbG(b)
rkG=rwGTherefore, two extreme cases need to be calculated for each combination ofF and β, as shown in Figure 16.Figure 16
Distribution of two extremes with all kinds ofF andβ.By integrating the probability distribution of each region for the two extreme cases under eachF, six probability distribution curves can be obtained as shown in Figure 17, where max is rkG=rwG and min is rkG=rbG. Both zone B and zone C belong to invalid operation zone and are combined into one consideration in probability distribution.Figure 17
Probability density distribution of different areas which the mutation vectors fall into.By taking the intermediate value ofmax and min, the characteristic curve of the probability distribution of variation vectors falling into each region as F changes as shown in Figure 18 can be obtained.Figure 18
Synthesis of probability density distribution.As can be seen from the figure,F is the best at the junction of blue and green lines (0.7~0.75) from a single solution, and both local optimization and local escape can be considered at this time. Since the DE algorithm is based on population mutation operation, the probability distribution of population operation results needs to be calculated. The probability of a population searching for a local better region after mutation operation and the probability of mutation jumping out of a local region are shown in equations (20) and (21). According to formula (20), the probability P↑NP of the population searching for a locally better region after mutation operation can be calculated when NP and F are combined, as shown in Table 1.
(20)P↑NP=1−1−P↑1NP,(21)P↗NP=1−1−P↗1NP.Table 1
ProbabilityP↑ of searching to better local area with different combinations of NP and F.
FNP=10NP=20NP=30NP=40NP=500.250.99800.9999961110.50.98780.9998510.9999980.99999810.6250.97970.9995880.9999920.99999210.750.94990.9974860.9998740.99987411.00.83400.9724470.9954260.9954260.9998741.50.64280.8723830.9544110.9544110.9941822.00.50710.7570130.8802230.8802230.970896FNP=60NP=70NP=80NP=90NP=1000.25111110.5111110.625111110.75111111.00.999870.9999970.999999111.50.994180.9992580.9997350.9999112.00.970890.9929280.9965140.998281According to the data in Table1, the probability surface composed of P↑ can be drawn, as shown in Figure 19.Figure 19
Probability distribution of combinations of NP andF effect on local optimizing performance (top view and side view).The parameter combination of NP andF can better guarantee the local optimization performance of algorithm mutation operation in this area. Similarly, according to formula (21), the probability P↗ of the population jumping out of local neighborhood after mutation operation can be calculated for various combinations of NP and F, as shown in Table 2.Table 2
ProbabilityP↗ of escaping from local neighborhood with different combinations of NP and F.
FNP=10NP=20NP=30NP=40NP=500.25000000.50.668570.890150.963590.9635930.9960010.6250.906610.991280.999190.9991850.9999930.750.965810.998830.999960.9999611.00.997190.999991111.50.9999911112.011111FNP=60NP=70NP=80NP=90NP=1000.25000000.50.996000.999560.9998540.9999510.6250.9999911110.75111111.0111111.5111112.011111According to the data in Table2, the probability surface composed of P↗ can be drawn, as shown in Figure 20. The combination of NP and F can better guarantee the local escape ability of the mutation operation of the algorithm in this region.Figure 20
Probability distribution of combinations of NP andF effect on local escaping performance.


### 3.1.2. Analysis Based on Noncircular Local Region Model

In the first section, the probability analysis of variation performance of the algorithm in regular circular local area is carried out by numerical analysis. However, to simplify the local area into a regular circle, the prerequisite condition is too strong for the application background of the algorithm. It can be seen from Figure21 that most of the ellipse neighborhoods with nonzero eccentricity exist in the surface, so this section carries out probability analysis on the local regional variation performance of ellipses under typical parameters following the thought in the previous section.Figure 21
Elliptic valley model with eccentricity being 0.5.The elliptic valley region can be obtained by centrifugal deformation of the circular local single-valley region model in Figure3, as shown in Figure 21, that is, the elliptic model with an eccentricity of 0.5.When the valley region is ellipse (eccentricity is not 0), the vector representation of variation operation of the DE algorithm is shown in Figure22. It can be seen that the coordinate positions of contour lines in this time domain no longer have equal geometric properties, and the possibility of geometric expression process of differential variation operation is more complex than that of circular operation. The optimal solution can be divided into three different situations: the minimum curvature position, the maximum curvature position, and the middle curvature position of the ellipse. The basis points of difference vectors are also divided into the minimum curvature position, the maximum curvature position, and some curvature position in the middle of the ellipse. Figure 23 shows nine typical geometric collocations of mutation operations.Figure 22
Geometric process of mutation operation with eccentricity being 0.5 (withF=1 as example).Figure 23
Nine typical geometry matches.According to the geometric collocation scheme shown in Figure23, numerical integration and analysis are carried out on the important probabilities (the probability of seeking local better and jumping out of local region) of two representative algorithm performances under three typical cases F=0.5, F=1.0, and F=1.5, respectively. Take F=1.0 as an example. The green and blue areas in Figure 24 are shown.Figure 24
Integral areas withF=1.0.According to the integration region of differentF, the curves of two probability values changing with the change of F can be obtained by integration, as shown in Figure 25.Figure 25
Probability density distribution of different areas which the mutation vectors fall into.Comparing Figure25 with Figure 17, it can be seen that the probability distribution of variation vectors falling into “A” and “D” regions almost does not change with the change of eccentricity. In other words, local optimization performance and local jump performance are mainly related to scaling factor and population size when solvers perform differential variation operation in such single-valley/peak region.Mutation operation is the core of the whole differential evolution operation, which determines the main performance of the DE algorithm. The subsequent crossover operation is a further probability optimization operator consistent with the DE algorithm.


## 3.1.1. Analysis Based on Regular Circular Local Area Model

(1) Model Establishment of Local Optimal Neighborhood. Figure 2 shows three kinds of classical test function surfaces and contour lines. It can be seen that local single valley (peak) areas are mostly close to circle or ellipse, which can be simplified as “bowl”-shaped surfaces in modeling.Figure 2
Surfaces and contours of three classical test functions.
(a)
Rastrigin function(b)
Ackley function(c)
Egg holder functionIn this paper, binary space is taken as an example, and similarly, high-dimensional space can be speculated. Figures3(a) and 3 (c) indicate that when the algorithm progresses to a certain stage, all the individuals in the population fall into a local single valley (peak) region in the search space. In order to facilitate modeling and analysis, the neighborhood was idealized into a circle centered on the local extreme point and made with a radius that could include the single valley (peak) region, as shown in Figures 3(b) and 3 (d), where the gray area was a ring domain composed of fitness value contour lines where the contemporary optimal solution and the worst solution were located, so all the solutions in the population were in this ring domain.Figure 3
Idealized model of test function surface within local neighborhood.
(a)(b)(c)(d)The relative ratios ofrb, rw, and rlim in the figure above have a wide variety of possibilities and can be extrapolated from each other. It is neither feasible nor necessary to parse all cases. For the proportion, therefore, this article selects the typical case, evenrb:rw:rlim=0.5:0.8:1, at this point A, B, and C; three areas for the area ratio ofsbetter:swaste:sworse=0.25:0.4:0.35, population’s band from distance, or the area proportion of moderate, away from the local optimum location, and local area border have similar distance to their own band width, facilitate local optimization algorithm, had been studied, and jump out of the performance of the local neighborhood. Thus, the local neighborhood was divided into four parts: region A was the region to be searched for A better solution in the local neighborhood, region B was the ring region where the current population was located, region C was the region with worse fitness than the current population in the neighborhood, and region D was other regions outside the neighborhood.(2) Analysis of Probability Density Distribution of Difference VectorsXjG−XkG. Convention of concept name is shown in Figure 4.Figure 4
Annulus probability distribution of differential vector.The optimal sitting of the local neighborhood is labeled asXlocalmin.The coordinate system is established by taking contemporary random solution vectorXkG as the center coordinate and the line direction between this point and center Xlocalmin as the axis.The distance between the contemporary random solution vectorXkG and the optimal coordinate vector Xlocalmin in the local domain is denoted as rkG.The distance between contemporary optimal solution vectorXbG and optimal coordinate vector Xlocalmin in the local domain is denoted as rbG.The distance between the contemporary worst solution vectorXwG and the optimal coordinate vector Xlocalmin in the local domain is denoted as rwG.The difference vector modulus of contemporary random solution vectorsXjG and XkG is rG,r∈0,rwG+rkG.
(a)
Calculate the probability density of modulusrG of difference vector XjG−XkGThe probability density ofrG is the ratio of the line area srdr (line width micro element is dr) across the belt region to the band area So=π⋅rwG2−π⋅rbG2 with the subtraction vector XjG as the center of the circle and the modulus value of vector XjG−XkG as the radius rG, namely,
(10)Pr=srdrπ⋅rwG2−π⋅rbG2,where the expression of line length sr is
(11)sr=∫αβφ′2t+ψ′2tdt=∫αβdxdt2+dydt2dt.Based on the geometric properties of hollow bands during integration, the analytic analysis of probability density function is discussed in the following four stages:(1)
WhenrG∈0,rkG−rbG, as shown in Figure 5(a), the curve traversed by the difference vector at this stage is a closed circle:(12)PrG=srSo=2⋅∫0πdx/dθ2+dy/dθ2dθSo=2rGSo(2)
WhenrG∈rkG−rbG,rwG−rkG, as shown in Figure 5(b), part of the curve traversed by the difference vector at this stage is subtracted from the inner circle:(13)PrG=srSo=2⋅π−arccosrG2+rkG2−rbG2/2rGrkG⋅rGSo(3)
WhenrG∈rwG−rkG,rkG+rbG, as shown in Figure 6(a), part of the curve traversed by the difference vector at this stage is subtracted by both the inner circle and the outer circle:(14)PrG=2⋅π−arccosrG2+rkG2−rbG2/2rGrkG−arccosrG2+rkG2−rwG2/2rGrkG⋅rGSo(4)
WhenrG∈rkG+rbG,rkG+rwG, as shown in Figure 6(b), part of the curve traversed by the difference vector at this stage is subtracted from the outer circle:(15)PrG=srSo=2⋅arccosrkG2+rG2−rwG2/2rkGrG⋅rGSoFigure 5
Calculation schematic (1) ofrG probability density.
(a)(b)Figure 6
Calculation schematic (2) ofrG probability density.
(a)(b)According to the above analytical formula, the probability density distribution curve ofW can be obtained by numerical integration method, as shown in Figure 7. In the figure, the significant depression in the middle region of the curve is caused by the distribution characteristics of the ring region of the difference vector (the dotted line is the curve characteristic when the inner ring region is not affected). When the outer ring of the population is relatively unchanged, the size of the inner ring (the population is relatively far from the local optimum in the neighborhood) is positively correlated with the width of the middle depression interval.
(b)
Calculate the probability density of angleαG of difference vector XjG−XkG(1)
Whenα∈0,arccosrbG/rkG, as shown in Figure 8(a), the vector does not pass through the inner region of the ring; then,
(16)PαG=sαSo=rkGsinαG+rwG2−rkGcosαG2dsSo(2)
Whenα∈arccosrbG/rkG,π−arccosrbG/rkG, as shown in Figure 8(b), the vector passes through the inner region of the ring; then,
(17)PαG=sαSo=rkGsinαG+rwG2−rkGcosαG2dsSo(3)
Whenα∈π−arccosrbG/rkG,π, symmetric with case (1), the vector does not pass through the inner region of the ring: then,
(18)PαG=rkGsinαG+rwG2−rkGcosαG2−2⋅rbG2−rkGcosαG2dsSo(4)
Whenα∈π,2π, as shown in Figure 9, the vector direction deviates from the inner region of the ring; then,
(19)PαG=sαSo=rwG2−rkGcosαG2−rkGsinαGdsSoFigure 7
Probability density distribution curve of differential vector modulus.Figure 8
Calculation schematic (1) ofαG probability density.
(a)(b)Figure 9
Calculation schematic (2) ofαG probability density.According to the above analytical formula, the probability density distribution curve ofαG is obtained by numerical integration method, as shown in Figure 10.Figure 10
Probability density distribution curve of differential vector angle.It can be seen that, similar to the probability density distribution ofrG, αG is also affected by the size of the inner ring and is reduced in some regions. It can be seen that the modulus value and direction of difference vectorXjG−XkGhave no equal probability distribution, and they are highly coupled with each other. The probability distribution of the vector synthesized with the basis vector will be quite complex, which needs to be discussed in detail.(3) Synthetic Analysis of Difference VectorXjG−XkG and Basis Vector XbG. The combination of basis vector XbG and difference vector XjG−XkG belongs to vector superposition. The process of vector superposition is shown in Figure 11. XbG is randomly distributed on the inner circle with uniform probability, and XkG is randomly distributed on the circle with radius rkG. In the mutation operation, vector XbG and XjG−XkG are superimposed at any included angle β, which is evenly distributed in the interval of 0°,360°. Therefore, as shown in Figure 12, coordinate XkG in the probability distribution area of difference vector XjG−XkG and coordinate XbG in the coordinate system where basis vector XbG reside should be overlapped and rotated once along this point to calculate the probability distribution of the variation vector. The probability sum of the four regions A, B, C, and D is the probability of the variation vector moving into each region.Figure 11
Geometric process of mutation operation (withF=1 as example).Figure 12
All possibilities of superposition of base vector and differential vector.In view of the complexity of the area integral of the overall rotation process in the distribution region of difference vectors, it is neither practical nor necessary to integrate all the distribution probabilities of variation vectors at all included angles. Then, 8 typical samples can be taken at equal intervals (every 45°) in the rotation process, as shown in Figure13. Integrate the probability area occupied by each region of A, B, C, and D in each sample, and then, perform equal-weight averaging to obtain the required results. In the figure, green is area A (better area), and its relative area is ΔS↑. Yellow is zone B (wandering zone) with a relative area of ΔS↕. Red is region C (degraded region), with a relative area of ΔS↓. Blue is area D (jump out area); its relative area is ΔS↗.Figure 13
Equal interval sampling of superposition of base vector and differential vector.Figure13 takes F=1.0 as an example to show the variation of probability distribution in each region when β is from 0 to 360°. In this chapter, the influence of scaling factor F on algorithm search performance is studied. The probability distributions of F=0.25, 0.5, 0.625, 0.75, 1.0, 1.5, and 2.0 are shown in Figure 14 (β=0° and β=45° are taken as examples). It can be seen from several literatures that the classical value of F is between 0.5 and 1.0 and is mostly in the range of 0.5 to 0.7. Therefore, 0.625 is added between 0.5 and 0.75 in this chapter to better study the influence of F on variation performance.Figure 14
Four areal situation for differentβ from F=0.25 to 2.0.
(a)
β=0°(b)
β=45°It can be seen from Figure13 that the areas of 180°~360° and 0°~180° are symmetric, so only the unilateral case of 0°~180° needs to be calculated. Figure 14 shows the general situation of any position in rkG∈rbG,rwG, while different values of rkG may cause significant differences in the probability distribution of variation vectors falling into different regions, as shown in Figure 15.Figure 15
Distribution changes of each area with differentrkG (with F=1 as example)
(a)
rkG=rbG(b)
rkG=rwGTherefore, two extreme cases need to be calculated for each combination ofF and β, as shown in Figure 16.Figure 16
Distribution of two extremes with all kinds ofF andβ.By integrating the probability distribution of each region for the two extreme cases under eachF, six probability distribution curves can be obtained as shown in Figure 17, where max is rkG=rwG and min is rkG=rbG. Both zone B and zone C belong to invalid operation zone and are combined into one consideration in probability distribution.Figure 17
Probability density distribution of different areas which the mutation vectors fall into.By taking the intermediate value ofmax and min, the characteristic curve of the probability distribution of variation vectors falling into each region as F changes as shown in Figure 18 can be obtained.Figure 18
Synthesis of probability density distribution.As can be seen from the figure,F is the best at the junction of blue and green lines (0.7~0.75) from a single solution, and both local optimization and local escape can be considered at this time. Since the DE algorithm is based on population mutation operation, the probability distribution of population operation results needs to be calculated. The probability of a population searching for a local better region after mutation operation and the probability of mutation jumping out of a local region are shown in equations (20) and (21). According to formula (20), the probability P↑NP of the population searching for a locally better region after mutation operation can be calculated when NP and F are combined, as shown in Table 1.
(20)P↑NP=1−1−P↑1NP,(21)P↗NP=1−1−P↗1NP.Table 1
ProbabilityP↑ of searching to better local area with different combinations of NP and F.
FNP=10NP=20NP=30NP=40NP=500.250.99800.9999961110.50.98780.9998510.9999980.99999810.6250.97970.9995880.9999920.99999210.750.94990.9974860.9998740.99987411.00.83400.9724470.9954260.9954260.9998741.50.64280.8723830.9544110.9544110.9941822.00.50710.7570130.8802230.8802230.970896FNP=60NP=70NP=80NP=90NP=1000.25111110.5111110.625111110.75111111.00.999870.9999970.999999111.50.994180.9992580.9997350.9999112.00.970890.9929280.9965140.998281According to the data in Table1, the probability surface composed of P↑ can be drawn, as shown in Figure 19.Figure 19
Probability distribution of combinations of NP andF effect on local optimizing performance (top view and side view).The parameter combination of NP andF can better guarantee the local optimization performance of algorithm mutation operation in this area. Similarly, according to formula (21), the probability P↗ of the population jumping out of local neighborhood after mutation operation can be calculated for various combinations of NP and F, as shown in Table 2.Table 2
ProbabilityP↗ of escaping from local neighborhood with different combinations of NP and F.
FNP=10NP=20NP=30NP=40NP=500.25000000.50.668570.890150.963590.9635930.9960010.6250.906610.991280.999190.9991850.9999930.750.965810.998830.999960.9999611.00.997190.999991111.50.9999911112.011111FNP=60NP=70NP=80NP=90NP=1000.25000000.50.996000.999560.9998540.9999510.6250.9999911110.75111111.0111111.5111112.011111According to the data in Table2, the probability surface composed of P↗ can be drawn, as shown in Figure 20. The combination of NP and F can better guarantee the local escape ability of the mutation operation of the algorithm in this region.Figure 20
Probability distribution of combinations of NP andF effect on local escaping performance.


## 3.1.2. Analysis Based on Noncircular Local Region Model

In the first section, the probability analysis of variation performance of the algorithm in regular circular local area is carried out by numerical analysis. However, to simplify the local area into a regular circle, the prerequisite condition is too strong for the application background of the algorithm. It can be seen from Figure21 that most of the ellipse neighborhoods with nonzero eccentricity exist in the surface, so this section carries out probability analysis on the local regional variation performance of ellipses under typical parameters following the thought in the previous section.Figure 21
Elliptic valley model with eccentricity being 0.5.The elliptic valley region can be obtained by centrifugal deformation of the circular local single-valley region model in Figure3, as shown in Figure 21, that is, the elliptic model with an eccentricity of 0.5.When the valley region is ellipse (eccentricity is not 0), the vector representation of variation operation of the DE algorithm is shown in Figure22. It can be seen that the coordinate positions of contour lines in this time domain no longer have equal geometric properties, and the possibility of geometric expression process of differential variation operation is more complex than that of circular operation. The optimal solution can be divided into three different situations: the minimum curvature position, the maximum curvature position, and the middle curvature position of the ellipse. The basis points of difference vectors are also divided into the minimum curvature position, the maximum curvature position, and some curvature position in the middle of the ellipse. Figure 23 shows nine typical geometric collocations of mutation operations.Figure 22
Geometric process of mutation operation with eccentricity being 0.5 (withF=1 as example).Figure 23
Nine typical geometry matches.According to the geometric collocation scheme shown in Figure23, numerical integration and analysis are carried out on the important probabilities (the probability of seeking local better and jumping out of local region) of two representative algorithm performances under three typical cases F=0.5, F=1.0, and F=1.5, respectively. Take F=1.0 as an example. The green and blue areas in Figure 24 are shown.Figure 24
Integral areas withF=1.0.According to the integration region of differentF, the curves of two probability values changing with the change of F can be obtained by integration, as shown in Figure 25.Figure 25
Probability density distribution of different areas which the mutation vectors fall into.Comparing Figure25 with Figure 17, it can be seen that the probability distribution of variation vectors falling into “A” and “D” regions almost does not change with the change of eccentricity. In other words, local optimization performance and local jump performance are mainly related to scaling factor and population size when solvers perform differential variation operation in such single-valley/peak region.Mutation operation is the core of the whole differential evolution operation, which determines the main performance of the DE algorithm. The subsequent crossover operation is a further probability optimization operator consistent with the DE algorithm.


## 3.2. Probability Analysis of Algorithm Crossover Performance

After mutation operation, the DE algorithm performs crossover operation to get the test vector. The essence of crossover operation is to selectively merge the vectors before and after mutation in all dimensions, so that the test vector after operation not only has the jump of mutation but also partly retains the traces before mutation; that is, the mutation operation is carried out by multidimensional carrier processing. Obviously, whenCR=1, crossover operation actually does not work, and the crossover result is the variation vector before crossover. When CR=0, crossover operation maximizes to play a conservative role, and only one dimension in crossover result comes from the variation vector, while the other dimensions are the same as the original vector before mutation. When CR value is between 0 and 1, the crossover results are abundant, which requires quantitative discussion.As shown in Figure26, coordinate points UixG and UiyG corresponding to various possible test vectors can be found according to the coordinates of crossvector XiG and variation vector ViG=XbG+XjG−XkG. Figure 26 shows the Cartesian coordinate system of three different angles, and three groups of different test vectors are obtained. Since the direction of the mutation operation has the property of isoprobability in each direction, the included angle between the coordinate axis and the mutation direction of the rectangular coordinate system has equal probability distribution between 0° and 360°. At the same time, all test vectors are evenly distributed on the circular line with the difference vector as diameter because the triangle surrounded by the coordinate point UixG (or UiyG), XiG, and ViG of the test vector must be a right triangle. The actual variation step is modXbG+XjG−XkG−XbG⋅sinγ, where γ is the included angle between the variation vector and the Y-axis of the coordinate system.Figure 26
“Probability circle” which the solution individual falls into by crossover operation.The crossover operation determines the value of the test vector according to the crossover probability factorCR. The value probabilities are shown in equations (22) and (23), respectively:
(22)pXbG+XjG−XkG=50%⋅CR+50%⋅CR=CR,(23)pUbxG=pUbyG=1−pXbG+XjG−XkG2=1−CR2.In order to establish the probability model of crossover operation, after mutation operation, the local region model was rotated around the geometric center until the mutation vector fell on the horizontal dashed line at the right end of the geometric center of the model. As shown in Figure27, the red dot is the position of ViG. The white point is the representative coordinate of XiG, and the dotted blue line shows where UiG might fall after the crossover operation.Figure 27
Various possibilities of crossover operation.In this chapter, the positions of variation vectors that may be generated under different scaling factorsF and the representative coordinates of XiG are traversed numerically, so as to obtain the probabilities of test vector UiG falling into the local better region and jumping out of the local region, respectively. The results are shown in Figure 28.Figure 28
Probability distributions of falling into local better area (a) and escaping from the local area (b).
(a)(b)According to Figure28(a), when the variation vector is in the locally better region (1#~3#), the probability of the test vector falling into the locally better region is large and increases with the increase of CR. When the variation vector is outside the locally preferred region (4#~14#), the probability of the test vector falling into the locally preferred region is very small and decreases with the increase of CR. According to Figure 28(b), when the variation vector is in the locally better region (1#~3#), the probability of the test vector jumping out of the local region is 0. When the variation vector is located outside the local preferred region (4#~14#), the probability of the test vector jumping out of the local region increases rapidly with the variation vector deviating from the center of the local region and increases with the increase of CR.


## 3.3. Probability Analysis of Comprehensive Performance of the Algorithm

Combined with the quantitative results obtained in the previous two sections, this section carries out a comprehensive probability analysis of algorithm performance and jointly examines the impact of scaling factor and crossover probability factor on two important performance indicators of the algorithm. According to the quantitative analysis of individual variation behavior (see Figure18 for details), the probability of variation vectors falling into different regions can be obtained as shown in Table 3.Table 3
Probability of mutation vectors falling into different areas.
The zoom factorThe probability of a mutant falling into a different regionLocally superior regionLocally nonsuperior regionJump out of local areaF=“0.25”43.22%56.78%0.00%F=“0.50”37.50%56.50%6.00%F=“0.625”34.04%50.32%15.65%F=“0.75”29.92%45.96%24.12%F=“1.00”18.59%43.62%37.79%F=“1.50”10.97%25.80%63.24%F=“2.00”7.65%17.13%75.22%According to the quantitative analysis of individual crossover behavior (see Figure28 for details), the probabilities of test vectors falling into different regions can be obtained as shown in Table 4.Table 4
Situations of crossover operation after mutation into different areas.
The condition after crossover operation of variation vectors falling into different regionsCRThe probability of jumping out of a local area after a crossover operationThe probability of falling into a locally better region after crossover operationThe vector of variation is locally superiorThe vector of variation is in a locally nonsuperior regionThe vector of variation is out of placeThe vector of variation is locally superiorThe vector of variation is in a locally nonsuperior regionThe vector of variation is out of place00.00%10.39%63.57%56.32%20.25%8.59%0.10.00%9.35%67.21%60.69%18.22%7.73%0.20.00%8.31%70.85%65.06%16.20%6.87%0.30.00%7.27%74.50%69.42%14.17%6.01%0.40.00%6.23%78.14%73.79%12.15%5.15%0.50.00%5.20%81.78%78.16%10.12%4.29%0.60.00%4.16%85.43%82.53%8.10%3.43%0.70.00%3.12%89.07%86.90%6.07%2.58%0.80.00%2.08%92.71%91.26%4.05%1.72%0.90.00%1.04%96.36%95.63%2.02%0.86%1.00.00%0.00%100.00%100.00%0.00%0.00%Joint calculation of the probability distribution in Tables3 and 4 can obtain the probability of a single solution individual jumping out of the local area and falling into the local better area under different F and CR, as shown in Figure 29.Figure 29
Probability of single solution individual falling into local better area (a) and escaping from the local area (b).
(a)(b)It can be seen from Figure29(a) that the probability of the solution individual jumping out of the local area is highly positively correlated with the scaling factor F, highly negatively correlated with the crossover factor when F is small, and gradually highly positively correlated with the increase of F.It can be seen from Figure29(b) that the probability of solution individuals falling into local optimization is highly inversely correlated with the scaling factor F, highly positively correlated with the crossover factor when F is small, and gradually highly negatively correlated with the increase of F.It can be seen from the above that ifF and CR are both constant values, in order to give consideration to the two performance of the algorithm, it is appropriate to select the middle value of the value range of F, which can be F∈0.5,1. IfFand CR are set by algorithm for adaptive adjustment, it can be set up in the early stage of the evolution of the largerFand CR to emphasize particularly on the performance of the play to jump out of local enhancement algorithm, namely, initial exploration ability; thus improve the odds of finding the global optimal area; and take smaller F in the late evolution and CR, focusing on the local optimization ability, namely, the improved algorithm later production ability, in order to perform a fine search near the optimal solution. At the same time, it can be found that CR always tends to take a larger value, and CR∈0.7,0.9 can be selected.According to formulas (20) and (21), the population size and the comprehensive characteristics of the two search performance are considered in fusion, as shown in Figure30. The geometric average method is used to synthesize the local optimization and local escape, which is better than the arithmetic average method to show the equilibrium degree of the two performances. It can be seen from Figure 30 that when the population size is small (10 solution individuals), F at about 0.75 is obviously superior to other values. Comparatively speaking, the comprehensive search performance is not very sensitive to the value of CR. When the population size reaches 40, the search performance of the algorithm is significantly improved in both aspects. The data conclusions in this section are derived from the circular local region model. It has been verified that the elliptic local region model can obtain a similar result curve.Figure 30
Probability of solution individual escaping from local area and falling into local better area when the population size is 10 (a) and 40 (b).
(a)(b)


## 4. Conclusion

Not only the DE algorithm, the performance of almost all optimization algorithms is mainly concerned with the search speed of global optimal solutions and the ability to get rid of local optimal solutions. The simplest method is to observe the iterative performance of all populations after they fall into a local neighborhood in the search space. On the one hand, observe how quickly individuals of the populations can find the optimal solution in this neighborhood through cooperation. On the other hand, observe how likely the population can jump out of the local neighborhood and continue to test other regions of the solution space.Based on the above situation, this paper puts forward for the first time the probability analysis method of DE algorithm mutation and crossover operation transient process, established a visual geometric probability model, quantitatively analyzed the behavior ability of DE algorithm in local areas, and obtained the characteristic curve of the impact of population size and scaling factor on algorithm mutation performance and the quantitative curve of the influence of the probability crossover factor on the crossover performance. Finally, the influence of parameters on the comprehensive performance of the algorithm is investigated, and a certain rule of parameter setting is obtained. In general, this chapter provides a theoretical reference for the application of the algorithm through the establishment and analysis of the algorithm probability model.This analysis model is also applicable to the performance analysis of other intelligent optimization algorithms. In the future, the study on how to analyze the solution performance of the algorithm in high-dimensional solution space through this model will be further studied.



---

*Source: 3936999-2022-06-03.xml*