MSGM: A Markov Model Based Similarity Guide Matrix for Optimising Ordered Problems by Balanced-Evolution Genetic Algorithms

Where traditional genetic algorithms tend to prematurely converge on local optima, adaptive strategies aim to maintain a healthy level of population diversity by introducing randomness to the population. Often times this is done through adjusting control parameters according to diversity measurements. While these approaches introduce diversity, they do not aid in focusing or directing the search effort. Meanwhile, other works in the literature propose creating individuals designed to improve the population’s health and quality but their effectiveness is limited outside of general problems. This article proposes novel sequence-wise approach to designing and editing genotypes for ordered problems. A Markov model based similarity guide matrix (MSGM) is used to determine the relationships between gene nodes in order to produce new genotypes that focus on improving fitness and increasing population diversity. The proposed MSGM based approach is implemented in a balanced-evolution genetic algorithm framework in order to investigate its characteristics with encouraging results demonstrating its effectiveness when solving combinatorial ordered optimisation problems.


I. INTRODUCTION
In the field of biotechnology, genetic engineering is used to directly manipulate an organism's genes. This can be done through isolating and copying genetic material or by artificially synthesising the DNA itself. As genetic algorithms (GAs) are inspired by natural evolution, many studies have investigated concepts borrowed from biotechnology to improve its various operators and strategies. One research problem that stands to gain from this is the ability to direct a GA's search in a direction to improve its diversity, solution quality, or both. While a number of studies propose the creation of individuals from known qualities, these studies are largely limited to general optimisation problems [1]. Many of these approaches are incompatible or inappropriate for optimising ordered problems due to their constraints and their sequential nature of the encoded solution [2]- [4].
Combinatorial optimisation problems are often generalised forms of problems frequently encountered in the fields of The associate editor coordinating the review of this manuscript and approving it for publication was Radu-Emil Precup . operations and manufacturing. The process of solving these problems involves finding an optimal solution in a countably infinite set of possible solutions. Problems such as the Travelling Salesman Problem (TSP) and the Capacitated Vehicle Routing Problem (CVRP) are generalisations of problems in operations and logistics where the order in which the salesman visits the cities or a vehicle makes deliveries directly impact the cost of the operation [5]. These ordered problems are also often found in operations [6] and robotics [7]. Given a combinatorial ordered optimisation problem with N nodes, there are N ! solutions, thus finding the sequences and the relationship between each node in the different combinations is an NP-hard problem. Studies into the fitness landscape of problems like the TSP highlight the difficulties in solving these problems [8]- [10].
Given the difficulty of finding the optimal solution, heuristic approaches are often used to find high quality solutions within a reasonable time-frame. Innovations to adaptive GAs have enabled GAs to introduce a sense of intelligence in order to improve the efficiency and effectiveness of the search process. Continual improvements to these techniques FIGURE 1. Generalised framework for adaptive genetic algorithms, where population is a set of solutions for a problem which evolves through crossover and mutation. Common adaptive strategies include parameter tuning for the tournament size (t size ), crossover probability (p c ) and mutation probability (p m ).
enable GAs to further improve their abilities in finding optimal solutions to NP-hard ordered problems. Originally theorised as being a hill-valley landscape, recent studies into the fitness landscape of combinatorial optimisation problems have identified clusters or funnels in the fitness landscape [11], [12]. These clusters of local optima can be found throughout the fitness landscape with the difficulty of finding the global optima being related to the number of clusters in the fitness landscape and the size of the cluster that the global optima resides in [13], [14]. As understanding of the fitness landscape improves, more advanced techniques to direct the search of GAs are introduced. While several frameworks introduce techniques inspired by genetic engineering to help direct an adaptive GA's search pattern, the constraints and the characteristics of ordered problems have not been taken into consideration. In particular, the relationships between nodes, how these relationships contribute to both fitness and diversity, and the sequential nature of the problems themselves. With existing works demonstrating how a Markov model can effectively establish the relationships between gene nodes [15], [16], the following research questions can be raised: RQ1 How can we use Markov chains to introduce genome editing to direct the search? RQ2 What balance between fitness and diversity is needed to maintain a healthy level of diversity? With the success of driving genotype editing through the similarity guide matrix (SGM) in existing works [1], RQ1 aims to investigate the appropriateness of a Markov model in comparison to traditional approaches, such as measuring the Hamming distance. By measuring the likelihood of a node proceeding another node, a genotype can be generated according to the likelihood of its sequence rather than the likelihood of gene values appearing at each gene position. However, when considering directing the GA's search and population diversity, a purely diversity focused approach may not be the most appropriate. As the GA should aim to search for the clusters of local optima, RQ2 aims to investigate different strategies to discover and investigate the multiple funnels of the fitness landscape.
In this article, we present a method for generating new genotype solutions that balance between introducing diversity and improving fitness using a Markov model based similarity guide matrix (MSGM). The proposed method is implemented into an existing balanced-evolution genetic algorithm (BEGA) framework to demonstrate its effectiveness for optimising ordered problems and compared against the original framework for a range of ordered problems from the TSP and CVRP benchmark instances.
The remainder of this article is organised as follows: Section II highlights the existing works in the literature for adaptive GAs with Section III discussing the balanced-evolution genetic algorithm (BEGA) framework and the limitations of current methods. Section IV describes our proposed approach to adapting the similarity guide matrix (SGM) to generate new genotypes for the TSP and CVRP problems. A discussion and analysis of the results are included in Section V with concluding remarks in Section VI.

