XOR Binary Gravitational Search Algorithm with Repository: Industry 4.0 Applications

: Industry 4.0 is the fourth generation of industry which will theoretically revolutionize manufacturing methods through the integration of machine learning and artiﬁcial intelligence approaches on the factory ﬂoor to obtain robustness and speed-up process changes. In particular, the use of the digital twin in a manufacturing environment makes it possible to test such approaches in a timely manner using a realistic 3D environment that limits incurring safety issues and danger of damage to resources. To obtain superior performance in an Industry 4.0 setup, a modiﬁed version of a binary gravitational search algorithm is introduced which beneﬁts from an exclusive or (XOR) operator and a repository to improve the exploration property of the algorithm. Mathematical analysis of the proposed optimization approach is performed which resulted in two theorems which show that the proposed modiﬁcation to the velocity vector can direct particles to the best particles. The use of repository in this algorithm provides a guideline to direct the particles to the best solutions more rapidly. The proposed algorithm is evaluated on some benchmark optimization problems covering a diverse range of functions including unimodal and multimodal as well as those which suffer from multiple local minima. The proposed algorithm is compared against several existing binary optimization algorithms including existing versions of a binary gravitational search algorithm, improved binary optimization, binary particle swarm optimization, binary grey wolf optimization and binary dragonﬂy optimization. To show that the proposed approach is an effective method to deal with real world binary optimization problems raised in an Industry 4.0 environment, it is then applied to optimize the assembly task of an industrial robot assembling an industrial calculator. The optimal movements obtained are then implemented on a real robot. Furthermore, the digital twin of a universal robot is developed, and its path planning is done in the presence of obstacles using the proposed optimization algorithm. The obtained path is then inspected by human expert and validated. It is shown that the proposed approach can effectively solve such optimization problems which arises in Industry 4.0 environment.


Introduction
It is highly desirable to prepare products to be served and used for customers with their own dominant needs, interests, energy and use standards. Industry 4.0, the forth generation of industry, is a response to such varied requirements which is predicted to lead to highly customized and re-configurable production. Such technology rely on big data and benefits from machine learning and artificial intelligence approaches. The use of such approaches may even make it possible to update production procedure to give appropriate response to fast changes in the design of product or make product customization possible [1]. To implement this on the factory floor, it is required to use robust intelligent techniques including advanced machine learning and artificial intelligence approaches to make the production line capable of implementing such varied production modifications. To avoid safety issues as well as wasting resources, a digital twin [2], a precise electronic replica of the factory floor, can be utilized. Such digital environments make it possible to design a robust system at lower cost and times to implement.
There exist different Industry 4.0 applications, in which the decision variables are binary rather than real values. Optimal task planning [3], flow shop scheduling [4] and assembly line design problems [5,6] are example applications of binary optimization algorithms in a digital manufacturing setup. The competition between different factories emerges mass production, while maintaining robustness and giving appropriate response to continuous design changes. In such an environment, it is required for the production elements, including robots and other automated machines, to perform task planning automatically to minimize the costs and time consumed, while maximizing production rates and profit. More specifically, rapid changes in products shortens the life-cycle time of production systems [7]. For certain product assembly lines, long-term shut down may be required to redesign the line and update recipes. Automated optimization of such processes are highly commercially beneficial to reduce these inefficiencies and improve production. Intelligent optimization algorithms may be beneficial for automatizing process design on the factory floor [2] which can consider energy, time and other processes. The requirement for optimization algorithms motivate us to develop a new algorithm, study its mathematical analysis and use it to perform an automated optimization process in Industry 4.0 applications.
Although classical optimization algorithms such as linear programming approaches, calculus of variations and gradient based optimization methods have proven to be efficient optimization algorithms, there are some barriers against use of these algorithms in all applications. When the function to be optimized is not explicitly defined and/or its partial derivative with respect to its parameters is hard to obtain because of high dimension and discontinuity, classical optimization processes fail. Intelligent optimization algorithms may be good candidates for high dimensional optimization problems. Intelligent optimization algorithms benefit from multiple start point, they are derivative free and are capable of finding multiple global optimums. Moreover, since these optimization algorithms are initiated from multiple points, they are less sensitive to initial conditions. Evolutionary optimization algorithms [8], swarm intelligence [9] and physics-inspired optimization methods [10] are the main categories of intelligent optimization methods applied in these cases.
Each solution in swarm intelligence approaches is defined in terms of a d-dimensional vector representing position. A velocity vector updates the positions for the next generation. Velocity vector benefits from exploration as well as exploitation terms. The third category of optimization algorithms is physics-inspired approaches in which, other than position and velocity, each particle benefits from more physical properties such as force, mass, electrical charge, etc. Examples of physics-inspired optimization algorithms are Big Bang-Big Crunch [11], central force optimization algorithm [12] and gravitational search algorithm (GSA) [13].
Among physics-inspired optimization algorithms, GSA is one of the most successful optimization algorithms [14]. Each particle in this algorithm benefits from mass, position vector, velocity vector and acceleration vector. The value of mass assigned to each particle is proportional to its merit with highest value being assigned to the best particle and the lowest value being assigned to the worst one. The gravitational force acting on these particles make the worst particle accelerated towards a low number of "best" particles, while searching the whole space for a possible "better" solution.
To apply real-valued optimization algorithms to binary optimization problems, they are required to be modified considerably. In swarm-based and physic-inspired optimization algorithms, a randomly weighted velocity vector is added to position vector to obtain new position. The well-known modification to such problems is to change the definition of velocity to probability of change in bits of each dimension [15,16]. However, the analysis in this paper shows that in the case of binary gravitational search algorithm (BGSA), the probability vector does not always behave as expected.
To alleviate this problem, the velocity vector and consequently the probability of change is modified in this paper to include an exclusive or (XOR) operator.
In this paper, the way mass values are assigned to each object, in the current BGSA, acceleration equation and formula for update of the change probability and position are fully analysed. The acceleration term is modified to use an XOR term. A repository is added as well to improve the performance of the optimization algorithm. To test the efficacy and performance of the proposed optimization algorithm, it is validated in simulation against a number of benchmark optimization problems. It is shown that the proposed algorithm is capable of optimizing such benchmark optimization problems with superior performance when it is compared to BGSA [16], improved binary particle swarm optimization (IBPSO) [17], binary particle swarm optimization (BPSO) [15], binary grey wolf optimization (BGWO) [18] and binary dragonfly optimization algorithm (BDA) [19].
Next, two Industry 4.0 applications namely assembly task planning and industrial robot path planning are considered. Assembly task planning is then validated on a digital twin first and then implemented on a physical robotic system. For the assembly task planning, a mock industrial calculator is designed and 3D printed. The studied task is to put bolts in their place using a robot manufactured by universal robot company namely UR5. Such an optimization process makes the production line robust as any new configuration for nuts and bolts can be easily given to the robot and it will optimally solve the sequence. It is shown that the proposed optimization algorithm is capable of planning this task and the length of path travelled by the robot is decreased by almost 17% in comparison to the human generated original pathway which leads to reduction in time and cost. Furthermore, path planing of the same robot in the presence of obstacles in its digital twin is considered. It is shown that the proposed approach is successful in finding an optimal path. In real-world practice the path generated in the digital twin would then be sent for human expert approval before transfer to physical robot. Such an automated optimization task reduces the need for intervention, yet keeps human elements involved as a safety measure. This paper is organized as follows: In Section 2, real-valued GSA is discussed. The current version of BGSA investigated in this paper is presented and analysed in Section 2. In Section 3, the proposed novel XOR BGSA is introduced and analysed. Moreover, a repository is added in Section 3 to improve the performance of XOR BGSA. Full theoretical analysis of the proposed method is given in Section 4. The proposed algorithm is then compared with several existing optimization algorithms for the optimization of several benchmark examples in Section 5. Implementation of the proposed algorithm on UR5 is presented in Section 5 as well. Finally, the concluding marks are presented in Section 6.

