A Novel Monarch Butterfly Optimization with Global Position Updating Operator for Large-Scale 0-1 Knapsack Problems

: As a significant subset of the family of discrete optimization problems, the 0-1 knapsack problem (0-1 KP) has received considerable attention among the relevant researchers. The monarch butterfly optimization (MBO) is a recent metaheuristic algorithm inspired by the migration behavior of monarch butterflies. The original MBO is proposed to solve continuous optimization problems. This paper presents a novel monarch butterfly optimization with a global position updating operator (GMBO), which can address 0-1 KP known as an NP-complete problem. The global position updating operator is incorporated to help all the monarch butterflies rapidly move towards the global best position. Moreover, a dichotomy encoding scheme is adopted to represent monarch butterflies for solving 0-1 KP. In addition, a specific two-stage repair operator is used to repair the infeasible solutions and further optimize the feasible solutions. Finally, Orthogonal Design (OD) is employed in order to find the most suitable parameters. Two sets of low-dimensional 0-1 KP instances and three kinds of 15 high-dimensional 0-1 KP instances are used to verify the ability of the proposed GMBO. An extensive comparative study of GMBO with five classical and two state-of-the-art algorithms is carried out. The experimental results clearly indicate that GMBO can achieve better solutions on almost all the 0-1 KP instances and significantly outperforms the rest.


Introduction
The 0-1 knapsack problem (0-1 KP) is a classical combinatorial optimization task and a challenging NP-complete problem as well.That is to say, it can be solved by nondeterministic algorithms in polynomial time.Similar to other NP-complete problems, such as vertex cover (VC), hamiltonian circuit (HC), and set cover (SC), the 0-1 KP is intractable.In other words, no polynomial-time exact algorithms have been found for it thus far.This problem was originated from the resource allocation involving financial constraints and since then, has been extensively studied in an array of scientific fields, such as combinatorial theory, computational complexity theory, applied mathematics, and computer science [1].Additionally, it has been found to have many practical applications, such as project selection [2], investment decision-making [3], and network interdiction problem [4].Mathematically, we can describe the 0-1 KP as follows: where n is the number of items, pi, and wi denote the profit and weight of item i, respectively.C represents the total capacity of the knapsack.The 0-1 variable xi indicates whether the item i is put into the knapsack, i.e., if any item i is selected and belongs to the knapsack, xi = 1, otherwise, xi = 0.The objective of the 0-1 KP is to maximize the total profits of the items placed in the knapsack, subject to the condition that the sum of the weights of the corresponding items does not exceed a given capacity C. Since the 0-1 KP was reported by Dantzig [5] in 1957, a large number of researchers have focused on addressing it in diverse areas.Some of the main early methods in this field are exact methods, such as the branch and bound method (B&B) [6] and the dynamic programming (DP) method [7].It is a breakthrough to introduce the concept of the core by Martello et al. [8].In addition, some effective algorithms have been proposed for 0-1 KP [9], the multidimensional knapsack problem (MKP) [10].With the rapid development of computational intelligence, some modern metaheuristic algorithms have been proposed for addressing the 0-1 KP.Some of those related algorithms include genetic algorithm (GA) [11], differential evolution (DE) [12], shuffled frog-leaping algorithm (SFLA) [13], cuckoo search (CS) [14,15], artificial bee colony (ABC) [16,17], harmony search (HS) [17][18][19][20][21], and bat algorithm (BA) [22,23].Many research methods are applied to the 0-1 KP problem.Zhang et al. converted the 0-1 KP problem into a directed graph by the network converting algorithm [24].Kong et al. proposed an ingenious binary operator to solve the 0-1 KP problem by simplified binary harmony search [20].Zhou et al. presented a complex-valued encoding scheme for the 0-1 KP problem [22].
As a novel biologically inspired computing approach, MBO is inspired by the migration behavior of the monarch butterflies with the change of the seasons.The related investigations [42,43] have demonstrated that the advantage of MBO lies in its simplicity, being easy to carry out, and efficiency.In order to address the 0-1 KP, which falls within the domain of the discrete combinatorial optimization problems with constraints, this paper presents a specially designed monarch butterfly optimization with global position updating operator (GMBO).What needs special mention is that GMBO is a supplement and perfection to previous related work, namely, a binary monarch butterfly optimization (BMBO) and a novel chaotic MBO with Gaussian mutation (CMBO) [42].The main difference and contributions of this paper are as follows, compared with BMBO and CMBO.
Firstly, the original MBO was proposed to address the continuous optimization problems, i.e., it cannot be directly applied in the discrete space.For this reason, in this paper, a dichotomy encoding strategy [44] was employed.More specifically, each monarch butterfly individual is represented as two-tuples consisting of a real-valued vector and a binary vector.Secondly, although BMBO demonstrated excellent performance in solving 0-1 KP, it did not show a prominent advantage [42].In other words, some techniques can be combined with BMBO for the purpose of improving its global optimization ability.Based on this, an efficient global position updating operator [16] was introduced to enhance the optimization ability and ensure its rapid convergence.Thirdly, a novel two-stage repair operator [45,46] called the greedy modification operator (GMO), and greedy optimization operator (GOO), respectively, was adopted.The former repairs the infeasible solutions while the latter optimizes the feasible solutions during the search process.Fourthly, empirical studies reveal that evolutionary algorithms have certain dependencies on the selection of parameters.Moreover, certain coupling between the parameters still exists.However, suitable parameter combination for a particular problem was not analyzed in BMBO and CMBO.In order to verify the influence degree of four important parameters on the performance of GMBO, an orthogonal design (OD) [47] was applied, and then the appropriate parameter settings were examined and recommended.Fifthly, generally speaking, the approximate solution of an NP-hard problem can be obtained by evolutionary algorithms.However, the most important thing is to obtain higher quality approximate solutions, which are closer to the optimal solutions more profitably.In BMBO, the optimal solutions of all the 0-1 KP instances were not provided.It is difficult to judge the quality of an approximate solution obtained by an evolutionary algorithm.In GMBO, the optimal solutions of 0-1 KP instances are calculated by a dynamic programming algorithm.Meanwhile, the approximation ratio based on the best values and the worst values are provided, which clearly reflect the degree of the closeness of the approximate solutions to the optimal solutions.In addition, the application of statistical methods in GMBO is one of the differences between GMBO and BMBO, CMBO, including Wilcoxon's rank-sum tests [48] with a 5% significance level.Moreover, boxplots can visualize the experimental results from the statistical perspective.
The rest of the paper is organized as follows.Section 2 presents a snapshot of the original MBO, while Section 3 introduces the GMBO for large-scale 0-1 KP in detail.Section 4 reports the outcomes of a series of simulation experiments as well as to compare results.Finally, the paper ends with Section 5 after providing some conclusions, along with some directions for further work.

