Global dynamics of alcoholism epidemic model with distributed delays

1.   Introduction

Global optimization problems widely exist in engineering application and scientific research. For linear and convex function optimization problems, the theories and corresponding optimization algorithms are well-developed. However, these algorithms often fail to solve non-convex and non-differentiable problem with many optima. Since 1980s, intelligent optimization algorithms or so-called mate-heuristic algorithms have been developed. They can obtain near-optimal solutions even optima in reasonable time. Classic meta-heuristic algorithms include Simulated Annealing (SA) ^[1], Genetic Algorithm (GA) ^[2], Particle swarm optimizer (PSO) ^[3], Differential Evolution (DE) ^[4,5], Ant Colony Optimization (ACO) ^[6], Electromagnetism-like Mechanism Algorithm 7 etc. They have been used in extensive applications in areas such as engineering optimization ^[8,9,10] and scheduling optimization ^[11,12,13] etc.

Particle swarm optimizer (PSO) was first proposed by Eberhart and Kennedy in 1995 3, and since then, it has become a very popular algorithm since PSO shows many advantages such as simple implementation, good optimization capability and minimum mathematical processing. In the review paper ^[14] about nature-inspired intelligence, PSO was considered one of the most popular nature-inspired algorithms according to related publications. However, the original PSO algorithm may get trapped easily in local optima for multimodal optimization problems. In order to improve the search ability of PSO, an increasing number of modified PSO versions have been proposed, including PSO with inertia weight (PSO-w) ^[15], PSO with constriction factor (PSO-cf) ^[16], fully informed particle swarm (FIPS) ^[17], unified particle swarm optimization (UPSO) ^[18], comprehensive learning particle swarm optimizer (CLPSO) ^[19], cellular particle swarm optimization (CPSO) ^[20], and adaptive division of labor PSO (ADOLPSO) ^[21] etc.

Among these modified versions of PSO, the comprehensive learning particle swarm optimizer (CLPSO) proposed by Liang et al. 19 is a quite famous one which shows a superior ability for solving multimodal problems. Since CLPSO was published in 2006, it has been cited for 2500+ times in google scholar until now. Some researchers have developed new algorithms based on CLPSO. Peng and Lu ^[22] proposed a hierarchical PSO based on CLSPO. Gülcü and Kodaz ^[23] proposed a better variation of CLPSO which is named parallel comprehensive learning particle swarm optimizer (PCLPSO). These studies show that CLPSO has a strong vitality in the area of metaheuristic algorithms but there is still room for improvement.

Adding other local search methods to PSO is a very popular way to improve its performance. Bao, Hu and Xiong ^[24] combined PSO with pattern search for SVMs parameters optimization. Petalas, Parsopoulos and Vrahatis ^[25] used the random walk with directional exploitation local search in PSO. Jia et al. ^[26] combined the canonical PSO with a chaotic and Gaussian local search procedure. The cellular particle swarm optimization (CPSO) ^[20] essentially used a local search based on cellular automata. Wu et al. ^[27] proposed a superior solution guided particle swarm optimization (SSG-PSO) combined with four gradient-based or derivative-free local search methods. In these algorithms, the local search usually focuses on enhancing the ability of intensification. Is there a local search method can also improve the capability of global search? To answer this question, in this work, Latin hypercube sampling with a control method is combined as the local search method.

Latin hypercube sampling (LHS) is a stratified sampling method which was proposed by Mckay in 1979 ^[28]. Compared with the Monte Carlo sampling, Latin hypercube sampling is more efficient and able to reflect the population more accurately from fewer sample points. Since the LHS has strong global search ability, symmetric Latin hypercube design, a variant of LHS, was used for population initialization in differential evolution ^[29]. Recently, LHS was successfully applied in model updating of a historic concrete bridge at two different levels ^[30]. On one hand, LHS was implemented with a sensitivity analysis phase to define parameters which mostly influence the FEM modal response. It was also used to define initiation points for another optimization algorithm called Trust-Region. The application also takes advantage of its global search ability. LHS was also used as local search method for a hierarchical PSO ^[22]. In their work, the length of each dimension of the hypercube is two times the length of the corresponding dimension of the selected particle ^[22]. The setting method will make the size of the hypercube be determined by the location of the selected particle.