Gravitational Search Algorithm
GSA is an optimization algorithm inspired by the physical forces between objects due to Newtonian forces between masses. Compared to PSO, in this optimization algorithm, each object benefits from more features such as acceleration term and the value of mass which is assigned with respect to the merit of each particle [13].
In this optimization algorithm, better particles are assigned a higher value of mass which causes other objects to be accelerated towards these better objects while still scanning the whole search space. In this case, if a particle come up with a better solution, it replaces the position of previous good solution found and attracts other objects. After a few iterations, all solutions are converged towards the heaviest mass resulting in termination of the algorithm. The algorithm is explained mathematically in the upcoming subsection [13].

Basic Gravitational Search Algorithm
Each solution in d-dimensional search space is represented by the position of an object and is represented as follows: where x j i is the jth component of the position of object and N is the total number of objects. The mass of particles are updated and normalized at tth iteration as follows [13].
where m i (t) is non-normalized mass value corresponding to ith particle at iteration number t which represents the quality of solution and is defined as follows [13].
where f worst (t) is the fitness function value corresponding to the worst particle and f best (t) represents the best fitness function observed by any of the particles. Hence, in this case m i (t) satisfies m i (t)∈[0, 1] and the value of mass of best solution is equal to one while the value of mass corresponding to worst solution is equal to zero. The parameters f worst (t) and f best (t) are updated at every iteration as follows: It is to be noted that in this case the best and worst particles are defined based on minimization problem. In maximization problems, the updating equations for f worst and f best swap their definition. Furthermore, k b number of best solutions are selected at each iteration as the masses which attract other masses.
The overall gravitational force acting on the ith particle is obtained as follows: where k b represents the number of best solution selected, . stands for Euclidean norm, ε is a small value added to prevent division by zero, r p is the power considered for the Euclidean distance between two particles, G(t) is the gravitational constant and r j ∈[0, 1] is a uniform random value. The gravitational constant is updated at each iteration using the following equation.
where G 0 has a constant real value and t max is the maximum value of the iterations of the algorithm. According to Newton's second law of motion, the acceleration of particles are calculated by dividing the applied force to ith mass by its value of mass as follows: where A i (t)∈R d is d-dimensional acceleration of the particles. It is further possible to update the velocity of objects using the following equation.
where V i (t)∈R d is d-dimensional acceleration of the particles at tth iteration and p i ∈[0, 1] is a uniform random number. Finally, the newer position of each particle is updated as follows:

Binary Gravitational Search Algorithm
The modification made to GSA to solve binary optimization problems is similar to BPSO [9]. In this case, the velocity vector in BGSA changes its meaning to the probability of change. Hence, a large absolute value for velocity in a dimension results in high possibility of change in such direction. Based on such a description, the function P(V j i ), the probability of change in jth bit of i th object, is considered to be a hyperbolic tangent of the velocity vector as follows [16].
where tanh(.) represents the hyperbolic tangent function as follows: Calculating the P(V j i ), the probability of change in every dimension, the next state in each dimension is calculated as follows [16].
where x j i (t) represents the complement of x j i (t) and r j i ∈[0, 1] has a uniform random value.