Monarch Butterfly Optimization
Animal migration involves mainly long-distance movement, usually in groups, on a regular seasonal basis.MBO [33,43] is a population-based intelligent stochastic optimization algorithm that mimics the seasonal migration behavior of monarch butterflies in nature.It should be noted that the entire population is divided into two subpopulations, named subpopulation_1 and subpopulation_2, respectively.Based on this, the optimization process consists of two operators, which operate on subpopulation_1 and subpopulation_2, respectively.The information is interchanged among the individuals of subpopulation_1 and subpopulation_2 by applying the migration operator.The butterfly adjusting operator delivers the information of the best individual to the next generation.Additionally, Lévy flights [49,50] are introduced into MBO.The main steps of MBO are outlined as follows: Step 1. Initialize the parameters of MBO.There are five basic parameters to be considered while addressing various optimization problems, including the number of the population (NP), the ratio of the number of monarch butterflies in subpopulation_1 (p), migration period (peri), the monarch butterfly adjusting rate (BAR), the max walk step of the Lévy flights (Smax).
Step 2. Initialize the population with NP randomly generated individuals according to a uniform distribution in the search space.Step 3. Sort the individuals according to their fitness in descending order (Here assumptions for the maximum).The better NP1 (p*NP) individuals constitute subpopulation_1, and NP2 (NP-NP1) individuals make up subpopulation_2.
Step 4. The position updating of individuals in subpopulation_1 is determined by the migration operator.The specific procedure is described in Algorithm 1.
Step 5.The moving direction of the individuals in subpopulation_2 depends on the butterfly adjusting operator.The detailed procedure is shown in Algorithm 2.

Step 6. Recombine two subpopulations into one population
Step 7. If the termination criterion is already satisfied, output the best solution found, otherwise, go to Step 3. where dx is calculated by implementing the Lévy flights.It should be noted that the Lévy flights, which originated from the Lévy distribution, are an impactful random walk model, especially on undiscovered, higher-dimensional search space.The step size of Lévy flights refers to Equation (2).

A Novel Monarch Butterfly Optimization with Global Position Updating Operator for the 0-1 KP
In this section, we give the detailed design procedure of the GMBO for the 0-1 KP.Firstly, a dichotomy encoding scheme [46] is used to represent each individual.Secondly, a global position updating operator [16] is embedded in GMBO in order to increase the probability of finding the optimal solution.Thirdly, the two-stage individual optimization method is employed, which successively tackles the infeasible solutions and then further improves the existing feasible solutions.Finally, the basic framework of GMBO for 0-1 KP is formed.