II. RELATED WORKS
Adaptive GAs aim to prevent premature convergence on sub-optimal solutions by introducing a degree of intelligence to the search process. These approaches often utilise feedback from population diversity measurements to adjust parameters in on online manner [17]. In doing so, adaptive GAs attempt to manage their focus between exploring the solution space and exploiting known solutions according to whether the population has converged too much or has too much diversity. A generalised approach to implementing an adaptive GA is shown in Fig. 1. Adaptive tuning for the tournament size (t size ) enables the GA to apply an appropriate amount of selective pressure while the crossover (p c ) and mutation (p m ) probabilities enable adaptive GAs to manage a balance between VOLUME 8, 2020 FIGURE 2. Different approaches to measuring the genotypic similarity between a pair of genotypes. local (exploitation) and global (exploration) search. The works in adaptive GAs can be classified as techniques and strategies for measuring the similarity between individuals, measuring the diversity of the population and mechanisms for controlling the amount of diversity [17]- [19].

A. GENOTYPIC SIMILARITY
Gene-wise measures generally consider the similarities and differences between genotypes according to the absolute positions of the genes. The two most common approaches are the Euclidean and Hamming distances [18], [19]. While several works have applied these to solving ordered problems [20], [21], they do not consider the relationship between genes when measuring population diversity. The limitations of these approaches are further highlighted by works implementing approaches that focus on these relationship between nodes as shown in Fig. 2a.
The broken pairs [22] approach (Fig. 2d) considers the relationships between each pair of neighbouring genes and measures the differences between two genotypes as the number of pairs that have been separated. Numerous works [23], [24] demonstrate the effectiveness and improvements that this approach introduces in comparison to a gene-wise approach. However, the considerations of the sequence-wise nature of the problem is limited to pairs of genes.
In order to consider the relationship between a wider range of genes, Nagata and Kobayashi [25] proposed the use of a Markov model to measure the similarity between genotypes in the population and demonstrates how increasing the scope of the relationship between genes can improve on a GA's ability to maintain diversity. This is further demonstrated with the variable-order Marokov model [15], [16]. Another method for measuring the sequence-wise diversity of a population is the use of the longest common subsequence length (LCS) distance [2]- [4]. Similar to the Hamming distance, it measures the number of gene nodes that share a common non-contiguous subsequence between two genotypes as shown in Fig. 2c. This is more effective than measuring the contiguous subsequence (Fig. 2b) as a non-contiguous measurement would also include contiguous subsequences [2].
While there are different approaches to measuring the similarities and differences between two genotypes, each one has different costs and benefits when applied to ordered problems. How they are applied to measuring population diversity greatly affects their contribution to an adaptive GA's overall performance.

B. POPULATION DIVERSITY MEASUREMENTS
While research into determining the similarities and differences between genotypes is an active and ongoing field, a wide range of strategies for applying these measurements to maintain population diversity exists in the literature [17]. Early works aimed to measure the Hamming best and worst performing solutions [26] in order to minimise the computational costs. However, it was limited in that it could not provide feedback on the state of the population.
While hardware was a limiting factor to the practicality in measuring the diversity from all individuals in the population, Shimodaira [27] proposed a compromise by measuring the Hamming distance between the best, worst and a selection of the population. While this increased its scope, it was still limited by design with later works [28], [29] measuring the distance between all individuals. Another common strategy is to measure the distance between a reference individual and the population. These methods would compute a solution that represents a central point in the code space. While this often leads to solutions that are not valid within the solution space, they still demonstrate a utility in measuring diversity.
Mc Ginley et al. [30] uses a genotype consisting of average values at each position as a reference point to measure the Euclidean distance of the population. As this results in a non-integer value, approaches that measure the Hamming distances [1] or LCS distance [2]- [4] must come up with other methods such as the use of a mode or median genotype. Another method includes measuring the entropy rate of the population [15] that measured the probability of a sequence of nodes occurring in the population. With many strategies available for measuring population diversity, how diversity is introduced into the population is an active and on-going research area.

