An improved immune algorithm with parallel mutation and its application

1.   Introduction

Due to their potential to solve complex optimal problems, nature-inspired algorithms (NIAs), such as evolutionary algorithms, artificial immune systems, swarm intelligence, and simulated annealing, have received a lot of attention ^[1,2]. One of the recently developed population based optimization methods is the immune algorithm (IA), which is based on the biological immune system theory^[3,4,5]. IA simulates the mechanisms of the biological immune response when a biological immune system is exposed to an antigen. Clonal proliferation ^[6], negative selection ^[7], immune network ^[8], danger theory ^[9], and dendritic cell model ^[10] are the most commonly used mechanisms of the biological immune response. The proposed technique illustrated in this paper is mainly based on the clonal selection and proliferation principle ^[11,12,13]. IA has been successfully applied to many aspects such as function optimization ^[14], control ^[15], learning ^[16], pattern recognition ^[17], engineering optimization ^[18,19], protein structure prediction ^[20,21], and multiple sequence alignment ^[22].

In general, IA is used to solve combinatorial optimization problems due to its learning capability and self-organizing. Certain important issues should be considered to design an excellent IA. One such issue is maintaining the population diversity as the search progress. IA begins with various candidate solutions and manipulates them traversing the search space to achieve better solutions by some operations, i.e., mutation, selection and clonal proliferation. If the diversity of population is too limited to facilitate the algorithm visiting more unexplored search regions, the premature convergence of algorithm occurs, which is likely to result in local optimal trapping problems. To alleviate this problem, Gao et al.^[23] has proposed an IA with feedback mechanisms to effectively handle the population size and thus maintain the diversity well. Other attempts on resolving this problems may refer to ^[24,25].

Another key issue while designing IA is to balance exploitation and exploration of the search in decision space. Exploration is the process of visiting the entirely new regions of a search space, whereas exploitation is the process of visiting those regions of a search space through the neighborhood of previously visited points ^[26]. In order to design a successful algorithm, it is necessary to establish a good ratio between exploration and exploitation. After realizing this, Khilwani et al. ^[27] proposed a fast immune algorithm (FCA) by combining cauchy mutation (CM) with gaussian mutation (GM). CM enables the search to carry out long distance jumps, thus specializing in exploring unvisited regions ^[28]. Besides, GM is more capable of exploitation in small regions of a smooth decision space. Nevertheless, due to the randomness of the search dynamics both in CM and GM, the performance of the algorithm may not achieve better convergence speed. One possible method to improve the convergence speed of the algorithm is the usage of more enriched search dynamics. In ^[29], Lee et al. has proposed a new mutation operator, which is based on Lévy distribution to generate an offspring that is far away from its parent, aimed to search the decision space much faster. Moreover, information exchanges (i.e., crossover) between two different solutions may lead a fast assimilation for the evolution. Most of the solutions undergone information exchanges will possess some common information (or sub-solution), therefore accelerating the convergence speed. In ^[30], Gong et al. has proposed a Baldwinian learning operator to enable the current solution learn from the co-adapted alleles in certain environments. Motivated by the mechanism of idiotypic network ^[31], a lateral interaction mutation (LM) has been put forward to effectively realize the knowledge exchange during the period of getting solutions ^[32].

Essentially, the evolution of solutions in IA is driven by the selection and mutation operators. Selection promotes the searching towards the regions of the best individuals. Heavy selection pressure urges the search towards more exploitation, and little selection pressure pushes the search towards more exploration ^[33,34]. Mutation enables the original solutions to move toward different regions. On the one hand, mutation randomly modifies individuals with a given probability, and thus increases the structural diversity of a population. From this point of view, a mutation operator is more of an exploration operator. Such an operator (CM, GM, etc.) facilitates the recovery of genetic diversity which would be lost during the selection phase, and explores new solutions. On the other hand, mutation can also refine an individual by conditionally incorporating the information derived from other individuals into it, aiming to generate the possibility of better offspring and faster convergence. From this perspective, such a mutation operator (LM, etc.) is more of an exploitation operator. As stated above, it is crucial to balance exploration and exploitation of the search. As a result, apart from the selection operator, two key questions should be considered to design a positive IA: 1) which kind of mutation operators should be used in the algorithm, and 2) what ratio of each mutation should be implemented by properly managing the control-parameter settings.