Dichotomy Encoding Scheme
KP belongs to the category of discrete optimization.The solution space is a collection of discrete points rather than a contiguous area.For this reason, we should either redefine the evolutionary operation of MBO or directly apply a continuous algorithm to discrete problems.In this paper, we prefer the latter for its simplicity of operation, comprehensibility, and generality.
As previously mentioned, each monarch butterfly individual is expressed as a two-tuple <X, Y>.Here, real number vectors X still constitute the search space as in the original MBO, which can be regarded as a phenotype similar to the genetic algorithm.Binary vectors, Y, form the solution space, which can be seen as a genotype common in the evolutionary algorithm.It should be noted that Y may be a valid solution because 0-1KP is a constraint optimization problem.Here we abbreviate the monarch butterfly population to MBOP.Then the structure of MBOP is given as follows: 1,1 The first step to adopting a dichotomous encoding scheme is to transfer the phenotype to genotype.Therefore, a surjective function g is used to realize the mapping relationship from each element of X to the corresponding element of Y. where is the sigmoid function.The sigmoid function is often used as the threshold function of neural networks.It was applied to the binary particle swarm optimization (BPSO) [51] to convert the position of a particle from a real-valued vector to a 0-1 vector.It should be noted that there are other conversion functions [52] can be used.Now assume a 0-1 KP problem with 10 items, Figure 1 shows the above process, in which each [ ] -5.0,5.0 ) is randomly chosen based on the uniform distribution.

Global Position Updating Operator
The main feature of particle swarm optimization (PSO) is that the particle always tends to converge to two extreme positions viz. the best position ever found by itself and the global best position.Inspired by the behavior of swarm intelligence of PSO, a novel position updating operator was recently proposed and successfully embedded in HS for solving 0-1 KP [16].After that, the position updating operator combines with CS [14] to deal with 0-1 KP.
It is well-known that the evolutionary algorithm can yield strong optimization performance under the condition of the balance between exploitation and exploration, or attraction and diffusion [53].The original MBO concentrates much on the exploration ability but weak exploitation capability [33,43].With the aim of enhancing the exploitation capability of MBO, we introduce a global position updating operator mentioned above.The procedure is shown in Algorithm 3, where "best" and "worst" represent the global best individual and the global worst individual, respectively.r4, r5, and rand are uniform random real numbers in [0, 1].The pm parameter is mutation probability.

Begin
for i=1 to NP (for all monarch butterflies in the whole population) for j=1 to d (all the elements in ith monarch butterfly) where rand ~ U(0,1) end if end for j end for i End.

Two-Stage Individual Optimization Method Based on Greedy Strategy
Since KP is a constrained optimization problem, it may lead to the occurrence of infeasible solutions.There are usually two major methods: Redefining the objective function by penalty function method (PFM) [54,55] and individual optimization method based on the greedy strategy (IOM) [56,57].Unfortunately, the former shows poor performance when encountering large-scale KP problems.In this paper, we adopt IOM to address infeasible solutions.
A simple greedy strategy, namely GS [58], is proposed to choose the item with the greatest density pi/wi first.Although the feasibility of all individuals can be guaranteed, it is obvious that there are several imperfections.Firstly, for a feasible individual, there is a possibility that the corresponding objective function value may turn to be worse by applying GS.Secondly, the lack of further optimization for all individuals can lead to unsatisfactory solutions.
In order to overcome the shortcomings of GS, the two-stage individual optimization method is proposed by He et al. [45,46].A greedy modification operator (GMO) is used to repair the infeasible individuals in the first stage.It is followed by the application of the greedy optimization operator (GOO), which further optimizes the feasible individuals.The method proceeds as follows.
Step 1. Quicksort algorithm is used to sort all items in the non-ascending order according to pi/wi, and the index of items is stored in an array H [1], H [2]..., H[n], respectively.Step 2. For an infeasible individual , GMO is applied.
Step 3.For a feasible individual , GOO is performed.
After the above repair process, it is easy to verify that each optimized individual is feasible.The significance of GMO and GOO seems particularly prominent while solving high dimensional KP problems [45,46].The pseudo-code of GMO and GOO can be shown in Algorithms 4 and 5, respectively.

Begin
Input: End.

The Procedure of GMBO for the 0-1 KP
In this section, the procedure of GMBO for 0-1 KP is described in Algorithm 6, and the flowchart is illustrated in Figure 2. Apart from the initialization, it is divided into three main processes.
(1) In the migration process, the position of each monarch butterfly individual in subpopulation_1 is updated.We can view this process as exploitation by combining the properties of the currently known individuals in subpopulation_1 or subpopulation_2.
(2) In the butterfly adjusting process, partial genes of the global best individual are passed on to the next generation.Moreover, Lévy flights come into play owing to longer step length in exploring the search space.This process can be considered as exploration, which may find new solutions in the unknown domain of the search space.
(3) In the global position updating process, we can define the distance of the global best individual and the global worst individual as the adaptive step.Obviously, the two extreme individuals differ greatly at the early stage of the optimization process.In other words, the adaptive step has a larger value, and the search scope is broader, which is beneficial to the global search over a wide range.With the progress of the evolution, the global worst individual tends to be more similar to the global best individual, and then the difference becomes small at the late stage of the optimization process.Meanwhile, the adaptive step has a smaller value, and the search area narrows, which is useful for performing the local search.In addition, the genetic mutation is applied to preserve the population diversity and avoid premature convergence.It should be noted that, unlike the original MBO, in GMBO, the two newly-generated subpopulations regroup one population at a certain generation rather than each generation, which can reduce time consumption.