Actually, an optimization algorithm can be treated as a sampling method from the view of sampling. It uses some rules to a sample from the whole space of independent variables and tries to find the global optimal point. Therefore, Latin hypercube sampling can be available in CLPSO. In this work, it is used as a local search method for the global best particle in CLPSO, and this derived method is been called as CLPSO-LHS. A control method to set the size of a hypercube is also proposed in this work. At the beginning, the size of the hypercube is set as a relatively large value, which means that the hypercube for sampling is large such that the local search can improve the ability of global search. The size decreases with the iteration. In this way, the global search ability will be weakened and the local search ability will be strengthened.

Since cylindrical shapes are very common in precision components, cylindricity error evaluation problem is an important issue in assembly automation which is critical part of intelligent manufacturing. This problem belongs to minimax problems which are not differentiable and it is tough to be solved by traditional optimization methods. The least squares method is often used in most current commercial software. However, it is prone to over-estimation ^[31]. Recently, many metaheuristic search methods have been applied in this problem, such as improved Genetic algorithms (GA) ^[32], particle swarm optimization (PSO) ^[33], hybrid particle swarm optimization-differential evolution algorithm (PSO-DE) ^[31], and hierarchical PSO with Latin sampling based memetic algorithm (MA-HPSOL) ^[22]. The CLPSO-LHS will be applied in this problem and compared with the above-mentioned algorithms.

This paper is organized as follows: A detailed introduction of CLPSO-LHS is presented in Section 2. Then, benchmark testing for verifying the efficiency of the hybrid algorithm is shown in Section 3. Thereafter, the application of the hybrid algorithm in cylindricity error evaluation problem is presented in Section 4. And finally, Section 5 gives the conclusion of this work.

2.   The hybrid algorithm CLPSO with local search based on Latin hypercube sampling (CLPSO-LHS)

2.1.   PSO and comprehensive learning particle swarm optimizer (CLPSO)

In the original PSO, a particle swarm simulating a bird flock is used to search the whole space. The i^th particle in t^th iteration is denoted as ${X_{i, t}} = (X_{i, t}^1, X_{i, t}^1, ...X_{i, t}^D)$ and its corresponding flying velocity is denoted as ${V_{i, t}} = (V_{i, t}^1, V_{i, t}^1, ...V_{i, t}^D)$. The corresponding personal best position which memories the best position found by the i^th particle is denoted as $pbes{t_{i, t}} = (pbest_{i, t}^1, pbest_{i, t}^1, ...pbest_{i, t}^D)$. The global best position recording the best position discovered by the whole swarm is denoted as $gbes{t_t} = (pbest_t^1, pbest_t^1, ...pbest_t^D)$. Shi and Eberhart proposed the PSO-w 15 by introducing an inertia weight (w) into the original PSO to enhance the performance. Since it improves the performance significantly, many later researchers considered PSO-w as the standard PSO. The updating equations of velocity and position in PSO-w are shown as follows:

Where c₁ and c₂ are acceleration constants, r₁ and r₂ are random number generated by a uniform distribution in [0, 1].

The comprehensive learning PSO (CLPSO) 19 adopts a learning strategy in which all other particles' historical best information is used to update a particle's velocity. The strategy makes the particles own more exemplars to learn from and a larger potential space to traverse to. The velocity updating way is shown in formula (3):

Here, ${f_i} = [{f_i}(1), {f_i}(2), ..., {f_i}(D)]$ defines which particles' pbests the particle i should follow. A parameter Pc_i, referred to as the learning probability, is used to assign the value for f_i. If a random value (in the range [0, 1] in dimension d) is greater than Pc_i, i is assigned to f_i(d). Otherwise, a random integer from the particles' indexes is assigned to f_i(d). The Pc_i value is set by the following expression:

$pbest_{fi(d)}^d$ can be the corresponding dimension of any particle's pbest, including its own pbest.

2.2.   Local search based on Latin hypercube sampling

Only one sample is selected in each row or column of each sub hypercube in Latin hypercube sampling. The pseudocode of LHS is shown in Figure 1 ^[22]. Note that n is the dimension number of the original hypercube.

Figure 2 shows an example for sampling 5 points from a 2-dimensional square by LHS. It can be observed that only one sub hypercube is selected in each row or column.

2.3.   Control method for the size of a hypercube