C. DIVERSITY CONTROL MECHANISMS
Once feedback on the amount of diversity has been received, an adaptive GA must then activate its mechanisms for introducing and controlling diversity. Early methods included re-initialising the population [31] to reboot the population as well as adaptive parameter controls [32]- [34].
Other approaches included adaptive population sizes [35], [36] and fixed populations with dynamic subpopulations [30]. The latter allowed for an adaptive GA to balance its focus between exploration and exploitation without increasing the computational effort. Where these works focused on general, unordered problems, their performance was limited when applied to solving within the constraints of ordered problems [2], [3]. Zhang et al. [1] proposed a method for generating new individuals that either converged or diverged from the current population in their Balanced-Evolution Genetic Algorithm (BEGA) framework. Using the Hamming distance to measure the population's diversity allows for BEGA to determine how much the population needs to converge or diverge. Positive and negative perturbations of the similarity guide matrix (SGM) were used to generate individuals for exploitation and exploration subpopulations. This approach was novel in that it influenced the direction of the search area but was limited in its application due to its gene-wise approach to both measuring and maintaining diversity. While works have demonstrated how adaptive GAs can improve their performance by adapting to sequence-based approaches to diversity maintenance [2], [4], developing sequence-based approaches to directing the search using sequence-based genetic engineering approaches is an open research problem.

III. BEGA FRAMEWORK
The balanced-evolution genetic algorithm (BEGA) framework [1] was originally proposed for solving unordered problems such as the Knapsack problem and general mathematical problems. The framework balances the exploration and exploitation responsibilities of its search through two subpopulations. These subpopulations follow similar mechanics to one another but with one crucial difference. Unlike traditional GAs, the BEGA framework does not select two individuals for crossover but creates a temporary population of candidate individuals that are used to crossover with the real population. These candidate individuals are designed specifically to either encourage convergence or introduce diversity depending on which subpopulation they are designed for. This mechanism is controlled by the similarity guide matrix (SGM) and the linear diversity index (LDI). Table 1 summarises the symbols used in this article.

A. SIMILARITY GUIDE MATRIX AND LINEAR DIVERSITY INDEX
The SGM creates a two dimensional matrix that expresses the probability distribution of the gene values at each gene position. This can be seen in Fig. 3a. In the case of an ordered VOLUME 8, 2020 problem such as the TSP, the matrix is a size of N 2 and can be expressed as Eqs. 1-2.
The two axis of the SGM (S) corresponds to the gene positions (Eq. 1) and gene values or nodes (Eq. 2). Each element in the matrix (s ij ) represents the proportion of individuals in the population that has the node value i in the jth position. From the SGM, negative and positive perturbations are computed for the exploitation and exploration subpopulations. The differences in these perturbations are demonstrated in Fig. 3b. The control amplitude (CA) determines the degree of change in the perturbations. This is calculated using the LDI which is calculated as the average Hamming distance between a mode genotype and the rest of the population. This mode genotype can be computed from the SGM as the gene values at each position with the highest probability. The process for calculating the LDI is demonstrated in Eqs. 3-4.
The normalised diversity of each individual (D p ) is calculated as the Hamming distance between individuals in the population (X p ) and the mode geneotype reference point (X r ) normalised against the genotype length (N ). The LDI (D l ) is then calculated as the average, normalised Hamming distance of the population (P). The D l is used in conjunction with the diversity shift limit (D sl ) to determine what stage of the evolution BEGA is in and the level of diversity that to be injected into the population. During the first stage where D l > D sl , D l is used as the CA value as shown in Eq. 5.
However, as BEGA converges on a optima, the population's diversity (D l ) can become a very small value which prevents the GA from being able to effectively maintain diversity. D sl is designed to provide a minimum degree of diversity maintenance even in a maximally converged population with the authors recommending a value of 0.075.