Begin
Step 1: Sorting.Quicksort is used to sort all items in the non-ascending order by i i p w ,

The Time Complexity
In this subsection, the time complexity of GMBO is simply estimated (Algorithm 6).It is not hard to see that the time complexity of GMBO mainly hinges on steps 1-3.In

Simulation Experiments
We chose 3 different sets of 0-1 KP test instances to verify the feasibility and effectiveness of the proposed GMBO method.The test set 1 and test set 2 were widely used low-dimensional benchmark instances with dimension = 4 to 24.The test set 3 consisted of 15 high-dimensional 0-1 KP test instances generated randomly with dimension = 800 to 2000.

Experimental Data Set
The generation form of test set 3 was firstly given.Since the difficulty of the knapsack problems was greatly affected by the correlation between the profits and weights [59], 3 typical large scale 0-1 KP instances were randomly generated to demonstrate the performance of the proposed algorithm.
Here function Rand (a, b) returned a random integer uniformly distributed in [a, b].For each instance, the maximum capacity of the knapsack equaled 0.75 times of the total weights.The procedure is as follows:  Uncorrelated instances: In this section, 3 groups of large scale 0-1 KP instances with dimensionality varying from 800 to 2000 were considered.These 15 instances included 5 uncorrelated instances, 5 weakly correlated instances, and 5 strongly correlated instances.The dimension size was 800, 1000, 1200, 1500, and 2000, respectively.We simply denoted these instances by KP1-KP15.

Parameter Settings
As mentioned earlier, GMBO included 4 important parameters: p, peri, BAR, and Smax.In order to examine the effect of the parameters on the performance of GMBO, Orthogonal Design (OD) [47] was applied with uncorrelated 1000-dimensional 0-1 KP instance.Our experiment contained 4 factors, 4 levels per factor, and 16 combinations of levels.The combinations of different parameter values are given in Table 1.
For each experiment, the average value of the total profits was obtained with 50 independent runs.The results are listed in Table 2.  Using data from Table 2, we can carry out factor analysis, rank the 4 parameters according to the degree of influence on the performance of GMBO, and deduce the better level of each factor.The factor analysis results are recorded in Table 3, and the changing trends of all the factor levels are shown in Figure 3.
As we can see from Table 3 and Figure 3, p is the most important parameter and needs a reasonable selection for the 0-1 KP problems.A small p signifies more elements from subpopulation_2.Conversely, more elements were selected from subpopulation_1.For peri, the curve was in a small range in an upward trend.This implied individual elements from subpopulation_2 had more chance to embody in the newly generated monarch butterfly.For BAR and Smax, it can be seen from Figure 3 that the effect on the algorithm was not obvious.
According to the above analysis based on OD, the most suitable parameter combination is p2 = 3/12, peri4 = 1.4,BAR1 = 1/12, and Smax3 = 1.0, which will be adopted in the following experiments.

The Comparisons of the GMBO and the Classical Algorithms
In order to investigate the ability of GMBO to find the optimal solutions and to test the convergence speed, 5 representative classical optimization algorithms, including the BMBO [42], ABC [60], CS [61], DE [62] and GA [11], were selected for comparison.GA was an important branch in the field of computational intelligence that has been intensively studied since it was developed by John Holland et al.In addition, GA was representative of the population-based algorithm.DE was a vector-based evolutionary algorithm, and more and more researchers have paid attention to it since it was first proposed.Since then, it has been applied to solve many optimization problems.CS is one of the latest swarm intelligence algorithms.The similarity of CS with GMBO lies in the introduction of Levy flights.ABC is a relatively novel bio-inspired computing method and has the outstanding advantage of the global and local search in each evolution.
There are several points to explain.Firstly, all 5 comparative methods (not including GA) used the previously mentioned dichotomy encoding mechanism.Secondly, all 6 comparative methods used GMO and GOO to carry out the additional repairing and optimization operations.Thirdly, ABC, CS, DE, GA, MBO, and GMBO were short for 6 methods based on binary, respectively.
The parameters were set as follows.For ABC, the number of food sources is set to 25 and maximum search times = 100.CS, DE, and GA are set the same parameters as that of [15].For MBO, we take the same parameters suggested in Section 4.2.In addition, the 2 subpopulations recombined every 50 generations.GMBO and MBO have identical parameters except for that mutation probability pm = 0.25 is included in GMBO.
For the sake of fairness, the population sizes of six methods are set to 50.The maximum run time is set to 8 seconds for 800, 1000, and 1200 dimensional instances but 10 seconds for 1500 and 2000 dimensional instances.50 independent runs are performed to achieve the experimental results.
We use C++ as the programming language and run the codes on a PC with Intel (R) Core(TM) i5-2415M CPU, 2.30GHz, 4GB RAM.