When Latin hypercube sampling is used as a local search method, the size of a hypercube is a very important parameter which significantly impacts the efficiency of sampling. The parameter should decrease with the converge process. Based on the principle, the length vector of dimensions δ for a hypercube is set as formula (5). In formula (5), δ_max is the initial value which is set as a ratio of (UpB-LowB). UpB and LowB denote upper boundary and lower boundary, respectively. And $\tau = Itera/\max \_{\rm{gen}}$, where Itera is the current iteration number and max_gen is the maximum number of iterations. γ is a number which determines the magnitude order. p is a positive integer which is the polynomial order. k is a non-negative integer which represents the cycle for the sine function. Note that the value Itera is the only variable in the formula. The value δ varies with the current iteration number Itera. δ and δ_max are vectors in a practical problem since each dimension's UpB and LowB is different.

Formula (5) is composed of an exponential function, a sine function and the exponential part of the exponential function is polynomial. Actually, it is obtained by trial and error. Based on the principle that δ value should decrease with the converge process, δ value was set firstly as linear control formula (6). It works well for global search. However, it converges too fast at the end leading to bad local search. Exponential control formula (7) is designed to overcome the shortage of formula (6). Further, formula (8) was designed to remaining the global search ability. At last, a sine function was added in order to improve the converge speed and search accuracy in formula (5).

When δ_max is 1, γ is 10, p is 3, k is 6 and Max_gen is 5000, the plots of the four control methods are shown in Figure 3. It can be observed that the converge curve of the linear control way for formula (6) decreases to zero suddenly at the end of iteration process, which would lead to poor intensification capability, although the global search ability is strong at the beginning. The curve of the exponential control method for formula (7) decreases as a straight line. Actually, it decreases so fast at the beginning, such that it may not provide sufficient help for global search. In order to overcome this shortage, a polynomial exponential part for the exponential function is designed in formula (8). Based on formula (8), a sine function is introduced in formula (5) to improve the converge speed and search accuracy. Compared with formula (8), the δ value in formula (5) varies periodically. Since the value will reach a smaller value shortly, it can speed the convergence if the near-optimal solution has been found. At the end of iterations, the δ value can reach a smaller value, which will improve the search accuracy. Since the local search only are executed every ten iterations, τ can only get discrete values (e.g. 10/5000, 20/5000…). The value (k+0.5)(1-τ) cannot be pure integer. This is the reason why the plot of formula (5) cannot reach 0 in figure 3. In a word, formula (5) is the outcome of balance between global search ability and local search ability. These control methods will be compared in Section 3.2.

2.4.   Flow chart of LHS local search

The flow chart of LHS local search is shown in Figure 4. In CLPSO-LHS, the local search is executed for the gbest in CLPSO for every 10 iterations. The framework of CLPSO can be found from 19.

2.5.   Computational complexity of CLPSO-LHS

Computational complexity analysis is important for evaluating the efficiency of an algorithm. The count of key program statements within an iteration in big O notation is used to represent computational complexity here. The CLPSO-LHS is compared with the original CLPSO. Based on the analysis from 20, CLPSO and simple PSO share the same computational complexity which is O(ND) where N is population size and D is dimensionality. The local search in Figure 6 includes two parts. Since the complexity of Randperm(H) which is used for generating a random sequence of [1, …, H] is O(H), the complexity for the first part is O(HD). Clearly, the complexity for the second part is also O(HD). The complexity for the whole local search is O(HD). Assuming the local search is executed for the gbest for every G iterations (G is 10 in current work), the total complexity for the CLPSO-LHS is O(ND + HD/G). Usually, H/G is less than N (N is 10 or 20, H and G is 10 in the benchmark testing), so the final complexity is still O(ND). As a conclusion, the computational complexity of the proposed CLPSO-LHS is the same as the original CLPSO.

3.   Benchmark testing

In this section, eight benchmark functions with 10 and 30 dimensions in 19 are used to test the CLPSO with Latin hypercube sampling-based local search. The details of the test functions are attached in the Appendix.

Firstly, the proposed hybrid algorithm is compared with the original CLPSO and the pure Latin hypercube sampling-based local search in Section 3.1. The four control methods are compared in Section 3.2. The comparison with other algorithms in literature is shown in Section 3.3.

3.1.   Comparison with the original CLPSO and the Latin hypercube sampling-based local search