B. COMPUTATION OF PERTURBATIONS
With both S and D l , the perturbations for the exploration and exploitation subpopulations can be computed. The process for computing the negative perturbation is demonstrated in Algorithm 1. After the control amplitude (C A ) is calculated, the new perturbation matrix (M ) must be populated with new values. For each vector (j) in S, the vector has a chance to be inherited directly from S without any changes. This is dependent on the state of the population's search process.

Algorithm 1 Computing Negative Perturbation
Input: S: SGM, N : Genotype length, CA: Control amplitude Output: M : Perturbation matrix If vector j is selected to undergo changes in its perturbation, then the index (k) of the maximum value (s kj ) of the vector S vj is selected. A randomly selected index (l) from S vj is selected as a pair for the the perturbation where k = l. The negative perturbation of S reduces the maximum S vj value by r and increases the value of the randomly selected member by the same amount with all other members of S vj are directly inherited by M . This spreads the probability distribution for values at each gene position to decrease the likelihood of a common genotype being generated.
The positive perturbation for exploitation works similarly but with k being the index of the minimum value of the vector S vj . Once the positive and negative perturbations of S are computed, candidate individuals can be generated according to the probability of values at each gene position (Fig. 3c).

C. EXPLORATION AND EXPLOITATION
The first step in generating new individuals is to produce a set of candidate individuals for both the exploration and exploitation subpopulations. Candidates for both methods are generated from their respective perturbation matrix and is demonstrated in Alg. 2. M is the SGM with negative or positive perturbations applied. T is the computed threshold matrix that allows for a uniformly random number (rand()) between 0 and 1 to select node i for the jth position in the candidate genotype Q.
The goal of the exploration subpopulation is to inject diversity into the population in order to explore new areas of the fitness landscape. Creating new genotypes for the exploration subpopulation is outlined in Alg. 3. The SGM with negative perturbations (M ) is used to generate a new candidate genotype Q. For crossover, the original BEGA framework uses uniform crossover between the candidate genotype and the reference genotype. Uniform mutation is applied to the new genotype with a probability p m = (1.0 + m s × CA)/N . The exploitation method is similar to the exploration method except the positive perturbation matrix is used to generate Q and constant mutation probability is used (p m = 1/N ).

D. FRAMEWORK: PUTTING IT ALL TOGETHER
The BEGA framework manages two sub-populations in order to balance between exploration and exploitation and is demonstrated in Fig. 4. A sub-population of elite individuals (P e ) is selected with a minimum Hamming distance (m d ) between each individual as the exploitation sub-population. Individuals that are not selected for the

E. LIMITATIONS
As BEGA relies on a gene-wise approaches for the SGM, LDI, the positive and negative perturbations, and generating the candidate individuals, there are significant limitations when applied to ordered problems such as the TSP. This can immediately be seen in the candidate individuals themselves where the perturbation matrices can result in solutions that do not adhere to the problems' constraints.
While duplicate gene values can be replaced with missing gene values after the individuals are created, this can reduce the effectiveness of the diversity balance mechanisms introduced by BEGA [2]. Another limitation is its focus on absolute gene positions rather than the relationships between gene values and their positions in the sequence. This limits the ability of both the SGM and LDI to express the diversity to a position-by-position approach. Without being able to accurately measure and monitor the state of the population in its evolutionary process, BEGA is unable to effectively determine when diversity is needed or by how much. These two limitations results in BEGA having difficulty both monitoring and maintaining population diversity in an effective manner when solving ordered problems.

IV. OUR APPROACH
To generate effective candidate genotypes for ordered problems, the following constraints must be maintained: C1 Each node in the problem must appear in the sequence and cannot be repeated, and C2 A route can start at any node but must return to its original node or depot for the TSP and CVRP respectively. C1 requires genotypes to be generated where there are N unique nodes in the sequence where the values range from [1, . . . , N ], inclusively. The similarity guide matrix used by the BEGA framework can result in invalid genotypes where the duplicated nodes must be replaced with missing nodes retrospectively. A strategy for producing genotypes for these ordered problems must prevent duplicated or missing nodes in order to produce meaningful sequence-based genotypes. C2 describes the cyclical nature of ordered problems where the sequence of nodes is more important than their absolute positions. Thus, a method for designing and producing new genotypes must consider the ordered nature of sequence-based problems.
To model the above constraints of the ordered problems in the similarity guide matrix (SGM), we propose three new approaches to compute the SGM based on the Markov model, called the Markov model based Similarity Guide Matrix (MSGM). The first method uses the Markov model's transition matrix as the MSGM. The second method measures the average fitness of the subsequence of genes and can be considered as the contribution of a subsequence of nodes to a genotype's fitness. The final method is a hybrid approach that considers both the probability of a subsequence occurring and its contribution to a genotype's fitness.
Where the original BEGA framework uses a gene-wise approach (Hamming distance) to measure diversity in Eq. 4, in our approach we use the LCS distance based genotypic similarity measurement explored in our previous works [2]- [4] to measure diversity as shown in Eq. 6.