The Experimental Results of GMBO on Solving Two Sets of Low-Dimensional 0-1 Knapsack Problems
In this section, 2 sets of 0-1 KP test instances were considered for testing the efficiency of the GMBO.The maximum number of iterations was set to 50.As mentioned earlier, 50 independent runs were made.The first set, which contained 10 low-dimensional 0-1 knapsack problems [19,20], was adopted with the aim of investigating the basic performance of the GMBO.The standard 10 0-1 KP test instances were studied by many researchers, and detailed information about these instances can be taken from the literature [13,19,20].Their basic parameters are recorded in Table 4.The experimental results obtained by GMBO are listed in Table 5.
Here, "Dim" is the dimension size of test problems; Opt.value is the optimal value obtained by DP method [7]; Opt.solution is the optimal solution; "SR" is success rate; "Time" is the average time to reach the optimal solution among 50 runs; "MinIter", "MaxIter" and "MeanIter" represents the minimum iterations, maximum iterations, and the average iterations to reach the optimal solution among 50 runs, respectively."Best", "Worst", "Mean" and "Std" are the best value, worst value, mean value, and the standard deviation, respectively.(1,0,0,1,0,0,0) f8 23 9767 (1,1,1,1,1,1,1,1,0,0,1,0,0,0,0,1,1,0,0,0,0,0,0) As can be seen from Table 5, GMBO can achieve the optimal solution for all 10 instances with 100% success rates.Furthermore, the best value, the worst value, and the mean value are all equal to the optimal value for every test problem.Obviously, the efficiency of GMBO is very high for the considered set of instances because GMBO can get the optimal solution in a very short time.The minimum iterations are only 1, and the mean iterations are less than 6 for all the test problems.In particular, for f6, HS [18], HIS [63], and NGHS [19] can only achieve the best value 50 while GMBO can get the optimal value 52.
The second set, which includes 25 0-1 KP instances, was taken from references [64,65].For all we know, the optimal value and the optimal solution of each instance are provided for the first time in this paper.The primary parameters are recorded in Table 6.The experimental results are summarized in Table 7.Compared to Table 5 above, three new evaluation criteria, that is "ARB", "ARW", and "ARM", are used to evaluate the proposed method."Opt.value"represents the optimal solution value obtained by the DP method.Here, the following definitions are given: Opt.value ARW= Worst Opt.value ARM= Mean (10) Here, "ARB" represents the approximate ratio [66] of the optimal solution value (Opt.value) to the best approximate solution value (Best).Similarly, "ARW" and "ARM" are based on the worst approximate solution value (Worst) and the mean approximate solution value (mean), respectively.ARB, ARW, and ARM indicate the proximity of Best, Worst, and Mean to the Opt.value, respectively.Plainly, ARB, ARW, and ARM are real numbers greater than or equal to 1.0.
Thus, the conclusion is that GMBO had superior performance in solving low-dimensional 0-1 KP instances.In this section, in order to make a comprehensive investigation on the optimization ability of the proposed GMBO, test set 3, which included 5 uncorrelated, 5 weakly correlated, and 5 strongly correlated large-scale 0-1 KP instances, were considered.The experimental results are listed in Tables 8-10 below.The best results on all the statistical criteria of each 0-1 KP instances, i.e., the best values, the mean values, the worst values, the standard deviation, and the approximation ratio, appear in bold.As noted earlier, Opt and Time represent the optimal value and time spending taken by the DP method, respectively.
The performance comparisons of the six methods on the five large-scale uncorrelated 0-1 KP instances are listed in Table 8.It can be seen that GMBO outperformed the other five algorithms on the six and five evaluation criteria for KP1 and KP4, respectively.In addition, GMBO obtained the best values concerning the best and the mean value for KP3 and was superior to the other five algorithms in the worst value for KP2.Unfortunately, GMBO failed to come up with superior performance while encountering 2000-dimensional 0-1 KP instances (KP5).MBO beat the competitors on KP5.Moreover, an apparent phenomenon can be observed, which points out that ABC has better stability.The best value of KP2 was achieved by CS.Obviously, DE and GA showed the worst performance for KP1-KP5.Meanwhile, the approximation ratio of the best value of GMBO for KP1 equaled approximately 1.0.Additionally, there was little difference between the worst approximation ratio of the best value (1.0242) of GMBO and the best approximation ratio of the best value (1.0237) of MBO for KP5.Table 9 records the comparison of the performances of six methods on five large-scale weakly correlated 0-1 KP instances.The experimental results in Table 9 differ from that in Table 8.It is clear that GMBO had a striking advantage in almost all the six statistical standards for KP6-KP9.For KP10, similarly to KP5, GMBO was still not able to win out over MBO.It is worth mentioning that the approximation ratio of the best value of GMBO for KP6-KP7, and KP9 equaled 1.0.Moreover, the standard deviation value of KP6-KP7 and KP9 obtained by GMBO was much smaller than the corresponding value of the other five algorithms.
A comparative study of the six methods on five large-scale strongly correlated 0-1 KP instances are recorded in Table 10.Obviously, GMBO outperforms the other five methods for KP11-KP14 on five statistical standards except for Std.ABC obtains the best Std values for KP11-KP15.To KP15, GMBO can get better values on the worst.CS, DE, and GA fail to show outstanding performance for this case.Under these circumstances, the approximation ratio of the worst value of GMBO for KP11-KP15 was less than 1.0019.For a clearer and more intuitive measure of the similar level of the theoretical optimal value and the actual value obtained by each algorithm, the values of ARB on three types of 0-1 KP instances are illustrated in Figures 4-6.From Figure 4, the ARB of GMBO for KP1, KP3, and KP4 were extremely close to or equal to 1. GMBO had the smallest ARB for KP1, KP3-KP5, except for KP2, for which CS obtained the smallest ARB.Similar to Figure 4, from Figure 5, GMBO still had the smallest ARB values, which are 1.0 (KP6, KP7, and KP9) or less than 1.015 (KP8, KP10).In terms of the strongly correlated 0-1 KP instances, GMBO consistently outperformed the other five methods (see Figure 6), in which GMBO had the smallest ARB values except for KP15.Particularly, the ARB of GMBO was even less than 1.0015 for KP15.
Overall, Tables 8-10 and Figures 4--6 indicate that GMBO was superior to the other five methods when addressing large-scale 0-1 KP problems.In addition, if we look at the worst values achieved by GMBO and the best values obtained by other methods, we can observe that for the majority instances, the former were even far better than the latter.
In order to test the differences between the proposed GMBO and the other five methods, Wilcoxon's rank-sum tests with the 5% significance level were used.Table 11 records the results of rank-sum tests for KP1-KP15.In Table 11, "1" indicates that GMBO outperforms other methods at 95% confidence.Conversely, "-1".Particularly, "0" represents that the two compared methods possess similar performance.The last three rows summarized the times that GMBO performed better than, similar to, and worse than the corresponding algorithm among 50 runs.