In order to prove the efficiency of the hybrid algorithm, firstly, it was compared with the original CLPSO and the Latin hypercube sampling-based local search. Note that the results for CLPSO here are different from those in 19, since that the population sizes are different. In order to take advantage of LHS, smaller population sizes have been used here to enlarge the Max_gen. The results in 19 will be compared with CLPSO-LHS in Section 3.3.

The three algorithms were run 30 times randomly on all functions with 10 and 30 dimensions. The mean values (Mean) and standard deviation (Std) of the results were recorded. The parameters of the algorithms are shown in Table 1. CLPSO and CLPSO-LHS shared some common parameters. The initialization method in 19 was used here. For the pure LHS local search, only one solution was generated at the beginning, and formula (5) was used to set the size of hypercube. The total sampling number was set as MaxNFE value. The parameters in formula (5) were set as the same as in CLPSO-LHS except δ_max. The value was set as a bigger value in LHS local search such that the algorithm can search the whole solution space.

From Table 2, it can be observed that the LHS local search found good solutions for 9 out from 16 functions such as F₁, F₃, F₄ and F₅. However, for the other functions, the LHS local search performs worse than the original CLPSO. It can be inferred that the LHS local search does have strong local search ability and it also has some degree of global search ability. Otherwise, it cannot find good solutions for F₃, F₄, and F₅ which are multimodal problems. However, its global search ability is limited, so it failed for some other functions and failed on more functions with 30 dimensions which requires stronger global search ability because the solution space expands exponentially with the dimensionality. Since CLPSO has strong global search ability, it performed better than LHS local search on these functions. Since the hybrid algorithm CLPSO-LHS has the global search ability of CLPSO and the strong local search ability of LHS local search, it should perform better than the former two algorithms. It did find best solutions for all functions. Even for F₂ which has a narrow valley from the perceived local optima to the global optimum 19, the results were much better than the original CLPSO and the pure LHS local search. For F₂ with 10 and 30 dimensions, CLPSO-LHS only failed five times and once to find the global optima, respectively. This is the reason why the standard deviation was large for F₂. In a conclusion, it is worthwhile to combine those two algorithms to form CLPSO-LHS.

3.2.   Comparison of the four control methods

Here, the CLPSO with Latin hypercube sampling-based local search using formula (6) to formula (8) have been denoted as CLPSO-LHS1, CLPSO-LHS2 and CLPSO-LHS3, respectively. The testing method and parameters are the same as those in Section 3.1.

The mean values (Mean) and standard deviation (Std) of the results for the four algorithms are shown in Table 3. From the table, several items can be observed as follows. 1) Among the four algorithms, the algorithms with exponential control methods (CLPSO-LHS, CLPSO-LHS2, CLPSO-LHS3) perform better than that with the linear control method (CLPSO-LHS1). They get the global optima value (0) for 9 out 16 functions. From Figure 3, the reason can be inferred that the linear control method provides relatively large δ values in the almost whole search progress. Large δ values contribute the global search. As a result, CLPSO-LHS1 found near-optimum solutions for most functions. 2) CLPSO-LHS2 and CLPSO-LHS3 are comparable. However, CLPSO-LHS3 performs better for F₂. The difference between CLPSO-LHS2 and CLPSO-LHS3 locates at the p value in formula (8). The p value is 3 for CLPSO-LHS3 and it is 1 for CLPSO-LHS2. The p value determines the decrease speed of δ values and controls the balance between global search ability and local search ability. Since F₂ is a hard problem which needs stronger global search ability, CLPSO-LHS3 with larger p value performs better. 3) If any other algorithms got global optimal value for a function, CLPSO-LHS also obtained its global optimal value. It performs best on those functions whose global optimal value is unable to reach. CLPSO-LHS takes formula (5) to set the δ value for LHS local search. Compared with CLPSO-LHS3, a sine function is added for two purposes. On one hand, the sine function can make the search converge in advance. On the other hand, it can provide smaller δ values which helps find more accuracy solutions. From Table 3, the second purpose was verified. The first purpose will be verified in convergence curves.