A. MARKOV MODEL BASED SGM
The BEGA framework utilises the SGM to express the distribution of a population or subpopulation in the solution space. This is then used as a reference point in the coding space order to measure the density of the solutions around that point and to produce candidate individuals. However, as highlighted in Fig. 3, this gene-wise approach can result in individuals being produced that do not consider the constraints and characteristics of ordered problems.
To address the aforementioned issues, we propose the use of a Markov model's transition matrix as the Markov model based SGM (MSGM). Where the original SGM had axes of coding position and value, a sequence-based approach considers the matrix axes as nodes from and to. If we consider the likelihood of a node preceding another node in a solution as a Markov chain with one order, the relationship can be modelled as in Eq. 7. P r (q 1 , q 2 , . . . , q N ) P r denotes the probability of the nodes appearing in a given sequence. However, this does not accurately portray the relationships between solution sequences and the fitness landscape in ordered problems. As a population begins to converge on clusters of optima, each cluster will being to show characteristics in its genotype sequence that are shared with nearby solutions. When computing the MSGM, we can calculate the probability of a node proceeding another node to build a sequence-wise similarity guide matrix. We can express the MSGM (S M ) as a two dimensional matrix shown in Eq. 9. However, each member of S M vj represents the probability of node i following node j as shown in Eq. 10.
S M vj = [P r (1, j), . . . , P r (i, j), . . . , P r (N , j)] 1 < (10) Where the original SGM expresses on the absolute positions of each gene node, the MSGM expresses on the probability of transition between nodes. This allows the MSGM to describe the sequence of nodes rather than their positions. The sum for each S M vj will always be 1 in ordered problems like the TSP and CVRP as each node must appear once in the sequence. Furthermore, for these problems, sequences can wrap around the genotype and remain similar to one another. For example, a sequence of [1,2,3,4] should be considered similar to [3,4,1,2] where both have the same P r (4, 1) and P r (2,3).
11) a ij represents the probability that node i follows node j given the assumptions that the i and j are the values from N nodes and that it is a steady state system. Further more, as s M ij ≥ 0 and the genotype cannot repeat nodes in its sequence, it must eliminate the probability of j following i if j already appears in the genotype sequence. This can result  The sum of all the probabilities of nodes proceeding i must be N j=1 s M ij = 1. S M can be considered as a two dimensional matrix that expresses the distribution of the sequence of gene nodes. This can be seen in Fig. 6 where a first order Markov model can be used to calculate the subsequences of the population in Fig. 6a into the S M (Fig. 6b). The vertical axis represents the current gene node while the horizontal axis represents the probabilities of the nodes that follow.
The process of creating a new genotype from the S M is demonstrated in Fig. 5. An initial gene node can be selected according to its probability as a starting point with the proceeding node being selected according to its probability. Once a node has been added to the genotype, it is removed from the possible options in the vertical axis. This ensures that the constraints of the ordered problem are maintained. With the MSGM, the transition matrix can be used as the similarity guide and is shown in Algorithm 4.
For each individual in the population P, the integer genotype sequence is expressed as vector X . S C is the transition count matrix to store the number of genotypes in P where node i following node j. Thus, for each node in X , we increment the count at S C [i] [j]. Once the sum of transitions between all nodes for all individuals in P have been added to S C , each value in the matrix is normalised by N to give a value ranging from 0 to 1. This is stored in the Markov Guide matrix (S M ) and returned as the transition matrix.
As the transition matrix describes the probabilities of one node following another, producing new genotypes is a stochastic process. As the population converges on a set of solutions, common subsequences in the population increases the probability of these subsequences appearing in genotypes produced by the MSGM. However, the stochastic nature of the Markov chains allows for the MSGM to introduce minor mutations to the overall sequence and can assist in the process for local search. VOLUME