GMBO ABC CS DE GA MBO
Table 14 shows the statistical results of the six methods based on the worst values.The ranking order of the six methods was GMBO (1.20), MBO (2.07),ABC (2.60), CS (3.93), DE (5), and GA (6), which was identical with that in Table 12.
Then, a comparison of the six highest dimensional 0-1 KP instances, i.e., KP4, KP5, KP9, KP10, KP14, and KP15, is illustrated in Figures 7-12, which was based on the best profits achieved by 50 runs.From Figure 7, it can be easily seen that the best values obtained by GMBO far exceed that of the other five methods.Meanwhile, the two best values of CS outstripped the two worst values of GMBO.By looking at Figure 9, we can conclude that GMBO greatly outperformed the other five methods.The distribution of best values of GMBO in 50 times was close to a horizontal line, which pointed towards the excellent stability of GMBO in this case.With regard to numerical stability, CS had the worst performance.From Figure 11, the curve of GMBO still overtopped that of ABC, CS, DE, and GA, as illustrated in Figures 7 and 9.This advantage, however, was not obvious when compared with MBO.Figures 8,10,and 12 show the best values obtained by six methods on 2000-dimensional uncorrelated, weakly correlated, and strongly correlated 0-1 KP instances in 50 runs, respectively.As the dimension becomes large, space is expanded dramatically to 2 2000 , which represents a challenge for any method.It can be said with certainty that almost all the values of GMBO are bigger than that of the other five methods except MBO.Similar to Figure 11, the curves of MBO partially overlaps that of GMBO in Figure 12, which may be interpreted as the ability of GMBO towards competing with MBO.For the purpose of visualizing the experimental results from the statistical perspective, the corresponding boxplots of six higher dimensional KP4-KP5, KP9-KP10, and KP14-15 are shown in Figures 13-18.On the whole, the boxplot for GMBO has greater value and less height than those of the other five methods, which indicates the stronger optimization ability and stability of GMBO even encountering high-dimensional instances.In order to examine the convergence rate of GMBO, the evolutionary process and convergent trajectories of six methods are illustrated in Figures 19-24.It should be noted that six high dimensional instances, viz., KP4, KP5, KP9, KP10, KP14, and KP15, were chosen.In addition, Figures 19-24 show the average best values with 50 runs, and not one independent experimental result.From Figure 19, the curves of GMBO and MBO were almost coincident before 6 seconds, but afterward, GMBO converged rapidly to a better value as compared to the others.From Figure 20, it is indeed interesting to note that MBO has a weak advantage in the average values as compared to GMBO.From Figure 21, MBO and GMBO have identical initial function values, and the average values obtained by MBO were better than that of GMBO before 3 seconds.However, similar to the trend in Figure 19, 3 seconds later, GMBO quickly converged to a higher value.As depicted in Figure 22, unexpectedly, when addressing the 2000-dimensional weakly correlated 0-1 KP instances, GMBO was inferior to MBO.Figures 23-24 illustrate the evolutionary process of strongly correlated 0-1 KP instances.By observation of 2 convergence graphs, we can conclude that GMBO and MBO have similar performance.Throughout Figures 19-24, GMBO has a stronger optimization ability and faster convergence speed to reach optimum solutions than the other five methods.