The convergence curves for the five algorithms (CLPSO is also included) for all of the functions are plotted in Figure 5 and Figure 6. Note that the iteration number is used as horizontal ordinate here. It is fair for comparison among the four CLPSO-LHS algorithms, because that they have the same number of function evaluations (NFE) at each iteration. However, the NFE in CLPSO is smaller than that in the four CLPSO-LHS algorithms. According to the parameters setting, the NFE in CLPSO approximately equals to ps/(ps + H/10) times of that in CLPSO-LHS at the same iteration. Since H value is 10, the difference between the NFE and ps/(ps + 1) *NFE is not big. Based on the above consideration, the iteration number is used as horizontal ordinate. By the way, since the MaxNFE is the same, the maximal iteration is larger for CLPSO.

From the two figures, it can be observed that:

1): The plots of the four hybrid algorithms are under the plots of CLSPO for most functions, especially at the beginning. Compared with CLPSO, LHS local search is added into the four hybrid algorithms. Since the local search is executed for the gbest particle and it has global search ability, the hybrid algorithms can find better solutions.

2): The plots of the four hybrid algorithms found good solutions for most functions, especially for CLPSO-LHS. The reason is that the LHS has good local search ability.

3): The plots of CLPSO-LHS1, CLPSO-LHS2 and CLPSO-LHS3 have similar trajectories with the plots for corresponding δ values in Figure 3 for some functions such as F₂, F₃. It can be inferred that when a near-optimum solution is found, the LHS local search works well and the search accuracy is determined by the δ value.

4): The plots of CLPSO-LHS2 converge fast for most functions except for F₂ which needs strong global search ability at the beginning. Since there is a narrow valley from the perceived local optima to the global optimum for F₂, the original CLPSO cannot find its good solution. With the help of LHS, the hybrid algorithms can find better solutions. However, if the δ value decreases too fast, the global search ability is not strong enough and the hybrid algorithm cannot find good solutions. This situation happens on CLPSO-LHS2.

5): The plots of CLPSO-LHS are staircase-like and they are under the plots of CLPSO-LHS3 for most functions which indicate that the sine function in formula (5) works well. The sine function really contributes to the convergence speed.

6): For most functions, CLPSO-LHS can reach better values at the beginning, and then, CLPSO-LHS2 exceeds it. However, CLPSO-LHS can find smaller values at the end. The sine function makes CLPSO-LHS converge fast at the beginning and it provides smaller δ values at the end leading to better final solutions.

Based on the testing results and plots, some conclusions can be made as follows: 1) The LHS local search not only contributes to the local search ability but also global search ability for CLPSO-LHS. 2) The control method shown in formula (5) can balance the global search and local search ability well and the sine function part can accelerate the convergence and improve the final results.

3.3.   Comparison with other variants of PSO

In this section, CLPSO-LHS is compared with other 12 variants of PSO from literature. The 12 algorithms are listed as follows.

● PSO with inertia weight (PSO-w) 15

● PSO with constriction factor (PSO-cf) ^[34]

● Local version of PSO with inertia weight (PSO-w-local) 19

● Local version of PSO with constriction factor (PSO-cf-local) 16

● UPSO 18

● Fully informed particle swarm (FIPS) 17

● FDR-PSO ^[35]

● CPSO-H ^[36]

● CLPSO 19

● CPSO-inner 20

● CPSO-outer 20

● PCLPSO 23

The first nine algorithms' results are obtained from 19. For 10-dimentional functions, the nine algorithms' population size is 10 and the maximum number of function evaluation is 30,000. For 30-dimentional functions, the population size is 40 and the maximum number of function evaluation (MaxNFE) is 200,000. CPSO-inner and CPSO-outer hybridized cellular automata in PSO. The results are from 20. The CPSO-inner is the hexagonal version which performed best among the three versions. In the two versions of CPSO, the population size is 36 and the maximum number of iterations is set as 5,000 for both 10-dimensional and 30-dimensional functions. The MaxNFE in CPSO-inner is 180,000 (36 × 5,000) and the value is 230,000 ((36 + 10) × 5,000) in CPSO-outer. PCLPSO 23 is a better variation of CLPSO. It uses a novel parallel multi-swarm strategy. The setting of population size and MaxNFE is the same as in the first nine algorithms. The parameters for CLPSO-LHS are set as in Table 1. The independent running number is 30 for all functions.

The results of all the algorithms on 10-dimensional and 30-dimensional functions are shown in Table 4 and Table 5. The best results for each function are stressed by bold.