Analysis of Binary Gravitational Search Algorithm
An analysis of BGSA shows that in certain cases, the algorithm does not behave as expected. Focusing on acceleration term of BGSA as in (7), the following analysis is performed in different cases.
As the value assigned to parameter X k (t), k= 1, ...,N in each dimension may only be any binary value of zero and one, the result of subtraction X j (t)−X i (t) in each dimension may be any of following possible conditions. In the following cases, undesirable behaviour of velocity vector, where occurred, is shown under unchanged sign of the velocity vector condition.
(1) In the case when X k One has X k i (t)−X k j (t) = 0 which results in no acceleration for the dimension k which is reasonable. In this case, the value of bits in the kth dimension for the corresponding object and the best solutions are exactly the same and no change is required in that dimension.
(2) If X k j (t) = 1, X k i (t) = 0, X k j (t)−X k i (t) = 1 resulting in a positive acceleration term. In this case, if the velocity in kth dimension is positive, this acceleration term works fine and adds up to the speed. However, if the velocity in kth dimension is negative, the acceleration term decreases the absolute value of velocity and makes the probability of change smaller which is completely undesirable.
(3) If X k j (t) = 0, X k i (t) = 1, X k j (t)−X k i (t) = −1 results in a negative acceleration term. In this case, if the velocity in kth dimension is positive, this acceleration term decreases the velocity and consequently the probability of change. Since X k j (t) is different from X k i (t), this behaviour is not desirable. On the other hand, if the velocity in the kth dimension is negative, the acceleration term increases the absolute value of velocity and makes the probability of change larger, which is desirable.
Motivated by the cases in which the acceleration term for BGSA does not behave as desired, this algorithm is modified using an XOR operator. It is shown that using such acceleration term, the undesirable conditions will never be met again.

Proposed XOR Binary Gravitational Search Algorithm
To alleviate the problems mentioned in Section 2, the BGSA is modified to have an XOR operator in its update rule for the acceleration term as follows: where X rj (t) is the jth solution in the repository and ⊕ represents the XOR operator acting bit-wise on X rj (t) and X i (t) with its truth table being presented in Table 1. Furthermore, the equation used for G(t) is as follows.
where α 0 and α 1 are selected as 0.8 and 0.2, respectively. The probability function P(V j i ) is modified as follows: where tanh(.) is defined as in (11). The update equation for the velocity vector as well as position vector remains the same as (8) and (9). Moreover, the mass parameter in the proposed algorithm associated with particles in repository is updated as follows: and m ri (t) is non-normalized mass value corresponding to ith particle at iteration number t which represents the quality of solution and is defined using the fitness of particles in repository as follows: where f (X ri ) is the fitness function associated with particle X ri in the repository, f r,worst (t) and f r,best (t) represent the worst and the best fitness function of particles in repository. While the original version of BGSA is memory-less, in the proposed algorithm a repository is added to the algorithm. In each iteration, a tournaments selection procedure is performed between the current generation and old repository and the new repository is updated as a result of this selection procedure. Tournaments selection is widely known to be an algorithm to preserve the diversity in the solutions. The flowchart of the proposed XOR BGSA with repository is depicted in Figure 1, where V max refers to the maximum absolute value of velocity in each dimension.

General Analysis of the Proposed XOR BGSA
The proposed XOR BGSA is investigated statistically to show that it is more compliant with the expectations of acceleration and velocity vectors of BGSA. Since the parameter X k (t), k = 1, . . . , N can only accept binary values of zero and one, the result of XOR operation X j (t)⊕X i (t) may fall in each of the following categories.
It is expected from the acceleration term that since the corresponding bit in the best particle is the same as the bit in the current particle, the probability of change reduces. This expectation is fulfilled using (13).
the acceleration term introduces a positive value which contributes positively to the probability of change and increases its value.
(3) In the case when X k Similar to the previous case, since the corresponding bit in the better particle is different than that of the current particle, a positive acceleration term is obtained which results in higher probability of change. This is an expected behaviour from the system.
According to this analysis, the issues mentioned earlier for BGSA do not exist in the proposed XOR BGSA.

Full Statistical Analysis of the Proposed Method
To perform theoretical analysis of the proposed XOR BGSA, a single dimensional optimization problem is considered. The parameter k b in this case is equal to 1. Hence, in this section, we have a single particle whose behavior is studied and one best particle in the repository. The acceleration term, velocity term and the probability of change in the case of single dimensional optimization problem is obtained as follows: where a(t), v(t), P(v) are single dimensional vectors representing the acceleration, velocity and the probability of change corresponding to the studied particle. Furthermore, x (t) is the single dimensional vector representing the current position of studied particle and x rb is its corresponding value in the best particle in the whole population and M rb corresponds to the mass of the best particle in the repository. The minimum value of G (t) is equal to α 1 . Two cases are investigated: the case when p i = 1 and the case when 0 < p i <1 . Case I: Proof. Considering (13), we have the following equation for the acceleration term a(t): where D = 1 ε is the dimension of the optimization algorithm. Consequently, considering (18), one obtains the following inequality: This implies that v t is monotonically non-increasing function of time.
which means that ∀t 0 < T Lemma 2 holds. However, in the case when v≤v(t 0 − 1), we have: Since x(t) =x rb , we have the following equation: Considering the fact that minG (t) =α 1 , one obtains the following inequality: To find the appropriate T, the following inequality needs to hold.
Considering the fact that r (t) is a random variable uniformly distributed in [0, 1] , the following equation holds.
By applying expected value of r (t) as in (27) to (26), the following equation is obtained.
and finally T is obtained as follows: Since <v. This concludes the proof. Remark 1. Lemma 2 implies that in the case when the position in one dimensional optimization problem is the same as that of best particle, the velocity vector corresponding to that dimension grows in negative direction resulting in the probability of change to vanish to zero. This in turn results in no change in that dimension.
Considering (17), we have the following equation for the acceleration term a(t): Consequently considering (30), one obtains the following inequality: This implies that v t is monotonically non-decreasing function of time.
Proof. In the case when v < v(t 0 − 1) since according to Lemma 3, v(t) is monotonically non-decreasing function of time we have: which agrees with the Lemma 4. However, in the case when v(t 0 − 1)≤v, we have: Considering the fact that x(t) =x rb , the following equation is obtained.
Considering the fact that minG (t) =α 1 , the following equation is obtained.
Hence, to find T which satisfies Lemma 4, the following inequality needs to be hold.
Considering the fact that r (t) is a random variable with uniform probability distribution function, the following equation holds.
By using the property of r (t) as in (36), the following is obtained.
Hence, the value of T to satisfy Lemma 3 is obtained as follows: This concludes the proof.