B. FITNESS BASED MSGM
Where the MSGM approach considers areas of interests in the solution space as where the solutions are converging, it does not directly consider the fitness of the solutions. This is particularly important when taking into consideration the complexities in navigating a fitness landscape with many clusters of local optima.
The fitness-based MSGM computes the average fitness contribution of the gene sequence [i, j] from the population which results in candidate genotypes more likely to inherit the fittest subsequences instead of the most common. This process is demonstrated in Algorithm 5 where S F is a two dimensional matrix which contains the fitness distribution of node i following node j as opposed to the transition probabilities in S M . Similar to Algorithm 4, for each solution in P, the genotype sequence is stored as X . As this is a fitness based transition matrix, the summation of the fitness is required later and is stored in f . For each node in the genotype X , the fitness of X is added to the matrix S F [i][j] where the sequence transitions from node j to i. Where the S M in Algorithm 4 is the summation of transitions between nodes, S F is the summation of the fitness of genotypes where a transition between nodes j and i exists. The fitness guide matrix (S F ) consists of the values of S F normalised against the fitness sum of the individuals in P to give a value between 0 and 1 and illustrates the contribution of each pair of nodes to the average fitness.

C. BALANCED MSGM
A limitation of S F is that it considers a lower fitness individual as being more diverse than an individual with a higher fitness and can be considered as being similar to a phenotypic approach to measuring diversity. Where as the limitation of S M is its difficulty in conducting intensive search in multiple areas of interest. This highlights the necessity for a strategy that balances both converging towards good solutions and exploring new areas of the solution space. We propose a hybrid approach, called balanced MSGM that balances between genotypic diversity and fitness contribution. This allows BEGA to quickly converge on areas where known, highly fit solutions exists while also explore new areas of the solution space.
The balanced MSGM is given in Algorithm 6. Similar to Algorithms 4 and 5, the process iterates for all individuals x in P with the genotype sequence of x being expressed as X . For each node in the genotype, the transition from node j to i is added to the count matrix S C [i][j] while the fitness of the individual is added to the fitness sum matrix S F [i][j]. Once S F has accounted for all individuals in the population, the balanced guide matrix (S B ) is calculated as the average value of the normalised S F and S M values.

V. EXPERIMENTS
Here, we evaluate the effectiveness of our proposed Markov model based similarity guide matrix (MSGM) in optimising ordered problems by balanced-evolution genetic algorithms.

A. SETUP, ALGORITHMS AND PARAMETERS
This section describes the setup, algorithms, their parameters and evaluation metrics used in our experiments.

1) ENVIRONMENTAL SETUP
All algorithms are implemented using C#/.NET Core with experiments run on Windows 10 machines with AMD Ryzen 2600x CPU and 32GB of main memory. Each framework algorithm ran each instance until 20,000 generations with no improvements in the solution quality was reached. This was repeated for 20 runs. The C# implementations of the algorithms are made available online [37].

2) BENCHMARK ORDERED PROBLEMS AND OPERATORS
Instances of ordered problems were selected from the TSP benchmark library 1 according to their range in problem size. The new set of benchmark instances 2 proposed by Uchoa et al. [5] were selected for the CVRP tests. Due to the ordered nature of these problems, the modified ordered crossover operator and partially shuffled mutation operator were selected for their performance demonstrated in empirical tests [38].

3) ALGORITHMS AND THEIR PARAMETERS
The following algorithms have been implemented to conduct experiments in this article.
• GA: This is a standard genetic algorithm with no adaptive features.
• BEGA: The original BEGA framework (as described in Section III) is implemented in C# according to the authors' work [1]. This ensures for consistent comparisons between each implementation.
• BEGA MSGM: BEGA framework with LCS distance being used to calculate population diversity (D p ) as in Eq. 6. We also use MSGM (Alg. 4) instead of SGM.
• BEGA MSGM Fitness: BEGA framework with LCS distance being used to calculate population diversity (D p ) as in Eq. 6. We also use MSGM Fitness (Alg. 5) instead of SGM.
• BEGA MSGM Balanced: BEGA framework with LCS distance being used to calculate population diversity (D p ) as in Eq. 6. We also use the MSGM Balanced (Alg. 6) instead of SGM. The BEGA framework parameter settings were kept consistent between each variant as recommended by Zhang et al. [1]. The original BEGA framework was implemented with its original operators. The GA, BEGA MSGM, BEGA MSGM Fitness and BEGA MSGM Balanced were implemented with the Modified Ordered Crossover (MOX) operator and the Partially Shuffled Mutation (PSM) operator due to their performance in ordered problems [38].
Population size (|P|) for each algorithm is set to 90. The size of the elite sub-population |P e | = 15 with a minimum Hamming distance between individuals being m d = 3 for the original BEGA. The BEGA MSGM variants require a minimum LCS distance m d = 3, shift limit of diversity D sl = 0.075. For BEGA and all of its variants, the multiplier factor for the mutation operator m s = 5.