From Table 4, 5, it can be observed that CLPSO-LHS performed best for 10 out from 16 functions. CPSO-outer had best results for 7 10-dimensional functions since its MaxNFE was much larger. PCLSPO performed best for 6 functions. CLPSO archived 5 best results. CPSO-inner archived best results on 3 functions and CPSO-H performs best on 2 functions. Benefiting from the LHS local search, CLPSO-LHS still was able to achieve good results when almost all the other algorithms failed for F₂. Only CLPSO, PCLPSO and CLPSO-LHS, three CLPSO based algorithms, archived good results for F₈. For 10-dimensional F₄, only CLPSO-LHS and CPSO-outer obtained good results. CLPSO-LHS still obtained good results for the other six functions where it didn't perform best. PCLPSO and CLPSO performed very well on most functions. CPSO-outer shows good performance in 10-dimensional functions with the larger MaxNFE, but it is unable to perform well on most 30-dimensional functions. From an overall perspective, among all the compared algorithms, CLPSO-LHS outperformed other algorithms, especially for 30-dimentional functions.

4.   Cylindricity error evaluation problem

The calculation of cylindricity error is illustrated in Figure 7. Assuming that the axis direction is n(m, n, 1) and the radius is R. The aim of the problem is to find an optimal axis direction n which passes the point Q(x₀, y₀, 1) to minimize the cylindricity error. The variants are m, n, x₀ and y₀. Given a point set P = {P_k|k = 1, 2, …, M}, assuming that the length from point P₁ to the axis is longest and the length from P₂ to the axis is shortest, the cylindricity error is,

If both |E₁P₁| and |E₂P₂| are bigger than R. When |E₁P₁| is bigger than R and |E₂P₂| is shorter than R, the cylindricity error is,

When both |E₁P₁| and |E₂P₂| are shorter than R, the same formula is obtained. So the objective function is written as formula (9).

Where

Four variables (x₀, y₀, l, m) were initialized in [−1, 1]. Two data sets from ^[33] and ^[37] were used to test the hybrid algorithm. The part for the first data set is with the dimension of 39 mm in radius and 120 mm in length. The measurement data for the cylindrical outer surface contains 80 records. The radius and the length of the part for the second data set is 60 mm and 30 mm, respectively. There are 20 data records. The data records are formatted as (x_k, y_k, z_k) which is the input for the objective function as formula (9).

The population size of the CLPSO-LHS was set to be 6 and the MaxNFE was set to be 40000. The other parameters were the same as those in Section 3.2. The algorithms were run 30 times independently on the two data sets. The statistical results (best, worst, median, mean value and the standard deviation) are shown as in Table 6.

Table 7 shows the comparison of best solutions between the CLPSO-LHS and other algorithms. The compared algorithms include improved GA1 ^[32], improved GA2 ^[38], PSO ^[33], MA-HPSOL ^[22], and PSO-DE ^[31].

CLPSO-LHS and MA-HPSOL perform best for the first data set. The best solution obtained by them is much better than the one obtained by PSO. For the second data set, CLPSO-LHS is the best one and it is much better than others including MA-HPSOL.

5.   Conclusions

In this paper, a local search based on Latin hypercube sampling is combined with CLPSO. The local search shows a good performance when a control method is used to control the size of the hypercube. The numerical experiments show that, at the beginning, the local search method contributes to global search ability. Ultimately, it provides good fine search ability. The algorithm performs well on testing functions. The hybrid algorithm also shows a good performance on its application in cylindricity error evaluation problems.

Although the algorithm is very effective, there is still some room for improvement. First, the designed control method has lots of parameters. A self-adaptive control method can be developed and its parameters can be reduced in the future. Second, although the local search method helps to find a good solution at the beginning, it can be observed from convergence figures that the solution cannot be improved for a short time for some functions, such as F₁, F₅ and F₈. How to utilize the global search ability to accelerate the converge speed deserves further study.

Acknowledgement

This research work is supported in part by National Natural Science Foundation of China (NSFC) under Grant No. 61603145, in part by Natural Science Foundation of Hubei Province under Grant No. 2015CFB528, and in part by Open Project of State Key Lab of Digital Manufacturing Equipment & Technology under Grant DMETKF2018011.

Conflicts of interest

All authors declare no conflict of interest in this paper.