Remark 2.
It immediately follows from Lemma 4 that in the case when the position in one dimensional optimization problem is different than that of the best particle, the velocity vector corresponding to that dimension grows resulting in the probability of change to increase to one which would result in increase in the probability of change.
1. In the case when x(t)=x rb , the acceleration term is non-positive and the velocity term is a non-increasing function of time and the probability of change converges to zero. 2. In the case when x(t) =x rb , the acceleration term is non-negative and the velocity vector is a non-decreasing function of time and the probability of change converges to zero.
Case II: 0 < p i < 1: In the case when x(t) =x rb , ∀t > 0 the velocity vector v(t) will never be larger than p t i v(0), where v(0) is the initial value of the velocity.
In the case when x(t) =x rb , from (17), it follows that Remark 3. From Lemma 5, it is concluded that in the case when x(t) = x rb , the acceleration term is non-positive acting in the direction of decreasing the probability of change. However, the impact of p t i v(0) depends on the sign of v(0).

Lemma 6.
In the case when x(t) = x rb , ∀t > 0 the velocity vector v(t) will never be smaller than p t i v(0). In the case when x(t) = x rb , from (17), it follows that Considering (18), the following equation is obtained for v(t) Remark 4. From Lemma 6, it is concluded that in the case when x(t) = x b , the acceleration term is non-negative acting in favor of making probability of change in the corresponding bit higher. However, the impact of p t i v(0) depends on the sign of v(0).

Theorem 2.
In XOR BGSA, considering sufficiently large number for maximum number of iterations, following properties hold.
1. If x(t) = x rb , ∃T s.t. ∀t > T, the probability of change becomes smaller than 0.5: ∀t > T, the probability of change becomes larger than 0.5: In the case when p i = 1, from Theorem 1 it follows that when x = x rb , as mentioned earlier in Remark 1, the velocity value grows in negative direction which results in: Moreover, in the case when p i = 1, from Theorem 1 it follows that when x(t) = x rb , as mentioned earlier in Remark 2, the velocity value grows larger which results in: Hence in the case of p i = 1, Theorem 2 holds. When 0 < p i < 1, two cases are considered. The case when v(0) = 0 and v(0) = 0.
If v 0 = 0, in the case of x = x rb , the velocity term according to (18) is obtained as follows: Hence: In the case of x(t) = x rb , the velocity term according to (13) is obtained as follows: Hence: If v(0) < 0, from (49), it can be seen that v (t) is the result of addition of two negative terms resulting in a negative value for v (t). However, if v(0) > 0, considering the fact E [r(t)] = 0.5, we have: Hence, the time T after which we have E[v (t)] < 0, is obtained from Hence, after a sufficiently large value of iterations specified by T, v (t) becomes negative and stays negative onwards. After that considering (19), we have the following: If v(0) > 0, from (54), it can be seen that v (t) is the result of addition of two positive terms resulting in a positive value for v (t). However, if v(0) < 0, similar to the analysis which resulted in (50) and considering the fact E [r(t)] = 0.5, we have: Hence, the time T after which we have v (t) > 0, is obtained from the following inequality Hence after sufficiently large value of iterations specified by T in (57), v (t) becomes positive and stays positive on-wards. After that, considering the probability of change as in (19), it follows that: This concludes the proof.
Remark 5. It immediately follows from Theorem 2 that if the value of a particle in one dimension is different than that of the best particle, the probability of change if not immediately, after sufficiently large number of iterations becomes larger than 0.5 making change in the corresponding dimension more probable than not. Moreover, if the value of particle in a dimension is similar to that of the best particle, the probability of change in that dimension if not immediately, after sufficiently large number of iterations becomes smaller than 0.5 making change in the corresponding dimension less probable.

Remark 6.
The results obtained in Section 4 are valid for the case when k b = 1. However, these results can be easily extended to the case when k b has a greater value in which case the value of k b particle effect other particles with a gain proportional to their mass. Hence, the overall probability of change is influenced by each of k b particles disproportionately.