The Comparisons of the GMBO and the Latest Algorithms
To evaluate the performance of the proposed GMBO, the two latest algorithms, namely, moth search (MS) [67] and moth-flame optimization (MFO) [68], were especially selected to compare with GMBO.The following factors were mainly considered.(1) The literature on the application of MS and MFO to solve 0-1 KP problem was not found.(2) The GMBO, MS, and MFO were novel nature-inspired swarm intelligence algorithms, which simulated the migration behavior of the monarch butterfly, the Lévy flight mode, or the navigation method of moths.
For MS, the max step Smax = 1.0, acceleration factor ϕ = 0.618, and the index β = 1.5.For MFO, the maximum number of flames N = 30.In order to make a fair comparison, all experiments were conducted in the same experimental environment as described above.The detailed experimental results of GMBO and the other two algorithms on the three kinds of large-scale 0-1 KP instances were presented in Table 15.The best, mean, and standard deviation values in bold indicate superiority.The dominant times of the three algorithms in the three statistical values are given in the last line of Table 15.As the results presented in Table 15, the number of times the GMBO algorithm had priority in the best, the mean, and the standard deviation values were 8, 10, 5, respectively.The simulation results indicated that GMBO generally provided very excellent performance in most instances compared with MFO and MS.The two metrics, namely mean and standard deviation, demonstrated again that GMBO was more stable.The comprehensive performance of MFO was superior to that of MS.To summarize, by analyzing Tables 4-15 and Figures 4-24, it can be inferred that GMBO had better optimization capability, numerical stability, and higher convergence speed.In other words, it can be claimed that GMBO is an excellent MBO variant, which is capable of addressing large-scale 0-1 KP instances.

Conclusion
In order to tackle high-dimensional 0-1 KP problems more efficiently and effectively, as well as to overcome the shortcomings of the original MBO simultaneously, a novel monarch butterfly optimization with the global position updating operator (GMBO) has been proposed in this manuscript.Firstly, a simple and effective dichotomy encoding scheme, without changing the evolutionary formula, is used.Moreover, an ingenious global position updating operator is introduced with the intention of enhancing the optimization capacity and convergence speed.The inspiration behind the new operator lies in creating a balance between intensification and diversification, a very important feature in the field of metaheuristics.Furthermore, a two-stage individual optimization method based on the greedy strategy is employed, which besides guaranteeing the feasibility of the solutions, is able to improve the quality further.In addition, the Orthogonal Design (OD) was applied to find suitable parameters.Finally, GMBO was verified and compared with ABC, CS, DE, GA, and MBO on large-scale 0-1 KP instances.The experimental results demonstrate that GMBO outperforms the other five algorithms on solution precision, convergence speed, and numerical stability.
The introduction of a global position operator coupled with an efficient two-stage repairing operator is instrumental towards the superior performance of GMBO.However, there is room for further enhancing the performance of GMBO.Firstly, the hybridization of the two methods complementing each other is becoming more and more popular, such as the hybridization of HS with CS [69].Combining MBO with other methods could indeed be very promising and hence worth experimentation.Secondly, in the present work, three groups of high-dimensional 0-1 KP instances were selected.In the future, a multidimensional knapsack problem, quadratic knapsack problem, knapsack sharing problem, and randomized time-varying knapsack problem can be considered to investigate the performance of MBO.Thirdly, some typical combinatorial optimization problems, such as job scheduling problems [70][71][72], feature selection [73][74][75], and classification [76], deserve serious investigation and discussion.For these challenging engineering problems, the key issue is how to encode and process constraints.The application of MBO for these problems is another interesting research area.Finally, perturb [77], ensemble [78], learning mechanisms [79,80], or information feedback mechanisms [81] can be effectively combined with MBO to improve performance.