4) EVALUATION METRICS
To evaluate the algorithms, the best known solution (BKS) and problem size (N ) has been included in the problem instance descriptions. To indicate the relative difference between each GA approach, the average error (Avg. Err) between the BKS and the solution found is reported with the coefficient of variation (C v ) demonstrating its consistency. BKS Found is reported to indicate the number of runs in which the GA approach found the BKS (out of 20 runs). The average diversity (Avg. Diversity) is reported as the average of the diversity of the last generation for each run. The diversity of the final population is measured as the average LCS distance normalised against the length of the genotype.

B. PERFORMANCE EVALUATION
Tables 2-3 demonstrate the abilities of each of the GAs. The fitness is represented as the cost of the solutions where a lower value indicates a higher quality solution.

1) EFFECT OF PROBLEM SIZE
While the benchmark GA performs well on smaller problems, it can be seen to quickly converge on local optima as the problem size increases. The original BEGA implementation, while a significant improvement over the GA, also prematurely converges. However, applying a sequence-wise MSGM greatly improves the performance of BEGA in its ability to find better quality solutions for both smaller and larger problems. This improvement is also demonstrated when more constraints are introduced as the CVRP (Table 3). The fitness based MSGM (BEGA MSGM Fitness) appears perform worse than BEGA MSGM in both solution quality and consistency. While MSGM Fitness is effective at exploiting its found solutions, it can be seen to struggle with effectively searching the solution space for other peaks in the fitness landscape. However, by combining both of these approaches, the hybrid approach (MSGM Balanced) increases BEGA's ability to find better quality solutions and improve on its consistency in finding these solutions.

2) EXPLORATION VS EXPLOITATION
The three BEGA variants with sequence-wise approaches offer very significant improvements over the GA and many improvements over the original BEGA. However, the improvements of the MSGM Balanced approach over both the BEGA MSGM and BEGA MSGM Fitness offer more insight into the contributions of the two. The characteristics of each of the BEGA MSGM variants can be seen in Fig. 7.
While BEGA MSGM Fitness variant performs better in earlier generations, BEGA MSGM is capable of finding better solutions for more generations. This pattern is made clearer in the Figs. 7c-7d with the CVRP instance X-n10001-k43. With further constraints added by the CVRP, both BEGA MSGM Fitness and BEGA MSGM Balanced reach their completion criteria before BEGA MSGM. While BEGA MSGM can be seen to evolve for longer than BEGA MSGM Fitness, Tables 2-3 demonstrate that having the highest diversity does not directly result in the best solution. BEGA MSGM Balanced maintains a healthy level of diversity to evolve for longer and produces better solution quality. VOLUME 8, 2020 TABLE 2. Performance for TSP instances highlighting the average error, coefficient of variation (C v ) in solution quality, the number of runs a GA framework found the best known solution, and the average diversity of the last population of each run.

C. DIVERSITY MAINTENANCE
In Tables 2-3, the average diversity highlights each GA framework's ability to maintain sequence-wise diversity. It is measured as the average LCS distance of the population at the completion of the run. As the value is normalised against the genotype length, a maximally diverse population approaches a value of 1 while 0 indicates maximal convergence. The capabilities for each GA framework in maintaining sequence-wise diversity is further highlighted in Fig. 8.

1) EFFECT OF PROBLEM SIZE AND CONSTRAINTS
In the TSP instances, the basic GA and original BEGA implementations can be seen to struggle to maintain a diverse population with a trend towards greater convergence with larger problem sizes. This suggests that while BEGA focuses on maintaining gene-wise diversity, almost half the sequences are common to all the solutions in the population. BEGA MSGM maintains a higher level of diversity regardless of the problem size. The BEGA framework  [5] set from CVRP-Lib with instances grouped by N size highlighting the average error, average coefficient of variation (C v ), number of instances where the GA framework found the best known solution in its run and the average diversity of the final population for each run. supports exploitation and exploration subpopulations where the exploitation subpopulation naturally converges on similar solutions. The higher levels of diversity in the BEGA MSGM Balanced variant suggests that while the exploration subpopulation is sufficiently diverse, the exploitation subpopulation also has a higher level of diversity allowing it to exploit multiple known peaks in the fitness landscape.
In Table 3, the effects of the additional constraints of the CVRP can be seen on the diminished population diversity. In particular, the ability to maintain population diversity is reduced in the BEGA MSGM Fitness variant as the problem size increases while both the BEGA MSGM and BEGA MSGM Balanced variant are able to maintain a more consistent level of population diversity.