To address the above problems, a new parallel mutation (PM) framework is proposed in this paper. In PM, three distinct mutation operators (GM, CM and LM) are strictly selected to carry out the individual evolution. Due to its fine-grained search characteristic, GM is responsible for the exploitation within a small search neighborhood. Owing to its coarse-grained search property, CM is capable to perform hill climbing mutations, and thus it will be beneficial to make the algorithm jump out of the local optima. Assisted by the surrounding individuals, the individual undergone LM is manipulated to learn knowledge from the neighborhood, therefore a fast convergence speed would be realized. During the implementation of the algorithm, each single mutation operator in PM is assigned an operation probability and a randomly generated probability is used to control which operator will be carried out to perform the search. By making use of the introduced PM framework, an effective immune algorithm (IIA) is proposed. It is expected that IIA can well balance exploration and exploitation of the search, by means of not only the reduction of the number of useless or redundant iteration number, but also the guidance of the search to promising regions which are close to the global optima with a high probability.

Lennard-Jones potential problem (LJPP) is a potential energy minimization problem, where neutral atoms are subject to two distinct forces under the limit of large distance and short distance: a dispersion force at long ranges, and a repulsion force at short distance, resulting from the overlapping electron orbits. In LJPP, the objective function is non-convexity, and there are huge number of local minima which is growing exponentially with the number of the atoms. Due to the excessive number of local minima for LJPP, it make extremely hard to find a global minimum of the energy. In fact, LJPP belongs to the classical NP-hard problems ^[35,36] for which no algorithm is guaranteed to find the global optimal. Hence, the traditional local optimization methods may not be efficient for such problems. At the earliest, Northby firstly introduced a effective and successful methods to LJPP ^[37] and further refined in ^[38]. This method is beginning local optimization with initial configurations which are produced by randomly placing atoms in predefined points in space. Then, a large number of efficient approaches are used to solve LJJP, such as basin hopping method ^[39], the big-bang method ^[40], Leary's descent method ^[41], genetic algorithms ^[42,43,44], differential evolution ^[45] and particle swarm optimization ^[35]. A geometric crossover technique has been introduced in ^[46,47].

To verify the effectiveness and applicability of the algorithm, IIA is testing on a two-dimension function, some benchmark functions and a real-world LJPP. Benchmark functions are often used in optimization scenarios, while LJPP is more practical and is hard to be solved. Thus, IIA is tested from two different aspects. Via the function optimization, especially a multimodal function, the capability of jumping out of local minimum of IIA is well illustrated and demonstrated. While via resolving the real-world optimization problem, IIA has a certain potential to deal with the larger scale problems. Experiments are conducted based on the above two problems, and the illustrative and comparative results could clearly show the effects of PM on IIA.

The rest of the paper is organized as follows. Section 2 describes the the model of LJPP. Section 3 introduces the PM framework and presents the details of IIA. The experimental results and discussions of IIA whenever solving the benchmark functions and LJPP are shown in Section 4. Finally, Section 5 presents some general remarks to conclude the paper.

The contributions of this study are as follows: 1) We propose a novel algorithm for improving the performance of search algorithms. 2) We demonstrate the effectiveness of the proposed algorithm through extensive experiments on benchmark problems. 3) We provide a comprehensive analysis of the proposed algorithm, including its advantages and limitations. 4) We compare the proposed algorithm with existing state-of-the-art algorithms in the literature and show that it outperforms them in terms of search performance. 5) We discuss the potential applications of the proposed algorithm in various domains.

2.   Lennard-Jones potential problem

Lennard-Jones potential problem (LJPP) is a potential energy minimization problem and it involves the minimization of molecular potential energy associated with pure Lennard-Jones cluster ^[45,48]. LJPP is a multi-modal optimization problem comprised of an exponential number of local minima ^[48]. The lattice structure of Lennard-Jones cluster has an icosahedral core and a combination of surface lattice points. Most of the global minima which have structures based upon the Mackay icosahedrons can be found in the Cambridge Cluster Database. An algorithm can be tested over this function because of its capability to conform molecular structure, where the atoms are organized in such a way that the molecule has minimum energy. Lennard-Jones pair potential for $D$ atoms which is given by the Cartesian coordinates

is represented as follows:

where

with gradient

Lennard-Jones potential has the minimum value at a particular distance between the two points ^[49]. The variation of Lennard-Jones pair potential $V(r) = r^{-12}-2.r^{-6}$ with pair distance $r$ have been plotted in Figure 1. It can be clearly found that it has a unique minimum value of -1 for a particular value of $r$ = 1. Its value increases very fast when $r$ decreases from the optimum value and tends to infinity near $r = 0$. This tendency of pair potential curve makes it turn to be a very hard optimization problem.

3.   Immune algorithm with parallel mutation