Benchmark Optimization Problem
To evaluate the efficacy and performance of the proposed XOR BGSA with repository, it is tested against several benchmark problems. In optimization of benchmark problems, 20 bits are used to decode numerical values into the binary domain. The number of solutions considered (population size) is equal to 25, maximum number of generations are then selected to be equal to 1000. Two sets of comparisons are presented. In the first set of optimization, the proposed XOR BGSA is compared to BGSA [16], IBPSO [17], BPSO [15], BGWO [18] and BDA [19]. The main goal in these simulations is to demonstrate how the modifications proposed in this paper including the addition of XOR to the acceleration equation for BGSA as well as the repository effect its performance.
As in the real world, an optimization problem may be applied to a wide range of different functions including unimodal and multimodal problems, the proposed optimization algorithm is tested against both classes of optimization problems. The ability of the proposed optimization algorithm in dealing with problems with multiple local minimas is required to be investigated as well. For having such study, some of the optimization problems considered as the benchmark optimization problems include multiple local minimas. The benchmark problems considered are a collection of functions with single optimum solutions and functions with multiple local optimums. Other than these specifications, it is required to illustrate the performance of the proposed optimization algorithm in dealing with non-differentiable problems. The specifications of benchmark optimization problems are mentioned in front of them for ease of reference. The parameter N represents the dimension of these optimization problems which is considered to be equal to either 1, 2, 5 and 10 to test the optimization algorithms for different dimensions. These dimensions are multiplied by the number of bits taken for decoding real values, which results in optimization in 20, 40, 100 and 200 dimensional space. The mathematical formulas for the test functions as well as their specifications are as follows [20].
Sphere function [21]: a differentiable, separable, unimodal function without local minima Multiplication of squares: a differentiable, non-separable, unimodal function without local minima Griewark function [22]: a differentiable, non-separable, multimodal function with local minimas Schwefel function [23]: a differentiable, separable, multimodal function with local minima Rastrigin Function [22]: a differentiable, separable, multimodal function with local minima Styblinski-Tang [24]: a differentiable, separable, multimodal function with local minima Michalewicz function [22]: a differentiable, separable, multimodal function with local minima 20 (65) Ackley function [25]: a differentiable, non-separable, multimodal function with local minima Rotated hyperellipsoid function [22]: a differentiable, scalable, unimodal function without local minima Alpine 1 function [22]: a non-differentiable, separable, multimodal function with local minima Cosine mixture function [22]: a differentiable, separable, multimodal function with local minima The results of comparison between the proposed algorithm and five other existing binary optimization algorithms are presented in Table 2. Optimization process in each case is repeated for 30 times and the mean values are reported. It can be seen that the proposed optimization algorithm performs well in optimizing uni-modal as well as multi modal optimization problems. In 28 cases out of 44 cases, the proposed method outperforms other investigated optimization algorithms. It is further observed that when the proposed method is not the best optimization algorithm, it is still very close to the desired results. For instance, in the case of optimization function number 1, a uni-modal optimization problem, the proposed version of BGSA outperforms other algorithms in all different number of dimensions except for the case of n = 1 where it is the second best optimization algorithm after BDA. Optimization function number 11, f 11 (x), is a multi-modal optimization problem with several local minimums for which the proposed XOR BGSA outperforms the other five investigated optimization algorithms in all four dimensions. It can further be observed from the Table 2 that BDA ranks the second after the proposed optimization algorithm. IBPSO [17] BPSO [15] BGWO [18] BDA [19] f 1 (X)

Application to Knapsack
As the second class of benchmark optimization problem, the multi knapsack problem is considered. The problem formulation of 0/1 knapsack problem is discussed in summary in this part. Let there be n items with each having a positive profit p i > 0 and a positive resource consumption r i > 0. The main objective is to fill several knapsacks with resource capacity of C j such that the accumulated profit of the items is maximized. Considering the constraints imposed by capacity of knapsacks, this problem is a constrained binary optimization problem. Such problems arises in various real-world applications such as in a distributed computer systems, warehouse optimization, cargo loading and stock problems [5,26]. This problem is explicitly formulated as follows: x i ∈ {0, 1} ∀i∈{1, . . . ,n} where n represents the item number, m is the number of constrains, p i ≥ 0 represents the profit value for each item, resource consumption is expressed by r ij ≥ 0 and the capacity is denoted by another positive value C j ≥ 0. The cost function is multiplied by negative value of −1 to change maximization to minimization problem. The constraints are added as penalty terms to the cost function to obtain an unconstrained optimization problem as follows: where the parameter β is used as the penalty coefficient and assigned the value equal to 10 10 . Several standard databases for knapsack optimization problems exist to provide means for comparison between the algorithms which provide special values for p i 's, C j 's and r ij 's. Among these database of multi knapsack optimization available at the website of Queen Mary University of London (http://www.eecs.qmul.ac.uk/~jdrake/mkp.html) formerly hosted at University of Nottingham (http://www.cs.nott.ac.uk/~jqd/mkp/index.html) are selected for performance comparison of the proposed algorithms with several state-of-the-art binary optimization algorithms. The total number of populations are taken to be equal to 10 with maximum number of iterations selected to be equal to 1000. Table 3 presents the results obtained using the proposed optimization algorithm versus other optimization problems. The solutions are the optimization solutions to the function represented by (72) which are multiplied by minus one to reflect the solutions to the original problem of (70) and (71). There exist several negative values in the table which show that in some cases, algorithms could not fulfil the constraints of optimization at least once. Such problem does not occur for the proposed XOR BGSA with repository. The best result obtained for each dimension is made bold faced to demonstrate its superior performance. The values used for optimization are collected from the website of Queen Mary University of London and their file names are presented for ease of reference. Problems have also been given a number for further discussions. As can be seen from Table 3. The proposed XOR BGSA with repository outperforms all other investigated optimization problems in almost all cases.
The convergence property of the proposed approach versus previous optimization methods for optimizations number 7 to 12 are demonstrated in Figure 2. As can be seen from this figure, the convergence of the proposed algorithm towards the solution is well above other investigated algorithms. Although during the first few iterations, some of the algorithms rise faster towards optimal solutions, after 100 iterations, the proposed XOR BGSA speeds up and becomes the dominant one. It can be concluded from the figure that other than quality of solutions, the speed of the proposed approach is well above other investigated optimization methods.  IBPSO [17] BPSO [15] BGWO [18] BDA [19] 13 500 0. 25  In most real-world optimization problems, the cost function evaluation requires high volume of computation. Hence, it would be highly desired for an optimization algorithm to result in superior performance after its single run. Box plot is a graphical method which gives more information about the diversity of the results rather than just its mean value. A small standard deviation leads to a smaller interquartile range which can be interpreted as its ability to avoid local minima. Outliers can also be easily spotted from box plot. In Figure 3 while the interquartile is illustrated using a blue box, the minimum and the maximum associated with data are demonstrated using the black line. There is a 50 percent chance that for a new experience, the result falls within the interquartile range. The outliers in this figure are represented by "+" sign. The median values are demonstrated with red line. The box plot obtained using different optimization algorithms for optimization problems 7 to 18 are illustrated. It can be observed from the figure that not only the median values obtained for the proposed approach in 10 cases out of 12 cases are the best, interquartile range obtained for the proposed approach in most cases is smaller meaning that it is more repeatable than in other investigated optimization algorithms. Furthermore, for optimization problems numbers 7, 8, 11 and 12, it can be observed that the interquartile range obtained for the proposed approach is well above the other optimization algorithms. This means that for every single run of the proposed optimization algorithm, there is an appreciable chance we can obtain superior performance for the proposed optimization algorithm compared to other investigated optimization methods. Figure 2. Trend of Cost function obtained from proposed XOR BGSA vs. four other algorithms namely: BGSA [16], improved binary particle swarm optimization (IBPSO) [17], binary particle swarm optimization (BPSO) [15], binary grey wolf optimization (BGWO) [18] and binary dragonfly optimization algorithm (BDA) [19]. Knapsack optimization problems numbers 7 to 12.  Figure 3. Box plot obtained for the knapsack problem using the proposed XOR BGSA vs. three other algorithms namely: BGSA [16], BGWO [18] and BDA [19]. Knapsack optimization problems 7 to 18.