Figure 1 .
Figure 1.The example of the dichotomy encoding scheme.
Step 1, Quicksort algorithm costs time O (n log n).In Step 2, the initialization of NP individuals consumes time O (NP * n).The calculation of fitness has time O (NP).In Step 3, migration operator costs time O (NP1 * n), and the butterfly adjusting operator has time O (NP2 * n).Moreover, the global position updating operator consumes O (NP * n).It is noticed that GMO and GOO cost the same time complexity O (NP * n).Thus, the time complexity of GMBO would be O (n log n

Figure 3 .
Figure 3.The changing trends of all the factor levels.

Figure 7 .
Figure 7.Comparison of the best values on KP4 in 50 runs.

Figures 7 , 9 ,
Figures 7,9, and 11 illustrate the best values achieved by the six methods on 1500-dimensional uncorrelated, weakly correlated, and strongly correlated 0-1 KP instances in 50 runs, respectively.From Figure7, it can be easily seen that the best values obtained by GMBO far exceed that of the other five methods.Meanwhile, the two best values of CS outstripped the two worst values of GMBO.By looking at Figure9, we can conclude that GMBO greatly outperformed the other five methods.The distribution of best values of GMBO in 50 times was close to a horizontal line, which pointed towards the excellent stability of GMBO in this case.With regard to numerical stability, CS had the worst performance.From Figure11, the curve of GMBO still overtopped that of ABC, CS, DE, and GA, as illustrated in Figures7 and 9.This advantage, however, was not obvious when compared with MBO.Figures 8,10,and 12 show the best values obtained by six methods on 2000-dimensional uncorrelated, weakly correlated, and strongly correlated 0-1 KP instances in 50 runs, respectively.As the dimension becomes large, space is expanded dramatically to 2 2000 , which represents a challenge for any method.It can be said with certainty that almost all the values of GMBO are bigger than that of the other five methods except MBO.Similar to Figure11, the curves of MBO partially overlaps that of GMBO in Figure12, which may be interpreted as the ability of GMBO towards competing with MBO.

Figure 8 .
Figure 8.Comparison of the best values on KP5 in 50 runs.

Figure 9 .
Figure 9.Comparison of the best values on KP9 in 50 runs.

Figure 10 .
Figure 10.Comparison of the best values on KP10 in 50 runs.

Figure 11 .
Figure 11.Comparison of the best values on KP14 in 50 runs.

Figure 12 .
Figure 12.Comparison of the best values on KP15 in 50 runs.

Figure 13 .
Figure 13.Boxplot of the best values on KP4 in 50 runs.

Figure 14 .
Figure 14.Boxplot of the best values on KP5 in 50 runs.

Figure 15 .
Figure 15.Boxplot of the best values on KP9 in 50 runs.

Figure 16 .
Figure 16.Boxplot of the best values on KP10 in 50 runs.

Figure 17 .
Figure 17.Boxplot of the best values on KP14 in 50 runs.

Figure 18 .
Figure 18.Boxplot of the best values on KP15 in 50 runs.

Figure 19 .
Figure 19.The convergence graph of six methods on KP4 in 10 seconds.

Figure 20 .
Figure 20.The convergence graph of six methods on KP5 in 10 seconds.

Figure 21 .
Figure 21.The convergence graph of six methods on KP9 in 10 seconds.

Figure 22 .
Figure 22.The convergence graph of six methods on KP10 in 10 seconds.

Figure 23 .
Figure 23.The convergence graph of six methods on KP14 in 10 seconds.

Figure 24 .
Figure 24.The convergence graph of six methods on KP15 in 10 seconds.

end if end if end for j end for i End.
Set the generation counter g = 1; set migration period peri, the migration ratio p, butterfly adjusting rate BAR, and the max walk step of Lévy flights Smax; set the maximum generation MaxGen.Set the generation interval between recombination RG.Generate NP monarch butterfly individuals randomly {<X1, Y1>, <X2, Y2>, …, <XNP, YNP>}.Calculate the fitness of each individual, f(Yi), 1≤i≤NP.
and the index of items is stored in array H[1...n].Step 2: Initialization.

Table 1 .
Combinations of different parameter values.

Table 2 .
Orthogonal array and the experimental results.

Table 3 .
Factor analysis with the orthogonal design (OD) method.

Table 4 .
The basic information of 10 standard low-dimensional 0-1KP instances.

Table 5 .
The experimental results of 10 standard low-dimensional 0-1 KP instances obtained by GMBO.

Table 10 .
Performance comparison of large-scale five strongly correlated 0-1KP instances.

Table 11 .
Rank sum tests for GMBO with the other five methods on KP1-KP15.

Table 15 .
Performance comparison of three algorithms on large-scale 0-1 KP instances.