In the literature, several IA models inspired by the mechanism or function of immune cells in the human immune system have been developed ^[12,50,51]. Among them, most IA techniques emulate the clonal selection principle which depicts the response procedure when human being is exposed to the extrinsic substances such as the bacteria and viruses. In clonal selection principle, two key operations take place, i.e., the differentiation and proliferation of the immune cells. From the algorithmic perspective, differentiation can be viewed as mutation operator of individuals, while proliferation be regarded as reproduction of individuals. The generic computational framework of IA can be depicted on the left and the proposed PM mechanism on the right in the Figure 2. Intuitively, the quality of mutated populations is determined by the capacity of mutation operators ^[52,53], which directly influences the performance of the algorithm. In other words, mutation drives the population to evolve more promising individuals. In order to make the mutation operator possess more search dynamics, a parallel mutation framework is introduced in this paper. Three distinct single mutation operators are strictly chosen to be implemented in the PM.

3.1.   Mutation operators

To make the paper self-explanatory, three mutation operators including GM, CM and LM are introduced briefly and respectively. To make the description more concise, We combine the representation of the existing population of people to be $\{X_1, ..., X_i, ..., X_N\}$, where the $i$-th individual $X_i = (x_{i1}, ..., x_{ij}, ..., x_{iD})$, $N$ stands for the population's size, whereas $D$ stands for the current optimization problem's dimension. As to LJPP, $D$ is the number of atoms. Being manipulated by mutation operators, the mutated population is represented by $\{X_1', ..., X_i', ..., X_N'\} = P\{X_1, ..., X_i, ..., X_N\}$, where $X_i' = (x_{i1}', ..., x_{ij}', ..., x_{iD}')$, $P$ represents the mutation used in IIA, and $P\in$ {CM, GM, LM}.

Gaussian mutation (GM): the mean value $\mu$ and the standard deviation $\sigma$ are two GM characteristics that are utilized to calculate the mutation's step size. The one-dimension gaussian density function with the mean value $\mu$ and the standard deviation $\sigma$ is defined by:

Like many other researches ^[27,54,55], the origin-centered gaussian density function is the foundation of the simplified GM mutation strategy when $\mu = 0$, $\sigma = 1$. In a single iteration of the population, GM is implemented based on the following equation:

where $N_j(0, 1)$ is a gaussian random number that is normally distributed and has a mean of 0 and a standard deviation of 1. The population's serial number is represented by the symbol $i$, and the dimension is represented by the symbol $j$. $s_g$ is the step size of mutation, and it is generated by ^[27,55]:

where $w_g$ is a random number that is generated uniformaly across the range of $[0, f_{gau}(0)]$.

Cauchy mutation (CM): the cauchy distribution, as depicted in Eq (3.4), is the primary foundation of CM.

The absence of an expectation is the most significant feature of cauchy distribution ^[27,55]. The shape of $f_{cauchy}(x)$ is identical to that of the gaussian density function, but it moves very slowly toward the axis. The CM is implemented as follows ^[27]:

where $\delta_j$ is a cauchy random variable. $s_c$ is the step size of cauchy mutation which is generated by:

where $s_c$ is a uniformaly generated random number in the range $[0, f_{cauchy}(0)]$.

Lateral mutation (LM): the lateral mutation according to the idiotypic network theory ^[8,32], lateral interaction typically occurs in immune systems. In other words, not only a foreign antigen but also the external idiotopes can identify each paratope on a person. Motivated by this mechanism, LM is implemented as:

where $k\in \{1, 2, ..., N\}$ and $k\neq i$. The learning rate $\beta \in (0, 1)$ is a real number generated randomly.

The reason why we choose the forementioned three mutation operators to be involved in PM is their distinct search dynamics. In Figure 3, we summarize the characteristics of each mutation operator respectively. The optimization problem's solution space is depicted by the solid rectangle $S$. The dashed circles denote contour lines of affinity, and the inner circles indicate that they represent higher affinities than the outer ones. In this figure we can notice that individuals mutated by GM and CM only utilize random perturbation on the individual itself, while those mutated by LM can make use of the information from the environment. As remarked by Yao et al. ^[55], Due to its long, flat density function tails, CM is more likely than GM to produce offspring far from its parent. It is more likely to carry out large variations in individuals because it is possible to execute longer jumps with a great probability. As a result, it is frequently used as a learning strategy to address the local optima trapping issue. However, as stated in ^[30,56], learning from one's surroundings offers a positive alternative and, most likely, a simpler approach to improving one's search performance. In particular, the mechanism in LM directs the current search by utilizing the data for a randomly selected member of the population. Successful guidance necessitates a significant amount of redundant searching if the chosen guide is far from the global optimal solution, as the quality of the chosen guide is crucial. It is worth emphasizing that only one or two single mutation operators can not fully realize the exploration and exploitation of the search. It is suggested that a well-designed mutation operator should contain at least three properties: 1) a fine-grained search ability derived from GM, 2) a coarse-grained search aptitude rooted in CM, and 3) a fast convergence manipulated by LM.

3.2.   Framework of the PM mechanism

In order to incorporate the three above mentioned complementary search dynamics in the algorithm simultaneously, a parallel mutation (PM) framework is proposed in this paper. The pseudocode for implementing the mutation mechanism is shown in Table 1. Three parameters $Prob_{C\!M}$, $Prob_{G\!M}$ and $Prob_{L\!M}$ satisfying Eq (3.8) are used to control the probability of the implementation of each mutation operator.

In the process of solution evolution, the current individual $X_i$ may mutate to $X_i'$ by means of GM, in the next generation, the individual $X_i'$ may mutate to $X_i''$ by means of CM, then the individual $X_i''$ may mutate to $X_i'''$ by means of LM. Only one mutation operator is used in each iteration of the individual maturation, and the mutation operator that is used is chosen at random. However, PM will have all of the search dynamics of the embedded single mutation operators because, after a large number of iterations, each single mutation operator will be implemented multiple times.

In the PM framework, each embedded single mutation operator will be carried out if its corresponding cumulative probability is larger than a randomly generated number, namely $q$, thus it is deduced that this operator's search dynamic contributes to the search of the entire algorithm. On the other hand, rather than being implemented deterministically, the single operator is implemented at random due to the stochastic nature of $q$. In addition, the fact that the implementation probabilities for each mutation operator are adaptable suggests that resetting these probabilities can effectively coordinate the search dynamics of the resulting PM. Several remarks can be concluded regarding PM:

1) PM inherits GM's ability to exploit individuals and can conduct fine-grained searches within a relatively small neighborhood.

2) Due to the occasional large mutations in CM, PM is able to conduct coarse-grained searches in a wider area around the individual. In some cases, these larger mutations are advantageous in that they cause the algorithm to jump out of the local optima.

3) The environment of LM can provide PM with information, such as from a guiding individual. Such a mutation would be very effective in accelerating convergence if the guiding individual is already close to the global optima, and it would also be very possible to improve the average quality of solutions.

3.3.   The suggested probability of the mutation operators

It is important to set the values of the operation probabilities for each single mutation parameter. Such combination of setting values puts significant influence on the performance of the algorithm. Usually, the trial-and-error approach, which is a time-consuming and tedious method, can be performed in an ad-hoc manner to set these values. In this paper, two kinds of parameter setting strategies are given. One is a simpler approach, the suggested probability of each mutation operator is equal. The other is affinity-based rule to determine the values. These two approaches will be introduced and performed on benchmark functions and LJPP to compare the performance, respectively. After setting the initial probability of each mutation operator, the adaptive mechanism has been proposed in our algorithm. The learning probability values of each mutation operator are gradually changing according to the iterations. The general form is given as following:

in which, $Prob_{C\!M}, Prob_{C\!M}, Prob_{L\!M}$ represents the participation probability of each mutation operator respectively, $T_{max}$ is the number of maximum iterations, $t$ is the current iteration.

1) Simpler approach with a gradually decreasing learning values (PMGD).

This is a very simple method which is not related with the search ability and search dynamic mechanism of each mutation operator. In this article, we design the strategy only related to the number of iterations and the current iteration. At the beginning of the experiment, the probability of each mutation operator is set to equal ($1/3$).

2) Approach depending on the fitness values of each mutation operator (PMDF).

After preliminary experiments, the optimal performance of each single mutation operator is verified. The principle of the deterministic setting rule is that the more powerful the search capacity of the mutation is, the larger operation probability it possesses, and vice versa. It is expected that the exploration and exploitation of the search can be well balanced by setting operation probabilities in such a manner. Further experimental evidences are presented in Section 4. The process of probability decision is as follows:

In the preliminary experiment, the search capacities of each single mutation operator are observed based on 30 independent runs. The mean objective function values obtained by CM, GM, LM separately are represented as $f_{C\!M}$, $f_{G\!M}$, $f_{L\!M}$ respectively. The maximum value among them are computed to normalize the operation probability.

where $\left\lceil {f_{max}} \right \rceil$ rounds up to the closest integer that approaches the value of $f_{max}$. As for minimization problems, the adjusted function values $f'_{C\!M}$, $f'_{G\!M}$, $f'_{L\!M}$ are obtained by reducing $\lceil f_{max} \rceil$, where $abs\{.\}$ means the absolute value of its variable. The deterministic setting rule of these operation probabilities are based on the following equations.

It is clear that the above three operation probabilities belong to the interval of $[0, 1]$ while satisfying Eq (3.8). The more fit individuals a mutation operator generates, the larger operation probability it possesses. By doing so, the three distinct mutation operators can be well merged together in PM to achieve a better search performance. It should be noted that, although it is possible to make these operation probabilities self-adapted or evolved by algorithmic parameter's coordination ^[57], such advanced setting methods are too sophisticated to design, which are even more difficult than the design of the optimization problem itself. As a result, general guidelines on how to set these operation probabilities are provided based on Eqs (3.12)–(3.18) after the search features of each mutation operator are identified.

3.4.   Procedures of IIA

The general procedure of IIA consists of five components as shown in Figure 2. The procedure involving initialization, affinity evaluation, proliferation, mutation, and selection are regarded as one generation of the evolution of population. The complete algorithm with PM incorporated is stated here.

1) Initialization: on the part of the function optimization problems, randomly initializes $N$ individuals in the $D$-dimensional decision space $P_i$, $i = 1, 2, ..., N$, where $P_i = (p_{i1}, ..., p_{ij}, ..., p_{iD})$ is the $i$th individual in the $D$-dimensional decision space. It should be noted that the representation of solutions for LJPP is slightly different from that of the function optimization problems. For LJPP, the search dimension size is $3\times D$, where $D$ denotes the number of the atoms, and $3$ means the three Cartesian dimension that atoms locate. The upper and lower bounds of the $i$th dimension in the decision space is $[low_j, up_j]$. To make the randomly generated individual be a feasible solution for the problem, the following formula is used.

where $\mu \in [0, 1]$ is a random number with the uniform distribution.

2) Affinity Evaluation: compute the affinity values $f(P_i)$, $i = 1, ..., N$, for each individual in the population. In general, an individual with larger affinity value indicates a better solution for the problem. For minimization problems, the affinity values is reversely proportional to the function values (or energy values in Eq (2.2)):

3) Clonal Proliferation: for each individual $P_i$, asexually produces $m_i$ clones where

where $M$ is the clone size, $\lfloor.\rfloor$ is the operator that rounds its arguments towards the closest integer. $m_i$ denotes the clone number of the $i$th individual after all the individuals in an descending order of the affinity are sorted. It is obvious that the amount of generated clones of an individual is inversely proportional to its affinity. In other words, the better the individual is, the more clones it produces. Afterwards, the population of individuals is represented by $(P_{1, 1}, P_{1, 2}, ..., P_{1, m_1};...; P_{N, 1}, P_{N, 2}, ..., P_{N, m_N})$;

4) Mutation: for each clonal individual $P_{ij}$, $i = 1, ..., N$, $j = 1, ..., m_i$, mutates it according to Algorithm 1. After mutation, the population becomes $(P'_{1, 1}, P'_{1, 2}, ..., P'_{1, m_1};...; P'_{N, 1}, P'_{N, 2}, ..., P'_{N, m_N})$. At this point, it is worth indicating that the mutated individual goes out of the search bounds is re-initialized by using Eq (3.19) to make it to be a feasible solution.

5) Selection: evaluates the affinity of each individual in the mutated population to obtain the updated population $P_i$, $i = 1, 2, ..., N$, where

6) Termination condition: The above steps 2)–5) would be repeated until the termination condition (i.e., the iteration number reaches $T_{max}$) is satisfied. Finally the optimal individual is outputted.

4.   Experimental studies

4.1.   Parameter setting

The algorithm's performance is strongly influenced by parameter values. Before it can be used to optimize problems, this paper's user-defined parameters need to be well coordinated. Table 2 summarizes the parameters used in IIA. $N$ denotes the population size of the algorithm. Large value of $N$ results in more individuals involving in performing search in the decision space, thus might achieve better search performance. However, large value of $N$ also implies more implementation times of mutation and selection of the operators, thereafter it requires more computational time. As suggested in ^[11,58], a trade-off between the quality of solutions and computational time is to set $N$ as 50. The clonal size can be controlled by parameter $M$. Large value of $M$ indicates more clones produced in the algorithm, and the more fit the individual is, the more plentiful clones it generates. As all $m_i$ clones of each individual $P_i$ will undergo the progress of the parallel mutation, $P_i$ will have $m_i$ times of opportunities to mature (i.e., become a more fit individual for the problem). Nevertheless, too many times of mutation for an individual might guide the search implemented in the same region, thus would easily cause diversity loss and pre-maturation of the population. Preliminary simulation results suggest that a compromising value for $M$ is 10. Moreover, the parameter $D$ means the search dimension which is problem-oriented.

The three operation probabilities in PM ($Prob_{C\!M}$, $Prob_{G\!M}$ and $Prob_{L\!M}$) are difficult to determine. In the previous section, we have proposed two types of determined mechanism for operation probabilities. The first category is easy to perform and only depends on the current iteration, however, the second type is different from the first one and based on the search abilities of every operator. We will illustrate these two mechanisms in the following.

1) The initial participation probability of CM, GM, LM with PMGD mechanism

The initially participation probabilities of the three mutation operators in the PMGD mechanism are equal. Therefore, the probability is set 1/3.

2) The initial participation probability of CM, GM, LM with PMDF mechanism

The participation probabilities of the three mutation operators depend on the fitness values of search capability in the PMDF mechanism. Therefore, the fitness value of a single mutation operator must be obtained and the participation probabilities are calculated according to the fitness values proportion.

The deterministic values are given below. Because the search performance of the algorithm strongly depends on their values. We present some guidelines on how to set them after the search capacity of each single mutation operator is observed. The preliminary experiments are conducted based on the LJPP with two different groups of atoms. Table 3 records the final affinity values $f$ (Eq (3.20)) obtained by every single mutation operator, the resultant amendatory affinity values $f'$ (Eqs (3.13)–(3.15)), and the corresponding operation probabilities $Prob$ (Eqs (3.16)–(3.18)). From this table, it is obvious that the operator $LM$ possesses the best search performance, $GM$ is the second, and $CM$ is the worst. The results are consistent with the search properties of each single mutation operator. $LM$ enables individuals to exchange information with each other, thus accelerate the convergence speed to acquire better solutions. $GM$ modifies the individual to carry out a local exploitation within the neighborhood slightly, and improves the solution quality gradually. $CM$ is able to visit the search space by jumping to a faraway location suddenly, but it can not focus on exploiting a small region, therefore its search performance is the worst. Actually, the search dynamics of the three mutation operators in PM are complementary. After observing such a phenomenon, the general guidelines of setting the operation probabilities can be summarized as:

1) None of the single mutation operator can be omitted in PM which means the probability can not be set to zero.

2) According to the search capacity, the ranked order during operation probabilities is $Prob_{C\!M} < Prob_{G\!M} < Prob_{L\!M}$.

Any combination obeying the two guidelines can be regarded as a suitable setting, although the actual setting is problem-oriented and it might follow Eqs (3.16)–(3.18) in every iteration. In the experiment, these values are settled as $Prob_{C\!M}$ = 0.1, $Prob_{G\!M}$ = 0.3 and $Prob_{L\!M}$ = 0.6 in PMDF mechanism.

3) Other possible participation probability among three mutation operators

Tables 4 and 5 record the statistical values of solutions obtained by IIA with different learning probabilities in PM for the unimodal function $f_1$ and the multimodal function $f_3$, respectively. These learning probabilities are taken from eight representative design points, where a, b and c in "a/b/c" denote the assigned implementation probabilities of CM, GM and LM respectively. Tables show that the setting of the implementation probabilities in PM using "0.1/0.3/0.6" exhibits the best learning performance in terms of solution qualities. It is worth emphasizing that, although "0.1/0.3/0.6" might not be the optimal setting of the probabilities, it outperforms most of the other settings due to the usage of experimental design with mixtures. However, the initial participation probability is 1/3 for each mutation operator which is significantly better than single mutation operator or two mutation operators. Therefore, we can conclude that the IIA algorithm is the best when three mutation operators participate simultaneously, and outperforms single or two mutation operators.

4.2.   Identification of search characteristics on function optimization

4.2.1.   Benchmark functions

In this paper, some benchmark functions with different categories are selected to evaluate the performance of the proposed algorithm. $f_1$ and $f_2$ are unimodal functions, $f_3$ and $f_4$ are multimodal functions with many local minima, $f_5$ and $f_6$ are multimodal functions with a few local minima. The function formula of them is given in Table 6 and the more detailed explanation of each function is listed in Appendix in the reference [55]. The symbol "Dim.'' denotes the dimension of the function. The domain $S = [low, up]^D$ represents the search space of feasible solutions for the $D$ dimension problem. $f_{min}$ is the global optimal solution. It should be noted that the results obtained are under the settings for parameters by, $N = 30$, $M = 5$, $T_{max}$ = 2000 for high-dimensional functions and $T_{max} = 100$ for low-dimensional functions in all compared algorithms.

The first set of experiments aims to compare the search dynamics of the parallel mutation mechanism with the single mutation operators and the performance between the PMDF approach and PMGD strategy. To realize this, we constructed three variants of IIA using only one single mutation operator by turns. IIA-1, IIA-2 and IIA-3 represent the utilization of $C\!M$, $G\!M$ and $L\!M$ in the algorithm respectively. Due to the different deterministic mechanism of parameter setting, IIA and IIA-4 record the PMDF approach and PMGD strategy to determine the operation probability, respectively. Table 7 summarizes the average results of 30 independent runs obtained by IIA, IIA-1, IIA-2, IIA-3 and IIA-4. It is obvious that IIA with parallel mutation mechanism and the PMDF strategy has the best mean value and standard deviation. Figures 4–6 show the comparative results by means of convergence graphs and box-and-whisker diagrams of solutions over 30 runs for tested functions. From Figures 4–6, it is apparent that IIA performs better than the other three mutation operators in terms of convergence rate and solution quality, suggesting that PM has better dynamic searching characters than any single mutation operator. Hence we can conclude that the hybridization of different single algorithm is a feasible and promising method to enrich and promote the searching performance.

As for unimodal functions, the convergence speed is more interesting than the final optimization result. In Figure 4, IIA-3 has the fastest convergence speed at the beginning among all algorithms, indicating that $LM$ is able to use information from other antibody to speed up the convergence, but frequent information exchange causes the algorithm fall into the local optimum easily. IIA-1 and IIA-2 have slow convergence rates because $CM$ and $GM$ only use a random perturbation in the process of solution maturing. However, it has the ability to jump out of local optimum. IIA and IIA-4 contains the advantages of the three mutation factors, therefore it has the second and the third convergence speed at the beginning. At the same time, IIA has the fastest convergence rate and obtains the best optimization value at the end of the algorithm, in other words, IIA with PMDF approach is significantly better than IIA-4 with PMGD strategy for unimodal functions. For multimodal functions, the finial results are much more important since they reflect an algorithm's ability to escape from local optima. From Figures 5 and 6, it is evident that IIA has the best optimization results than the other three single mutation operators for all the functions. To sum up, parallel mutation operator is more effective than single mutation operator in searching better solutions and it has complementary search effects on one another no matter how inferior their own search capacities are. It is also obvious that the experiment result obtained by IIA with PMDF approach is more excellent than IIA-4 with PMGD strategy for all the multimodal functions. That is to say, the deterministic probability mechanism of PMDF is suitable in this algorithm.

The results compared with the existing IAs, i.e., SIA, FCA ^[27], PSO and CLONALG ^[56] are summarized in Tables 8 and 9. Each result of IIA is obtained from 30 independent runs. It is obvious that IIA outperforms FCA, SIA, PSO and CLONALG on multimodal functions $f_3$–$f_6$. For unimodal functions $f_1$ and $f_2$, IIA is superior to SIA, but it is inferiorer than FCA. So we can conclude that IIA has advantage in solving multimodal problems. In the following sub-section, IIA is applied on LJPP.

4.2.2.   Two-dimension function

Problem $G_1$ ^[59]:

In this section, the search characteristics of the proposed IIA are verified on a two-dimension function shown in Eq (4.1). This problem is a minimization problem. Figure 7 depicts the 3D plot of the function and initial solutions for problem $G_1$ generated by IIA. The global minimum of this function in the interval $0 < x, y < 10$ is located at (0.9039, 0.8668) and the minimum value is -18.5547. It is visible that $G1$ is a multimodal function which has many local optima.

In order to identify the search characteristics of every single mutation operator and the proposed PM, a comparative experiment is conducted. The above four mutation operators (i.e., PM, CM, GM, and LM) are adopted one by one within the same algorithmic framework of IA. The same initial fifty solutions (as shown in Eq (4.1), where the star marks represent the solutions) which are randomly generated in the search space are utilized as the initial population of individuals. Figures 8–11 illustrate the distribution of solutions obtained at the iteration of the 5th, 10th, 30th and 50th by using the four mutation operators respectively. At first glance, the solutions in all figures concentrate on the nearby locations in the search space along with the iteration number's increment, indicating that all the constructed algorithms converge to the promising search areas gradually. It is apparent that PM enables the algorithm to possess the best search performance. Observing from the contour graph in Figure 8, it is clear that all solutions move to the local minima fast even at the first 5th iteration. At the 10th and 30th iteration, the convergence has been maintaining, and finally at the 50th iteration all solutions achieve the global minimum point and a local minimum point. Contrarily, the final solutions obtained by CM, GM and LM are trapped by more than one local minima. Thus, it can be concluded that PM is capable to manipulate individuals so as to find the global optima with a very greater probability and to be trapped by local minima with a smaller probability than any other mutations which used single mutation operators.

Furthermore, the empirical results shown in Figures 9–11 support the above intuition in Section 3.1. Compared with CM and GM, LM exhibits better search performance in terms of solution quality and convergence speed, which sustains that information exchange among individuals can lead to a fast convergence. On the other hand, too much information exchange might also make the search premature, which means that LM is lack of exploration capacity to jump out of the current search region. Besides, CM's and GM's behaviors on $G1$ at the 10th iteration are quite illuminating. At the 10th iteration, both of the solutions obtained by GM and CM are widely distributed in the whole search space, which indicates their slow convergent ability. However, compared with their corresponding solutions at the 5th iteration, most individuals mutated by GM move slightly towards the local minima (about 17 clusters in the search space), while those mutated by CM make long distance jump and thus own better local minima climbing capacity (about 12 clusters).

Overall, the illustrative results of the search dynamics verify that PM which combines the three distinct search characteristics of GM, CM and LM is such a good methodology to construct the more powerful mutation operators.

4.3.   Performance verification on LJPP

To further demonstrate the effectiveness and efficiency of the proposed IIA, LJPP obtained from the Cambridge Cluster Database has been used as the test suit. Firstly, the proposed IIA algorithm has been compared with the IIA-1, IIA-2, IIA-3, IIA-4 on 2, 3, 10, 38 atoms and the optimal solutions of them are already known. The optimal energy (i.e., fitness) and the maximum iteration number of the search are also set in Table 10. With the number of atoms increasing, the superiority of IIA becomes more significant which suggests IIA is more suitable for solving large scale problems. Besides, it is worth reiterating that the experimental results of IIA-1, IIA-2, IIA-3 and IIA-4 on LJPP still support our previous statements. In order to evaluate IIA more overally, the convergence graph and the box-and-whisker diagrams of solutions are also depicted for 2, 10, 38 atoms in Figure 12. It is apparent that PM outperforms the other mutation operators in terms of convergence speed and solution qualities. Due to the small data set, we have expanded the dataset and the optimal solutions of some expanded atoms are not known. However, it does not affect the performance of comparing different algorithms. We compared IIA (PMDF) with IIA (PMGD) algorithms with atomic numbers of 2, 3, 10, 15, 20, 25, 30, 38 and 50, and the experimental results are shown in the table below. The statistical results obtained on the basis of 30 independent runs are summarized in Table 11. The number of atoms is up to 50 which means the dimension of the search space is 150. The experimental results recorded are the average fitness (Mean) and the standard deviation (Std) over 30 runs. It is obvious that IIA performs significantly better than its competitors for all the tested instances consistently. IIA can find the optimal solutions for the LJPP with 2 or 3 atoms, and the Std is the smallest which indicates the robustness of the search performance.

Table 12 gives the t-test values between IIA (PMDF) and the competitors including IIA-1, IIA-2, IIA-3 and IIA-4 (PMGD). All the values are less than 0.05 which means that the average values are equal for small probability. We believe that the average values of IIA (PMDF) are significantly different from IIA-1, IIA-2, IIA-3 and IIA-4. Thus, the proposed PM has been demonstrated to be an effective methodology which combines different single mutation operators together to be a more powerful one.

Finally, Table 13 gives the computation times of the algorithm were calculated and compared for the algorithm at different numbers of atoms. Although the algorithm's computation time is slightly longer than the GA and PSO algorithms, the longitudinal comparison shows that the algorithm does not increase the time complexity of the algorithm too much as the number of atoms increases.

5.   Conclusions

This paper proposes an effective immune algorithm (IIA) with parallel mutation (PM) and evaluates its performance on a two-dimension function, some benchmark functions and a number of Lennard-Jones potential instances. The experimental results show that IIA performs much better than its variant algorithms which use single mutation operators (i.e., CM, GM, LM) or a combination of CM and GM. The paper analyzes IIA in depth in terms of the user-defined parameters and the search dynamics of each incorporated mutation operator. The empirical evidences suggest that PM is an effective methodology which combines the three distinct but necessary search dynamics: 1) a fine-grained search ability derived from GM, 2) a coarse-grained search aptitude rooted in CM, and 3) a fast convergence manipulated by LM. The paper also investigates the setting for the ratios of operation probabilities. Some insights and guidelines in general have been concluded from the experimental results. The idea of PM can also be applied to the other population-based algorithms, such as evolutionary algorithms, swarm intelligence, etc. Future program on IIA will include the incorporation of more different mutation operators into PM and the self-adaption of the operation probabilities.

Conflict of interest

The authors declare there is no conflict of interest.