Industry 4.0 Applications
In current manufacturing environments, modifications made to a product often require redesign and retraining/programming of production elements on the factory floor. Optimal or sub-optimal solutions to design changes in a process is the goal of decision-making processes. Recent work has explored showing that systematic optimization is possible in a well-defined detailed digital twin in which searching algorithms can perform multiple trial and errors without safety and resource restrictions. Additionally, current robot programming techniques require highly talented professionals to perform programming, which increases response times and increases costs [27]. Real-time decision making and robot programming is therefore highly desirable and is a stated requirement of Industry 4.0 [28]. However, autonomous programming of robots will impact safety and certification requirements as well. To mediate such impacts, one solution is to tie the digital twin solutions to the real-production methods. Such a digital twin needs to be as accurate as possible to ensure compatibility of actions with real world. The simulated movements in such an environment can then be transformed to the real robot after full examination of resulting optimization solutions by a human expert. This framework makes the most of machine learning and artificial intelligence approaches. In this part, the proposed XOR BGSA with repository is used to find an optimal solution for assembly task planning and motion planning for an industrial robot namely UR5 that would then be passed to a human for certification before implementation. Previously the proposed XOR BGSA is tested against different types of optimization problems including uni-modal, multi-modal, those with single optimum solution, those which suffer from multiple local minimums and knapsack optimization problem which is a constrained one. The fact that the proposed optimization procedure obtained the best results among six binary optimization problems in most cases motivates us to use it in Industry 4.0 optimization problems.
It is required for robots to be able to optimally plan their own task autonomously otherwise massive resources may be used to redesign, reconfigure and reprogram them. To fulfil this goal, it is required to decode the planning process in the form of an optimization problem and use artificial intelligence approaches to solve the problem on a UR5 capable of handling 5 Kg loads to a precision of order of 0.1 mm .
The assembly task planning to be optimized in this application is to put bolts in their respective places while minimizing the distance. The bolts used in this study are from three different types of of M6 × 20 mm, M8 × 20 mm and M10 × 20 mm. Hence, we have three different type of nuts which must be fastened using the bolts. The place and the position of the nuts and bolts are depicted in Figure 4 with the z component of position for nuts being on 5.7 cm and for bolts being equal to 7 cm.
Such an assembly task is similar to the benchmark problem frequently solved in computer science projects called the travelling salesperson. The cost function to be optimized in this case is the total distance travelled with main difference between the assembly task designed to test the proposed XOR BGSA and the travel salesperson being that in the current experiment, we have certain bolts which can be used to fasten certain nuts; this imposes some constrains on the optimization problem.
To encode the optimization problem, 31 bits are considered for each individual. Inspired by the results obtained by the proposed optimization procedure, it is used to minimize the total distance travelled by UR5 during the assembly task of industrial calculator (See Figure 5). In this case, there exist 2 31 = 2.1475 × 10 9 different combinations for performing the assembly task. As can be seen from Figure 5, the first two bits determine which bolts are used to fasten nut #1, the bits numbers 3 and 4 are used to determine which bolts are used for nuts number 2, the fifth bit is for nut number 3. Since we have four M10 × 20 mm bolts, the remaining bolt goes to nut number 4. Bit number 6 determines which M8 × 20 mm is used to fasten nut number 5, the other one is used to fasten bolt number 6. The bit number 7 is used for bolts numbers 7 and 8. However, since the order of performing these tasks effects the cost function, it is required to determine the order of performing the tasks using the next 24 bits. All 3 of them making an integer value from the interval of [0,7]. In this case, if the ith value after converting to decimal value is equal to j, the ith and jth task swap their position. This results in avoiding forbidden conditions which is difficult to avoid otherwise. The cost function to be optimized is the total distance travelled by robot with its initial conditions being at point (10,30,10).
The evolution of total distance travelled by the robot is depicted in Figure 6. The total number of particles in this optimization is chosen to be 10 with the maximum number of iterations being equal to 200. As can be seen from this figure, the best random initialization for the total distance travelled by robot is 328.0 cm which decreases to 273.2 cm after optimization, resulting in a 16.7% improvement. The optimum solution of the problem is represented in Table 4 and the assembly process done by UR5 is illustrated in Figure 7. The total assembly time is 116 seconds. Where for our laboratory safety reasons, the joint speed of the robot is limited to 90/s. However, as the UR5 is capable of working with joint speed equal to 180/s, the total assembly time could be reduced if it works under full speed. Table 4. Optimal assembly solution obtained from the optimization algorithm.