2) EXPLORATION VS EXPLOITATION
In Fig. 8a both the basic GA and BEGA can be seen to converge very early in the search process. While the basic GA has no diversity maintenance processes, the original BEGA framework uses a gene-wise approach to measuring and maintaining diversity. This can result in a higher level of sequence-wise similarity in the population. BEGA MSGM has the highest level of overall diversity due to its bias towards exploration as seen in Fig. 8b. As the framework aims to generate sequences that maximises the differences in sequences, this can naturally lead to an effective exploration process. The BEGA MSGM Fitness approach focuses on the exploitation of the fittest individuals. This is done through the fitness based guide matrix where the fitness of each sequence influences the likelihood of a subsequence appearing in a genotype. This results in the population converging on certain subsequences and lowering the overall sequence-wise diversity. Furthermore, the three BEGA MSGM variants can be seen to have a high level of fluctuation in population diversity in the early stages of their evolution. As the population begins to converge, the fluctuations in diversity significantly reduces. This can be seen as the second stage where the shift limit (D sl ) is used to maintain diversity as the LDI (D l ) is too small due to the level of convergence. What should be noted here is the degree of the fluctuations in diversity, particularly with the constraints of the CVRP in Figs. 8c-8d. BEGA MSGM and BEGA MSGM Balanced both display similar fluctuations between the peaks and valleys in diversity that suggests that diversity is introduced in a very acute and effective manner. However, the BEGA MSGM Fitness variant is not as effective.
While the BEGA MSGM Balanced variant has less diversity than the BEGA MSGM, it maintains diversity levels closer to that of the BEGA MSGM. This suggests that the mechanisms for managing the subpopulations are better equipped to maintain a high level of diversity for exploration but a smaller, more intense subpopulation for exploitation. When this is applied to solving ordered problems, BEGA MSGM Balanced is able to maintain a healthy level of diversity while also being able to effectively find good quality solutions. These characteristics are what likely enable the BEGA MSGM Balanced variant to outperform the other benchmark GAs.

D. DISCUSSION
The main objective of the BEGA framework is to balance between intensive local search (exploitation) and diverse global search (exploration). By a novel implementation of the similarity guide matrix, BEGA designs and produces genotypes to encourage convergence for the local search subpopulation and introduces diversity to the global search to improve coverage of the fitness landscape.
In order to determine the statistical significance of any improvements demonstrated by the proposed methods, two-sample z-tests were conducted between a benchmark GA, the original BEGA implementation, as well as the MSGM, MSGM Fitness and MSGM Balanced implementations. The p values for the z-tests were 0.05 or less being considered a significant difference and 0.01 or less indicating a very significant difference. Table 4 highlights the statistical significance of the solutions found by BEGA MSGM Balanced.
In general, the MSGM approaches demonstrate a significant improvement in the BEGA framework's ability to find good quality solutions by adapting the gene-wise approaches to sequence-wise approaches. By implementing a Markov model to measure the probability of node transitions, BEGA is able to produce genotypes that better reflect the solution space and also reflect the ordered nature of the problems. By balancing between sequence-wise diversity and fitness contribution, our experiments demonstrate how our MSGM approaches are able to engineer genotypes for ordered problems and outperform the original BEGA framework.

VI. CONCLUSION AND FUTURE WORK
In this study we propose three methods for computing a Markov model based similarity guide matrix (MSGM) that takes into consideration the sequence of nodes for optimising ordered problems by the balanced-evolution genetic algorithms (BEGAs). We demonstrate how the MSGM can be used to genetically design and engineer genotypes to encourage convergence on good performing subsequences or encourage search in unexplored areas of the fitness landscape when implemented in an existing BEGA framework. Our experimental results highlights how a sequence based MSGM can introduce bias towards exploration while a fitness based MSGM biases towards exploitation. A balance between these two MSGM approaches demonstrated the ability to produce the best candidates for maintaining a healthy level of population diversity. Future research into improving the balance between exploration and exploitation for new genotypes should consider monitoring and maintaining this balance in an online manner.