Motion Planning of Industrial Robots Inside a Digital Twin Created by V-REP Software
Proposed XOR BGSA with repository is used to find the optimal path for an industrial robot in its digital twin environment. The industrial robot used in this part is of a cobot type namely UR5. Cobot is a slang which refers to collaborative robot, generally accepted as any robot designed to operate in close proximity to human workers without need for a separating cage. The problem considered here is to find the shortest collision free path for the robot in 3D space while undergoing a pick and place task. Figure 7 illustrates the digital twin of the manufacturing environment simulated in V-REP software, where there exist a human worker working closely to the robot.
Simulations have been carried out using V-REP software, a user-friendly software, which works under various operating systems including Windows. This software V-REP is designed around a versatile architecture with different independently working functions. The options embedded in this software can be easily enabled or disabled as required. The communication with this software is possible from ROS using rostopics as well as directly from Python or Matlab using respective API developed specially for this software. A wide range of robots including industrial robots from different companies including Universal Robots are available in this software. The driver included in this software for each of these robots possesses much of real robot features including stopping movement in the case of collision with an object. Real features of robot including maximum speed, maximum joint angles and maximum torque are available in the robot and may be modified at user desire making it replicate the real robot as much as possible. Another feature which exist under V-REP software is its ability to compute the inverse kinematic of robot and find the joint angles of the robot automatically which considerably facilitates robotic simulation.
The obstacles in this environment are possible obstacles in real manufacturing environments such as other machines, tools and devices. V-REP software is a robust software for simulating the physical behaviours of objects, including collision between them as well as forces such as gravity and friction that can all be included as part of the software environment. In this way, the user can easily define object positions and orientations while the V-REP software takes care of updating movement of the various other objects and collisions which may happen during motion. As the result is a game like 3D interface, examination of the movement can be easily achieved by the human expert. In this way the human expert can ascertain if there may be problems with the solution that are not specified in the programming. To perform the automated optimization, it is required that the robot movements be encoded in terms of an evolvable array which may then be translated to motions to evaluate cost function and make changes to optimize the problem. A similar approach has been previously done using Ant colony optimization in another software namely Gazebo [28]. However, the benefit of using V-REP versus Gazebo is that the former benefits from built-in inverse kinematic which makes it possible to perform the simulations in a timely manner. A comparative study has been performed previously in [29] which compared V-REP software and Gazebo. It is concluded that V-REP is a more user-friendly and intuitive software than Gazebo in which intelligent optimization can be carried out with more comfort.
While Figure 8 depicts the real robot, its digital twin developed in V-REP is depicted in Figure 9. As can be seen from these figures, digital twin matches with real robot in terms of dimensions and position of robot and obstacles, respectively. As can be seen from this figure, several objects as well as the UR5 robot exist in this digital twin. The obstacles are in the form of several cubic boxes which replicate the machines and devices which exist close to the robot.

Encoding Movements
To optimize movements, they are encoded in terms of an evolvable array according to the proposed XOR BGSA with repository which is implemented under Matlab software. In this paper, it is assumed that the initial orientation of the robot remains the same during optimization and its position needs to be changed to perform the displacement from the start point to the end point. Considering three x, y and z dimensions and the fact that they can be increased and decreased, there would be six movements which can be represented by any integer value from the set S = {0, 1, . . . , 5}. However, the proposed optimization algorithm would generate three bits for each movement. To avoid infeasible solutions, the set is considered to have all eight conditions and duplicate movements may be added. Moreover, during optimization duplicate coordinates may occur as a result of back and forth in the movements. To avoid such loops in robot movement, all coordinates are compared with each other and all intermediate movements which result in a loop are eliminated. The interpretation of movements are presented in Table 5. Table 5. Interpretation of values of array generated by the proposed XOR BGSA.

Kinematic and Inverse Kinematic of UR5
The UR5 robot benefits from 6-degrees of freedom. Position and orientation of this robot and its end-effector in 3D work-space is represented by x ∈ R 6 . While the first three dimensions of this vector present three translational variables, the next three dimensions represent its orientation. The joint angles, represented by q = q 1 q 2 . . . q 6 , are solutions to the inverse kinematic problem.
The forward kinematic function is the function used to compute the position and orientation of the robot corresponding to its assigned joint angles.
The standard Denavit-Hartenberg (DH) convention to represent the forward kinematics of the robot is used in this paper. The DH parameters associated with UR5 are presented in Table 6. Using such convention, inverse kinematic of the robot can be found as follows: where T i−1 i presents the homogeneous transformation matrix which makes the transform from frame i − 1 to frame i defined as follows: where s(.) and c(.) represent the sinus and cosine of their arguments. The values obtained from the movement array need to be translated in terms of robot movements. To enable this a transformation matrix that represents position and orientation of robot is obtained based on work by [30,31]. The transformation matrix for frame 6 with respect to frame 0 is as follows: where 0 6X , 0 6Ŷ and 0 6Ẑ are unit vectors defining the axis of frame 6 in relation to frame 0 and: in which the rotational matrices are obtained as follows: Inverse kinematics are then used to find joint angles of robot to move to a desired position given by proposed XOR BGSA. The angles of the robot are then obtained as [30,31]: where:

. Communication between Matlab and V-REP
Communication between Matlab and V-REP is done using API developed by Coppelia robotics. Figure 10 demonstrates the overall communication requirement between Matlab and V-REP. As can be seen from this figure, proposed XOR BGSA with repository generates the path in terms of a binary sequence. This binary sequence is then translated to real movement with required resolution. Such movements are transmitted to V-REP in terms of positions and orientations of robot which is translated to joint angles of robot using built in inverse kinematic software. When a collision occurs, the target position of robot needs to be modified to the last possible point before collision in the system. After performing simulation, the results including total simulation time and distances are reported back as the real-world perception of the binary code in terms of a cost function. After all particles are evaluated, the proposed XOR BGSA iterates for another iteration to generate the next solutions.
The parameters x i , y i , z i , a i , b i and c i are the position and orientation of robot after every movement. Since obstacles are static objects in V-REP, in case of a collision, the robot stops, and a collision flag is returned. Furthermore, the desired position given to the robot is not met because of collision. Hence, the difference between desired and obtained position of robot is a measure of collision. Such position difference can be used to penalize the cost function. The increment in the position suggested by the proposed XOR BGSA would be the actual position of robot plus the decoded increment from Table 5.

Cost Function to Be Optimized
As the goal of the optimization is to obtain a desired position and orientation while ensuring that no collision happens, total distance between the desired position and the actual position of the robot is as follows: where x d , y d and z d are coordinates of P D the desired position of the robot. Other than distance to the desired position, total time is used as another term in the cost function.
where T represents the total time, D C represents the distance between the actual position of robot and its desired position in every single movement due to collision with obstacles. The parameters K D and K T are two gains used for optimization which are selected as equal to 10,000 and 1, respectively. For each individual particle, after performing the movements given by the particle, P D is given as the final position of the robot. Meanwhile if there exist no obstacles between the last position obtained from the proposed optimization algorithm and P D , robot moves easily to the final position. Otherwise, obstacle stops robot from movement. It is further assumed that if the end effector is close to the desired point robot does not need to continue with the movements given by the proposed XOR BGSA and total time is considered to this point. When end effector is within obstacle free range of the desired position, the other movements given by the proposed XOR BGSA are bypassed and the desired position is used as the next step for robot.

Simulation Results
The proposed XOR BGSA with repository is iterated for 30 times in Matlab software connected to V-REP through API using 30 particles. The evolution of cost function during optimization is depicted in Figure 11 on a logarithmic scale. As can be seen from this figure, the best cost function value among the 30 particles starts from a high value of 3446 due to collisions and is reduced during optimization to 6.77, which is a collision free, time optimal solution to the problem which hits the desired point. This indicates that for the final collision free solution, the total time is equal to 6.77 s. The resulting movements are depicted with a red line in Figure 12, which can be inspected by human expert before final implementation on a real robot.
The 3D environment provided by simulation software provides complete visual demonstration of movements, which makes it is easy to perform the inspection and ensure suitability of the solution to real-world implementation before that happens. This can be an important safety feature within any working environment and such linking between the human and automated elements of production is a key factor in Industry 4.0 practice.

Conclusions
In this paper, an analysis of BGSA has been performed on mathematical update formulation for its acceleration term. Using such analysis, we show the fact that the existing acceleration and velocity term is not behaving as expected under certain conditions, and therefore conclude that it is required to modify the acceleration term of BGSA. We then propose the use of an XOR BGSA which benefits from repository. Analysis of this algorithm illustrates that the behaviour of acceleration and velocity terms are as expected. A repository is added to maintain best particle and provide a guideline for the movements of other particles. Complete statistical analysis is provided which demonstrate how the expected value of particles in each dimension converge to the particle selected from repository using the proposed version of XOR BGSA. To test the efficacy of the proposed optimization algorithm, a number of simulation studies are performed against different uni-modal and multimodal optimization problems. Simulation results indicate that the proposed modification outperforms other binary optimization algorithms namely BGSA [16], IBPSO [17], BPSO [15], BGWO [18] and BDA [19]. The proposed algorithm is used to solve a constrained maximization problem namely multi-knapsack problem. It is observed that the proposed algorithm outperforms BGSA [16], IBPSO [17], BPSO [15], BGWO [18] and BDA [19] for this optimization problem as well. Not only the mean value of results obtained during 30 times of optimization is better than other five optimization algorithms, but also the proposed algorithm benefits from a fast convergence rate. Moreover, repeatability of the proposed algorithm is better in most cases than in other investigated algorithms.
The conclusion made from the first part of simulation is that the proposed XOR BGSA with repository is an effective optimization algorithm to deal with binary optimization problems. Motivated by such results the algorithm is used for task planning of robot performing an assembly task. Industry 4.0 emerges the need for robust machine learning and artificial intelligence approaches to be used on the factory floor to determine responses to production changes resulting from continuous design changes. Assembly task planning and path planning in an industrial environment as two frequent tasks are considered. The assembly task planning problem scenario considered in this paper is to find the optimum order of bolts to be fastened to minimize the total distance travelled by robot. Using the proposed optimization algorithm the distance travelled by robot decreased by 16.7%. This is a practical application of the proposed optimization algorithm in flexible manufacturing. In addition, the proposed XOR BGSA with repository is utilized to assist in path planning in a manufacturing environment. Such application would allow the robot to find time optimal path while avoiding obstacles and moving to a desired position. This demonstrated the ability to successfully perform such an operation, while enabling human feedback on results for a final implementation of a real robot.