Adaptive multiple crossover genetic algorithm to solve workforce scheduling and routing problem

The workforce scheduling and routing problem refers to the assignment of personnel to visits, across various geographical locations. Solving this problem demands tackling numerous scheduling and routing constraints while aiming to minimise the operational cost. One of the main obstacles in designing a genetic algorithm for this problem is selecting the best set of operators that enable better GA performance. This paper presents a novel adaptive multiple crossover genetic algorithm to tackle the combined setting of scheduling and routing problems. A mix of problem-specific and traditional crossovers are evaluated by using an online learning process to measure the operator’s effectiveness. Best operators are given high application rates and low rates are given to the worse. Application rates are dynamically adjusted according to the learning outcomes in a non-stationary environment. Experimental results show that the combined performances of all the operators were better than using only one operator used in isolation. This study provided an important opportunity to advance the understanding of the best method to use crossover operators for this highly-constrained optimisation problem effectively.


Introduction
The workforce scheduling and routing problem (WSRP) is described as the assignment of personnel to visits, requested by customers, across different geographical locations.
A type of WSRP arises in home health care where nurses and care workers should be assigned to visit patients in their homes to carry out some tasks, e.g.administering medication, monitoring serious illness and unstable health status and injections (Castillo-Salazar et al. 2015).
This problem combines scheduling and routing problems, both of which are known to be NP-Hard (Lassaigne and Rougemont 2012).The scheduling aspect of the problem assigns personnel to visits in order to fulfil work demands and other requirements.The routing aspect of the problem consists of generating routes for the workers to service customers across various locations within given time-windows.The objective is to minimise operational costs while attending the additional requirements expressed by customers, workers and the business.
Genetic Algorithms (GAs) have been shown to be effective approaches to find solutions for problems combining scheduling and routing where exact methods are less effective, e.g.Cowling et al. (2006), Mutingi and Mbohwa (2014) and Algethami et al. (2017).
A baseline GA was proposed by Algethami et al. (2016) that identified the best set of operators and parameters for each instance of WSRP.Despite its success in obtaining good solutions, the GA performance was limited by a computationally expensive parameter tuning method.Thus, an adaptive parameter control approach was proposed by Algethami and Landa-Silva (2017).In the later, the method maintained population diversity to enhance the performance of the GA.On the other hand, in the study by Algethami et al. (2016), there was no evidence that the performance of one crossover operator was superior to a group of operators during the run.Therefore, a random crossover exchange was also proposed by Algethami and Landa-Silva (2017), in addition to the parameter tuning method.Further improvements were suggested, especially on operators selection.
This study proposes an adaptive multiple crossover GA that uses a learning process to enhance the performance.The idea is to use an adaptive allocation rules on a mix of problem-specific and traditional crossovers, which are evaluated to measure the operator's effectiveness.Based on their efficiency, a roulette wheel concept is used to select an operator.The contribution of this study is to present a dynamic operator selection, as the search progress, as part of the multiple crossover framework, to reflect the crossovers behaviour in one iteration, according to the learning outcomes in a nonstationary environment.
This study claims that using a multiple crossover framework can enhance the GA efficiency when tackling WSRP scenarios.To the best of our knowledge, this adaptive multiple crossover GA has not been investigated for WSRP.To this end, the specific objectives of this study are: • To present an adaptive multiple crossover GA that aims to enhance the baseline GA efficiency.• To analyse the impact of crossover operators collaboration on the quality of solutions as well as algorithm speed.• To carry out illustrative computational tests to show the utility of the proposed GA approach by comparing the performance of the various genetic operators consid-ered.A total of 42 problem instances from six different real-world home health care scenarios were used in the experimental study.
In what follows, Sect. 2 reviews related work.Section 3 describes the WSRP, its formulation and the instances used in this paper.Section 4 outlines the proposed algorithm.Section 5 presents the experimental study and discusses the obtained results.The paper is then concluded in Sect.6.

Related work
Recent research on the WSRP considered here is reviewed next.A mixed integer programming (MIP) with decomposition method (Laesanklang and Landa-Silva 2016) required considerable computation time (up to several hours) to solve larger problem instances with hundreds of tasks, indicating the need for faster solution methods.A Variable Neighbourhood Search (VNS) algorithm using problem-specific neighbourhood heuristics was presented in Pinheiro et al. (2016).The VNS obtained high-quality solutions and in less computation time for the same set of problem instances used in Laesanklang and Landa-Silva (2016).A number of studies have applied GAs to real-world problems where scheduling and routing are combined.Examples include Cowling et al. (2006) and Mutingi and Mbohwa (2014).In those works, the focus has been on algorithm design to obtain good solutions.
An investigation was presented by Algethami et al. (2017) comparing various genetic operators within two simple GAs, to tackle the subject problem.A more efficient GA was proposed by Algethami et al. (2016) with tuned parameters and using a customised solution representation to maintain feasibility of solutions.To enhance the GA efficiency, adaptive concepts were used by Algethami and Landa-Silva (2017) in addition to random crossover exchange.
Using a learning method, to adjust parameter values has proven to enhance the baseline GA performance (Algethami and Landa-Silva 2017).Nevertheless, selecting a random crossover, without any prior knowledge of its performance, was a straightforward procedure.It is still not clear which operator was the best for each problem set, during the run.Thus, this paper proposes a learning-based multiple crossover framework to tackle WSRP.
According to Maturana and Saubion (2008), multiple operators GAs has been used on real-world applications to benefit from the different performance of synchronised operators.Most researchers have utilised a group of operators, crossovers (Consoli and Yao 2014) or mutations (Contreras-Bolton et al. 2016), as part of the algorithm.For example, the study by Puljić and Manger (2013) that has investigated mixing eight crossover operators on routing problems.One crossover was selected at random, with all the operators having the same probability.This operator is then applied at the current iteration.The study suggested multiple crossovers method obtained better results than a traditional GA.Another study that proposed a multiple operator algorithm for routing problems was presented by Osaba et al. (2014).The aim was to test it for a variety of combinations of operators.As a result, the algorithm was proven to be efficient to solve routing problems, such as the Traveling Salesman Problem, the Capacitated Vehicle Routing problem, Vehicle Routing with Backhauls and Multi-Depot Vehicle Routing Problem.
Two main operations are applied when using multiple operators, operators selection and evaluation.Adaptive operator selection (AOS) is identified as the online adjustments of the crossover function (Belluz et al. 2015), where all crossovers are used as one operator.However, each crossover has an application rate, that is relative to the crossover performance.The better the crossover, the higher the chances to be implemented more than a poor crossover.Thus, providing a variety of performances in one iteration.The diverse set of operators explore the search space differently, thus, more areas are explored.
To measures a crossover effectiveness (performance), an operator evaluation process is used to analyse the impact of an operator application on the current search (Li et al. 2014).This method includes giving some reward to an operator according to the operator impact on the search.The next crossover is selected based on rewarded value.The best operator is scored higher in the credit registry and therefore this crossover is selected more.
In general, fitness improvement was used to measure a crossover performance, i.e. the better the new solutions, generated by an operator, the more chance of an operator to be applied.However, more performance measurements were explored and proven to provide good results, such as distance-based measurement (Consoli and Yao 2014) and an operator execution-time (Maturana and Saubion 2008).
Many considerations were taken into account when designing a multiple crossover GA, including the evaluation measurement used and the number of iterations required to obtain sufficient information on the operators' performance.
Existing research has focused on using a fixed number of iterations to determine "the best" operator.However, an AOS might converge to the best performing operator early in the search.Not to mention that different operators provide different results at different stages of the search.Thus, using preliminary information can have limited flexibility and leads up to a loss of performance.This observation was proved by Belluz et al. (2015), in which the use of a "dynamic" approach was established to be better than using a static one.
The intended contribution of this study is to present an efficient adaptive multiple crossover genetic algorithm (AMCAGA), as an adaptive variant of the GA presented in Algethami and Landa-Silva (2017).
The proposed GA uses allocation rules on a mix of problem-specific and traditional crossovers.Best operators are given high application rates, and low rates are given to the worse operators.Based on the learning outcomes, application rates are dynamically modified, to reflect the crossovers' behaviour at the current iteration.A roulette wheel concept is also used when selecting an operator, as an extension of the works by Tuson (1995) and Thierens (2005).However, in this study, the selection process is dynamic to provide a better reflection of the operator's performance, as the search progress.This study claims that using an active change of the application rates can enhance the GA efficiency when tackling WSRP scenarios.The objective is to show that the proposed method outperforms the GAs presented in Algethami et al. (2016) and Algethami and Landa-Silva (2017) in terms of results quality and run-time.

Problem description
A WSRP solution S is a daily plan of visits, i.e. a set of workers W = {w 1 , w 2 , . . ., w |W | } assigned to perform a set of tasks T = {t 1 , t 2 , . . ., t |V | } for customers at different locations, where V is the number of visit assignments.The assignment of a worker w to travel to a customer location in order to perform a task is called a visit.
Several features have been identified as important in solutions to WSRP scenarios, such as distance travelled and customers' and workers' requirements and preferences (Mankowska et al. 2014).Thus, a good quality solution should have low operational cost as well as all visits are assigned while satisfying the existing requirements.
For example, an illustration of a possible plan for a WSRP instance is presented in Fig. 1.In this example, 3 workers (w 1 , w 2 , w 3 ) are assigned to 15 different visits, located at different geographical areas.A path is plotted as a dotted line with a different colour for each worker, i.e., worker 1 is green, worker 2 is blue and worker 3 is orange.Worker w 1 is assigned to the set of visits (1,4,6,10,11,9).Worker w 2 is assigned to the set of visits (5,7,8,9,12,13,14,15).Worker w 3 is assigned to the set of visits (5, 2, 3).Note that two workers are assigned to visits 5 and 9, to perform different tasks for the same visit.

Problem constraints
Table 1 lists the objectives and constraints in the WSRP considered here as described in Laesanklang and Landa-Silva (2016).The binary decision variable x w i, j indicates if Assign preferred workers with a specific skill Assign workers to preferred areas a path connects two nodes (visit i and visit j) or not.The assignment x w i, j = 1 means that worker w travels from visit i to visit j, thus w makes both visits as given by constraint (1).
i∈T j∈T In WSRP scenarios, it is possible that some visits are left unassigned as there is not enough workforce or no worker has the required qualifications/skills.In such cases, an integer variable y j is used to indicate the number of unsatisfied assignments for visit j, i.e. visit may require more than one worker (Laesanklang and Landa-Silva 2016;Rasmussen et al. 2012).
If visit j is fully assigned then y j = 0, otherwise y j takes a positive integer value equal to the number of workers required to the visit.Constraint (2) ensures this requirement is met even for visits that are unassigned, r j is the number of workers required for visit j.
The path for each worker w should begin at a start location and terminate at an end location, e.g.their home or a central office.The start location and the end location of worker w are D w and D w , respectively.The condition is enforced by constraints (3) and (4).Workers may leave their start location and enter their end location at most once (although the start and end locations may be different) as expressed by constraints ( 5) and ( 6), respectively.
Workers are required to have suitable skills for every assigned visit.Let q w j be a binary parameter that represent a qualification parameter, where q w j = 1 when a worker w has the skills to take visit j, and q w j = 0 otherwise.Only qualified workers can make the visit as indicated by constraint (7).
Travelling between visit locations must be feasible in terms of travel time.Decision variable a w j takes a positive fractional value that gives the arrival time of worker w to the location of visit j.Note that the maximum arrival time value here is 1440 min, which is equivalent to the 24th hour of the day.Let a w i , a w j be the arrival times of worker w at the locations of visit i and visit j respectively.The arrival time at visit j must consider the time duration δ i spent on performing visit i and the travelling time t i, j between visit i and visit j.This is enforced by constraint (8) where M is a large constant number.
A worker w must arrive at visit j within the given time-window.For visit j, the earliest arrival time is κ L j and the latest arrival time is κ U j .This requirement is enforced by constraint (9).
If a visit j assignment has time conflicts with the visit i assignment, τ w j = 1.A time-conflict occurs when a worker is assigned to visits overlapping in time.
Time availability can be different for each worker according to their individual contracts.The time availability period for each worker is as follows.The shift starttime and shift-end time of the worker w are indicated by α w L and α w U respectively.However, in the scenarios tackled in here, visits can be assigned outside the workers shift but subject to a penalty cost.A binary decision variable ω w j = 1 is introduced to indicate such penalisation.The time availability constraints for worker w are given by expressions (10) and (11).
Another working regulation tackled here is not to exceed the maximum working hours for each worker.Each visit j requires δ j minutes to be completed.The maximum working hours for a worker w is given by h w .Constraint (12) enforces this regulation.
Each worker is associated to a set of geographical regions defined by the service provider.A geographical region contains several visit locations and a visit location may have several visits to be assigned.Ideally, a worker should only be assigned to visits in those geographical regions.However, if necessary, a worker can be asked to travel to locations outside their geographical regions subject to some penalty cost.A binary parameter γ w j = 1 is defined to indicate that visit j is located in the workers regions and γ w j = 0 otherwise.A binary variable ψ w j = 1 to indicate that visit j assigned to worker w is outside the workers regions, and ψ w j = 0 otherwise.Constraint (13) presents the relation between these binary variables for the different possible cases.
Some of the constraints expressed in the above MIP formulation are soft constraints in WSRP scenarios, including constraints (10) and ( 11) (workers may be asked to work outside their shift hours) and constraint (13) (workers may be asked to work outside their geographical regions).Later, preferences calculations are explained as they are used as part of the objective function.

Objective function
The objective function includes the operational cost and the penalty cost.The operational cost includes wages plus journey costs for all workers and is given by the accumulated cost d i, j + p w j , where d i, j is the distance travelled between visit i to visit j and p w j is the cost of assigning worker w to visit j.These costs are set by the service provider in the HHC scenarios used here.
Since feasibility is not guaranteed in a WSRP solution, the penalty costs are tackled as the accumulated penalty for violations of the constraints presented in Table 1.
An assignment of worker w to visit j is made to a path connecting between two nodes and can be written as a tuple x w i, j , y j , a w j , ψ w j , θ w j , τ w j .This assignment is composed of binary decision variables indicating violations occurrences on area availability (ψ w j ), time availability (θ w j ) and conflicting assignments (τ w j ).If a visit j violates time availability, then θ w j = 1.The same applies to area availability violation where ψ w j = 1, when violation occurs.Likewise, τ w j = 1, if the assignment has time conflicts with the other assignment.
The non-satisfaction of preferences is also included in the penalty cost.There are three types of preferences including preferred worker-customer pairing, worker's preferred region and customer's preferred skills.There is a degree of satisfaction for these preferences when assigning a worker w to a visit j and is given by ρ w j which has a value that ranges between 0, 3 , where 0 means no penalty charged and 3 is full penalty.
For each assignment, the satisfaction value for each preference ranges between 0, 1 , from not satisfied to satisfied.There are four-level preferences: low (0.2), neutral (0.5), preferred (0.5), and most preferred (1.0).The satisfaction level is calculated by reverting it to a penalty, by subtracting it from the full satisfaction score, which is 3r j for a visit j.The best solution should have: the least operational cost and the least penalty cost.A weighted sum is proposed to combine the objectives into a single scalar value (Pinheiro et al. 2016;Algethami et al. 2016).The objective function is written as in equation ( 14), where weights λ 1 , . . ., λ 5 are defined to establish priority between objectives (more about the weights used here next).

Weights calculations
The weights associated with each objective are set to values that reflect the difference between the priority levels, as suggested by Rasmussen et al. (2012) and Castillo-Salazar et al. (2016).The main goal is to maintain a balance of priorities for each instance based on each problem specifications.Thus, the weights calculations presented in Table 2 differ from one instance to another.1 Hard constraint violations are not allowed in the WSRP solution.However, time conflicts constraint is more difficult to satisfy.Thus, the highest priority level is given to minimise the conflicting assignments τ w j , where the associated weight λ 5 is set to the highest value.In this work, a solution with conflicting assignments is an infeasible solution.
The second highest priority to be minimised is the unassigned visits y j where λ 4 is set to be very high, but still less than the value of λ 5 .This is due to the service provider requirements for completing as many visits as possible.However, in this study, a chromosome representation ensures all visits to be assigned to workers, hence no unassigned visits violations.This is explained later in the next section.
In practice, the service provider may ask workers to undertake visits that are outside their time availability and/or geographical region.Thus, the next objective priority is λ 3 , given to minimise the soft constraints penalty, i.e., the number of workers with timeavailability violations ψ w j and the number of workers with area-availability violations θ w j .As presented by Table 2, the weight values for λ 3 , λ 4 and λ 5 are relative to the number of the assigned visits, therefore only assigned visits are violated.On the other hand, if the service provider could not fulfil the highest preference level; the fourth priority is given to minimise the preferences penalty through λ 2 .
Finally, the lowest priority is given to minimise the operational cost through λ 1 .This weight value is presented in Table 2, where k is the operational cost.Hence, mileage is calculated for the assigned workers only.
Hard constraint violations are always penalised with a larger weight than soft constraints and therefore the total preference penalty cost.Accordingly, λ 3 + λ 4 + λ 5 > 1 while λ 1 < 1.Under those conditions, the objective is shifted to minimise the constraint violations, by moving the f (S) cost-value closer to the feasible region.

Problem instances
Problem instances from six UK real-world HHC scenarios are used as instances of WSRP2 in this study.There are 7 problem instances in each scenario for a total of 42 instances.
Table 3 shows the main features of each problem instance.Scenario A instances are considered the smallest, while instances in scenario F are the largest.Problem instances in scenario C are of disproportional nature as the number of workers is much larger than the number of visits.
There are some real-world problems that have similar characteristics to HHC, however, they are not tackled in here, as they are beyond the scope of the study.Such as security guard routing and rostering, worker allocation and maintenance personnel scheduling (Castillo-Salazar et al. 2016).
Figure 2 illustrates a group of box plots for each WSRP problem set, that present the overall total of the number of visits assignments |V | and the total number of workers |W |.These plots provide a useful way to visualise the problem size by the vertical distance between the smallest value and the largest value, including any outliers.
As it can be seen from Fig. 2, the gap between the number of visits and the number of workers increases through the problem sets.Distributing assignments of visits over a set of workers is very significant in WSRP solutions, therefore any surplus between the number of visits assignments and the number of workers results in prioritising one over the other.
In some of the problem sets, the number of visits is larger than the number of workers available.Accordingly, any violations of visitations requirements constraints, such as conflicted visits, increases the value of the penalty function.On the other hand, in problem set C, the number of visits is less than the number of workers; hence, the opposite occurs with violations on workforce-related constraints.
The number of visits assignments and the number of available workforce determine the problem size of a WSRP.A small size WSRP in this study has around 50 tasks and 30 workers, as in the case of A and B. Whereas a mix size problem has around 80 visits and 800 workers or more, as in the case of C. A large problem size has more   than 400 visits and around 270 workers, as the case of D and E. Problem set F has an average of 1470 visits and 900 workers, which characterise it as even larger than large instances.Table 4 presents an overview of optimal solutions provided by a mixed integer programming (MIP) solver in Laesanklang and Landa-Silva (2016).The cost-values, i.e. objective function values, is denoted as f (S).The computational-times is denoted as C pt.
Researchers have reported that real-world instances of the WSRP are difficult to solve (Misir et al. 2011;Castillo-Salazar et al. 2016).This was proven further by the MIP, where it was unable to provide solutions for 24 instances denoted as n/a.
The problem size affects the algorithm performance, where small instances are less challenging computationally.Hence, MIP has provided optimal solutions for A, B and some of the C problem set.Nevertheless, the MIP ran out of memory for D, E and F problem sets.Thus, larger problems require more efficient algorithms with problem specific-information.

Adaptive multiple-crossover genetic algorithm (AMCAGA)
This section describes the proposed adaptive aspects of the AMCAGA as an extension of the GA approach in Algethami et al. (2016) and the diversity-based adaptive GA (AGA) presented in Algethami and Landa-Silva (2017).
Algorithm 1 shows the steps of the proposed adaptive multiple crossovers genetic algorithm (AMCAGA).The proposal here is a combination of the adaptive genetic algorithm (AGA), which includes a control mechanism for genetic operator rates P c and P m , and the adaptive multiple-crossover genetic algorithm (AMCGA).The reason for selecting AGA as part of AMCAGA is because it has provided the best results in comparison to other adaptive variations implemented for WSRPs, as stated in Algethami and Landa-Silva (2017).

Algorithm 1 Adaptive Multiple-Crossover Genetic Algorithm (AMCAGA)
Require: A crossover rate P c , a mutation rate P m , set of crossovers χ = x 1 , x 2 . . .x L , WSL for all visits.1: Create a population P of M individuals using WSL 2: repeat 3: Evaluate each individual in P with equation ( 14) 4: for 9M/20 times do 5: Select p 1 , p 2 ← P 6: end for 7: Select , for all pairs of parents 9: Score i ← Performance of x i 10: Update The proposed AMCAGA works as follows.First, an initial population P of M individuals (1-day plans) is created based on an indirect chromosome encoding to ensure solutions feasibility (line 1).At the start of each generation, 9M/20 pairs of parent individuals are selected using binary tournaments (line 4).With some probability P x i , a crossover x i is selected, using the roulette wheel concept (line 5).This crossover, with some probability P c , is applied at a time to each pair of parents to generate two offspring (line 6).The operators' performances are evaluated and recorded as scores(lines 7-8).With some probability P m , each offspring goes through a flat-cost flip (FCF) Fig. 3 An indirect chromosome encoding scheme illustration mutation operator (line 9).At the end of each generation t, the population is updated using an elitism strategy, where the M/10 best individuals from the current population are kept.Along with the 9M/20 offspring individuals generated, the new population of M solutions is formed.If there are no cost value improvements on the best so far solutions for a number of generations, the GA is considered to have stagnated.Thus, P c and P m are modified at every generation, depending on the results obtained in the previous one, to maintain a diverse population and therefore improve the effectiveness of the search (lines 11-22).The indirect chromosome encoding, genetic operators, adaptive parameter updates and multiple crossover mechanism are described next.

Indirect chromosome encoding
According to Osaba et al. (2014), the crossover process is not efficient for the optimization capacity of the technique when it is applied to routing problems using path encoding.Thus, indirect chromosome encoding scheme was proposed in our earlier study (Algethami et al. 2016).The aim is to construct feasible solutions by having suitable workers only.To do so, a worker suitability list (WSL) is created for each visit k, , where k = 1, 2 . . ., |V |.Suitable workers are ranked by a penalty score M k, j , i.e., the lower the score value, the better suited is a worker.The score is calculated for a worker w j by estimating the impact of assigning that worker to that visit, considering incurred operational cost and penalty cost due to preferences and availability restrictions.
Initial population is created randomly by generating a vector of |V | integers between 0 and L k − 1, where L k is the length of WSL for a visit k.When a worker is being considered for a visit, the solution is evaluated for time-conflicts.If it occurs, by the random assignments, the next worker in the WSL for that visit is considered until a suitable worker is found with no time-conflicts arising.Figure 3 illustrates an example of the indirect chromosome encoding for a day plan with seven visits.Each visit has a WSL of four suitable workers, with the best worker for that visit at the top, followed by the next best worker and so on.Below the chromosome, the decoded solution shows the actual worker assigned to each visit.On the right, the encoded solution is shown with an index of the workers as in the WSL for each visit.Time-conflicts are indicated with *.
For example, w 2 is assigned to both v 1 and v 3 while w 3 is assigned to v 2 .No timeconflict arises because v 1 and v 3 do not overlap.However, w 2 is assigned to v 4 and a time-conflict arise as v 3 overlaps with v 4 .Then, the next most suitable worker that does not provoke a time-conflict, in this case w 5 , is assigned to v 4 .The WSL decoder in this indirect chromosome encoding scheme helps to assign suitable workers to visits while avoiding time-conflicts.The penalty scores are shown in Fig. 3 are not used during the generation of initial solutions, they are utilized for the tailored genetic operators described in Algethami et al. (2016).
Note that the representation is designed to include all assigned visits |V |, thus chromosome length varies according to the problem-size.

Genetic operators
The AMCAGA incorporates various genetic operators including some problemspecific ones that utilise heuristics to generate improved offspring (Algethami et al. 2016).All operators used are taken from the study (Algethami et al. 2016).This is because they were proven to perform well under WSRP.They are applied as part of the AMCAGA, with probability P c .For an effective collaboration, a mix of traditional and problem-specific operators are used.Four well-known operators are considered, including one-point crossover (1PX) (Hartmann 1998), two-point crossover (2PX) (Hartmann 1998), uniform crossover (UX) (Hartmann 1998) and half uniform crossover (HX) (Hartmann 1998).Three problem-specific operators are considered, specially designed for the solution representation considered here.Including flat-costs crossover (FCX) (Algethami et al. 2016), partially-matched flat crossover (PMFCX) and flat-cost flip mutation (FCF) (Algethami et al. 2016).
The reason for using a diverse set of operators is to allow dominance of one crossover over the others.The problem focused operators, on the other hand, have the advantage of using heuristics on WSRP instances, where they perform better than traditional crossovers as stated in Algethami et al. (2016).These operators extend the WSL capabilities.The operators used in this study are as follows.
1. 1PX A point is selected at random between 1 and the chromosome length.The genes before this point are copied from one parent and the genes after are copied from the other parent.2. 2PX Two points are selected at random, between one and the chromosome length.
Alternating segments are swapped to get new off-springs.3. UX The number of crossing points are not fixed.Instead, a mixing ratio (50%) is used to choose a uniform random real number u from interval 0, 1 .If u ≤ 0.5; then genes are swapped from parent's individuals, to be copied into offspring.Otherwise, copy the same genes with their current positions from the parent's individuals (Chu and Beasley 1998).
4. HX Similar to UX, a mixing ratio is used.However, exactly half of the nonmatching gens are swapped.Thus, the Hamming distance (number of differing gens between the two parents) is calculated and divided by two. 5. FCX Uses penalty scores that were initially calculated in the WSL at the start of the GA.These values are denoted as 'flat-cost', where each worker w j has a penalty score according to their suitability to work in visit i. FCX goes through each of the V positions in the parent's chromosome.A gene-wise comparison is enforced, for each gene in i th position, with respect to the WSL estimated penalty scores M i, p i .The best suitable worker is given to offspring (o 1 ) and the other worker for offspring (o 2 ). 6. PMFCX This crossover selects a segment within two cutting points.Positions of genes are reversed between these points and the FCX is applied to fill the rest of the offspring.7. FCF Introduces new workers to the chromosome, even if the worker are not suitable for the corresponding visit.A random position i of the chromosome is replaced.Hence, FCF increases the diversity by the random process, even if the best worker for visit i is not selected.p i = (0, 1, 0, 0, 2, 1, 0).If positions 3 is selected at random, gene 0 is replaced by a random number within the list of visit i in WSL, to generate the child chromosome o j = (0, 1, 3, 0, 2, 1, 0).

Adaptive parameter rates
Initial parameter values, P c and P m , were selected after an offline tuning method (Algethami et al. 2016).After the new population is created, if there are no cost value improvements on the best so far solutions for a number of generations, the GA is considered to have stagnated.This means that the search process could be trapped in a local optimum, or that the population could be concentrated in the same region of the solution space.Thus, P c and P m are modified at every generation, depending on the results obtained in the previous one, to maintain a diverse population and therefore improve the effectiveness of the search.Such feedback from the search is used to generate the next population.To avoid early convergence, a diversity-based scale is used to calculate the required change in the adaptive features.Population diversity in here is measured from two aspects.
• The genotype space (Morrison and De Jong 2002), denoted as Diver sit y g .It is the distribution of pairwise differences between individuals in a population.• The phenotype space, denoted as Diver sit y p .It is the population variance, i.e., how far each individual in the population is from the mean value.
The two measurements accommodate different views of the loss of diversity.Combining these methods can overcome those difficulties provided by one measurement used in isolation ( Črepinšek et al. 2013).
Population diversity has an effect on the setting of parameter values (Algethami and Landa-Silva 2017).For example, if instances are relatively small, a large population size M is required as a result of the large degree of similarity between individuals in P. When similar individuals are recombined, inbreeding occurs, and no additional diversity is added through crossover.For a large instance, there is less chance to have similar solutions, thus, a small M is sufficient.
Whenever the best solution found by the technique has not been improved in the last generation, rates updated, according to the calculations presented by Algethami and Landa-Silva (2017).This means that the search process did not evolve correctly and that it is necessary to diversify the population through adaptive parameter rates.
Diver sit y p −Diver sit y g Diver sit y p A resetting process is also used for the adaptive rates P c and P m to forgets the rates history if the rates values are out of range, where rates starts from P c = 45% and P m = 0.10%.This forgetting step means ignoring all the previous feedback process that led up to the inflation of the adapted values.More details of the AGA framework is described in Algethami and Landa-Silva (2017).

Adaptive multiple crossover framework
In relation to the multi-crossover feature, the proposed AMCAGA uses more than one crossover operator, from a set of crossovers χ = x 1 , x 2 . . .x L of a size L, which are alternated during the execution, similar to Osaba et al. (2014).However, the replacement strategy used in this paper differ to the one described in Osaba et al. (2014), by using allocation rules rather than random replacement.
At the beginning, a number of the crossover functions are applied, until the next generation is created.This number is indicated as the memory size σ .One operator is assigned at random and it is randomly replaced by another crossover .All crossovers are applied uniformly, to ensure all crossovers are used, while allowing repetitions.There is at least 1 L chance for each operator to be selected.The operators' performances are evaluated and recorded as scores for a number of iterations during the learning process (also denoted as a cycle of a size υ).For each crossover, the accumulated scores are stored into the crossover reward matrix (CRM), considered as a reward registry for that crossover in the current cycle.Scores are then transformed into application rates P x i , where P x i is the probability of applying a crossover x i .The operators' performances are evaluated and recorded for each cycle, where application rates are dynamically adjusted as the search progress.
Better operators have a higher score value, and therefore a higher probability to be utilised more, and week operators are not left without a chance.For example, given a set χ of six crossovers operators, with application rates as follows: x 1 = 23%, x 2 = 4%, x 3 = 11%, x 4 = 34%, x 5 = 17% and x 6 = 12%.Crossover x 4 has the highest application rate at this cycle, henceforth x 4 has more chance to be implemented.At some point of the search, when x 4 performance changes for the worse, its application rate decreases.Application rates of other crossovers also vary accordingly, causing a shift in the search.
Each operator is implemented independently; however, all operators share the same cycle until it resets.Hence, AMC is considered as a mix between natural systems (recombination operators) and synthetic systems (autonomous agents).The natural system aspect is reflected by the collective performances of the crossovers, onto a population (Eiben et al. 2007).The synthetic systems aspect is reflected by sharing a cycle as part of the learning process that directs their application (Panait and Luke 2005).Details for crossovers performance measurements and CRM construction is explained next.

Performance measurements
Scores are considered as performance indicators that enforce the rewarding mechanism, where ∀x i ∈ χ , a score is given at each iteration, where i = 1, 2 . . ., L. To calculate the scores, performance measurements are used to model the behaviour of different operators in the current stage of the search.
Different types of measurements are used in the literature, such as fitness-based measurement (Whitacre et al. 2006;Li et al. 2014), distance-based measurement (Consoli and Yao 2014) and combined measurement and operator execution-time measurement (Maturana and Saubion 2008).In this study, the later can be relative to the problem size for WSRPs.Hence, the time is not absolute in real-world settings, and therefore this measurement is excluded.
Parameter settings are affected by the population diversity (Algethami and Landa-Silva 2017).Therefore, distance-based measurement is considered here along with the fitness-based measurement.Performance measurements used in this study are as follows.
1. Fitness-Based, This method is selected to maximise the cumulative improvement, as a historical fitness record of the cost value of an offspring (o) and its parent ( p), where f (o) < f ( p).When there is an improvement on the overall cost value, the score increases by one.This value ensures the convergence by selecting crossovers with better offspring quality, that eventually improves the overall performance of a GA. 2. Distance-Based, This method is selected to ensure a level of diversity among generated solutions, by calculating the distribution of the pair-wise differences between an offspring and its parent.If the percentage is 60% or more, the score increases with respect to the hamming distance H between two individuals.This value is selected to prevent inbreeding among generations.The value of H is calculated as the number of assignments in which the corresponding genes are different.Figure 4 illustrates an example of the dissimilarities percentage calculations (presented in white cells) between a parent and an offspring.3. Hybrid Approach, This method is selected to harness the power of the two above measurements.When there is an improvement in the fitness or there exist differences between parents and offspring, the score increases.This means that with every crossover evaluation, the score value increases gradually.If only one measurement has proved an operator efficiency, the score value increases by one.On the other hand, if both performance measurements have proved an operator efficiency simultaneously, the score value increases by two.
Each crossover is evaluated in isolation, after each application, with a performance measurement.The purpose of such separate evaluation is to obtain their relative score, under the evaluation criterion.In the event of one crossover is selected repeatedly, the score value increases gradually, based on the feedback from the online learning mechanism in AMC.Details for application rates calculations is described next.

Application rates
At each iteration (cycle), one of three above mentioned performance measurements (fitness, distance or hybrid) evaluate all crossovers, and then gives a score value ∀x i ∈ χ , where i = 1, 2 . . ., L. For the next iteration, scores are updated, until all cycles are completed.The accumulated score value over different cycles is then transformed to an application rate P x i = Score i Score Sum , where Score Sum = L i=1 Score i .When recording the accumulated scores, a 2 × L matrix is used, noted as crossover rewards matrix (CRM), i.e., the score value in (1, i) is added to the value in (2, i).The 2 × L CRM is used to allow faster computation, rather than using a L size matrix.
As a result of scores updates, a gradual change in the P x i value exists in CRM, with no rapid increase or decrease in the overall application rates.This process is repeated, for each cycle.Note that the application rates are always modified before the next cycle (iteration).
Figure 5 shows an example of the CRM construction process, for three crossovers x 1 x 2 , and x 3 , with a cycle size υ = 3.For every cycle, two steps are executed.First, scores are retrieved for each crossover.In this case, the first cycle has a Score 1 = 3, Score 2 = 4 and a Score 3 = 5.Next, the scores are shifted to (2, i).Once a score is shifted to position (2, i), a value of 1 is stored at position (1, i), to ensure the application of all crossover functions, regardless of a crossover performance.If a crossover application has resulted in no improvement or even worse solution, it is still applied in the AMC.The accumulated scores are then calculated by adding up the scores at the position (1, i) to the scores at the position (2, i).The active calculation of the application rates maintains a balance between the most effective operators, which provide good results, and week crossovers, which provide poor results.

Experimental study and results
The implementation has focused on applying the AMCGA on the GA presented on Algethami and Landa-Silva (2017) to capture the advantages of this method on the GA performance.For these set of experiments, the main elements used were the indirect chromosome representation and the FCF mutation.In this paper, six operators were used to allow dominance of one crossover over another, including 1PX, 2PX, UX, HX, FCX and PMFCX.Since parameter settings and running times stated in Algethami et al. (2016) were successful in providing good results, each run in these experiments was executed with the same set of settings as shown in Table 5.Note that the crossover rate P c and the mutation rate P m values stated in Table 5 are used as initial settings for the AMCAGAs.However, these values are fixed when used with AMCGAs.
Similar to genetic operators rates, each WSRP problem set has different memory size σ .This value indicates the number of crossovers required to obtain enough information for the learning process, where P c × σ = const.Thus, the crossover rate P c affect the value of σ , to record the crossover performances on at least M/2 of the population.
For example, in smaller instances, P c = 50%, therefore, less crossovers are applied.Hence, larger σ is required to obtain sufficient information about the crossovers performances.In larger instances, P c = 100%, this means that more crossovers are applied, and a large memory size is not necessary.Different cycle settings were used here, where υ = {5, 10, 15, 25, 40, 50}.The reason for selecting these values was due to the preliminary experimentation.The testing revealed that a cycle of a size less than 5, had a slow learning process.Accordingly, all operators were applied equally.A cycle of a size larger than 50 had a low accuracy of the learning factor, in which crucial information was lost.For example, a crossover operator might be the best performer at the start of a cycle, yet it can get worse in the same cycle.
Results were grouped by the performance measurement method applied.Each measurement is noted as follows: fitness improvements (AMC f G A), dissimilarities between individuals (AMC d G A) and the hybrid approach (AMC h G A).The best performing method was selected to be included in AMCAGA and combined with adaptive parameter control method AG A, presented in Algethami and Landa-Silva (2017).
Each algorithm was executed 8 times (runs).This means that for each of the 42 problems instances, there were 8 × 6 (cycles) × 3 (methods) = 144 runs, seeded with the same initial population.

Performance of AMCGAs with different cycle sizes
This set of experiments aims to investigate the effectiveness of the feedback mechanism in the AMCGAs, when using different cycle sizes.To do so, different cycle sizes were examined with respect to the method applied.
Figure 6 shows the total number of the best solutions for all problem instances.Each bar at the x-axis represents a rewarding method applied with a cycle of size υ.Methods shown in the plot are: AMC f G A (grey bars with striped lines), AMC d G A (black solid bars) and AMC h G A (blue bars with dots).The higher the bar the better, i.e. the more "best" values.
It is clear that all methods provided relatively similar values with a slight increase on one of the methods over the others.Still, AMC f G A obtained the highest number of best solutions.For AMC f G A and AMC d G A, the highest number of best solutions was the when υ = 10.For AMC h G A, on the other hand, the highest number of best solutions was when υ = 25.This result is related to scores given to an operator throughout the search, in which the accumulated score increases or decreases gradually.An operator performance at the start of the cycle might changes drastically by the end of the cycle, especially when a large υ value is used.Nevertheless, υ > 25 was proven to be less effective for all methods.Interestingly, the highest number of best solutions was obtained while using AMC f G A, with υ = 10.Hence, there is no need for larger υ value to keep long historical records.
With one performance measurement, i.e., AMC f G A and AMC d G A, small changes in the score values occur.This provided enough time for updating the application rate, for each crossover, and therefore these crossovers were used in the current population.On the other hand, a large υ was required when the hybrid performance measurement was applied.This is because of a large and discrete change in the score values that occurs as the outcome of combining both performance measurements.Nevertheless, if the υ was small; there was not enough time to reward all operators.Thus, one operator will dominate the algorithm, resulting in an uneven distribution of the application rates.Henceforth, a larger cycle-size was required.
In summary, the larger the υ, the vast increase in the score value.This process has more effect on the search than having scored in small cycles, in which are ignored earlier in the search.Thus, υ value can affect a crossover dominance in the AMCGA method over another.
Figure 7 illustrates the overall average computational times in seconds, on problem sets A to F, for all methods.Each sub-figure corresponds to a cycle size υ and presents the average computational time used by AMC f G A (grey bars with striped lines), AMC d G A (black solid bars) and AMC h G A (blue bars with dots).The lower the bar the better, i.e. the less computational time.
On average, it was apparent from the plots that different methods provide solutions in approximately similar computational times.
For AMC f G A, it was more computationally expensive to make large adjustments in the score values, especially when smaller υ was used.On the contrary, more time existed when using larger υ values, providing a shift in the score values when AMC f G A was used.Thus, the largest computational times for AMC f G A, among all methods, was recorded when υ = 40 or 50.
For AMC d G A, more time was spent on calculating the scores, that resulted in more computational time, especially when υ was equal 5 to 25.This was the outcome of the diversity-based calculations, computed gene by gene, and executed V times for each individual to calculate the scores accurately.
Fitness-based calculations were faster than the diversity-based calculations.However, when υ was equal 40 to 50, the sensitivity of the distance-based has provided good solutions in less computational time.This is because that distance-based measurement method was more sensitive to the change of performances than fitness-based measurement.Note similar execution time to AMC d G A with AMC h G A was recorded, with υ = 15.
Contradictory to AMC d G A, AMC h G A was not affected by the convergence speed of the algorithm.The combinations between the diversity-based measurement and the fitness-based measurement in AMC h G A have obtained better computational times than the separate methods.This is because AMC h G A have averaged the performances of both measurements by allowing the use of the calculations within the time limitations.The reason for this was the large increase in the score values due to both calculations.As a result, the scores were computed in less time while using AMC h G A.
The next set of experiments is designed to use one cycle-size per method, selected according to the highest total number of best solutions (see Fig. 6).Thus, υ = 10 was chosen for AMC f G A as well as for AMC d G A, and υ = 25 was selected for AMC h G A.

Overall comparison between AMCGAs methods
The second set of experiments was designed to identify if there is a significant difference between the proposed AMCGAs variations or not.The work by García et al. (2008) has recommended the use of Friedman analysis as a non-parametric statistical test to establish statistical significance in EAs.Thus, a two-way Friedman analysis was used to measure the significant difference between groups of data when the dependent variable being measured was ordinal.
In this study, an IBM SPSS 22 two-way analysis was used to compare the variances of seven related-samples, with a significant level of α = 0.05 and 95% as a confidence interval.Table 6 provides the results generated by the Friedman analysis on the 42 problem instances including the mean value, the standard deviation, the minimum cost-value, the maximum cost-value and the mean rank.The results presented in the mean rank column show the methods order based on the statistical analysis, where a low rank indicates the best method while a high rank indicates the worse method ranked overall.All problem instances were used to set the sample size of one method as large as possible, this increases the probability of accepting or rejecting the null hypothesis.Three additional values were calculated and used to measure the performance of each algorithm Dev. is the average percentage deviation from the best-known value (best solution of all the algorithms applied), Best is the fraction of instances in a set for which an algorithm matches the bestknown value (best solution of all the algorithms applied).Additionally, the Score is the fraction of the instances for which a competing algorithm 'wins', i.e. produces better solutions than the configuration being scored.This score is calculated as , where p is the number of methods compared, q is the number of problem instances, and r is the number of instances in which the p − 1 competing configurations find a better result.Hence, the best score value is 1, when r = 0, and the worst score value is 0, when r = q × ( p − 1).The best results are highlighted in bold.
We applied the Friedmann non-parametric statistical to the data in Table 6 and obtained a p value of 0.02153 < 0.05, degrees of freedom = 2 and χ 2 = 7.677.This indicates the existence of significant performance differences among the three methods.
To examine where the differences actually occurs, an additional analysis was implemented.A post-hoc statistical analysis is needed in all cases.Holm's test was chosen to detect the significance difference among all variations.The Holm procedure is an example of a step-down procedure.Step-up procedures start testing hypothesis H m and step up through the sequence while retaining the hypotheses.The procedure stops at the first rejection (for example H i ), and H 1 , . . ., H i are all rejected.
In this case, Holm's method obtains the p-values higher than the significance level, that is to be interpreted in the sense that we do not have enough evidence to reject the null hypothesis.
However, the descriptive statistics and the measures explained above revealed that AMC d G A method had the lowest mean value, standard deviation value, mean rank value and Dev% value.This method had also the highest fraction of the number of best solutions and the score values.This finding suggests that different methods applied to different datasets can obtain different results and henceforth each problem set benefited from each method accordingly.
AMC d G A obtained the best values in the above table.Thus, it can be argued that all variations were suitable for WSRP while AMC d G A was slightly better than all the proposed multiple crossovers algorithms.Using the distance-based measurement has provided more accurate scores adjustments in comparison to other methods.

Crossovers dominance based on application rates distribution
This section explores the change in the crossovers' application rates (P x i ) to investigate the effectiveness of one crossover over the other, in each AMCGA.Thus, this set of experiments aims to identify the most used operator in each method applied, with respect to the performance measurement.Each bar at the x axis in Fig. 8 illustrates the average application rate P x i values, for all problem sets.All three methods were applied, with the set of crossovers, i.e. 1PX, 2PX, UX, HX, FCX and PMFCX.The AMC f G A is plotted in grey bars with striped lines, the AMC d G A is plotted as black solid bars and the AMC h G A is plotted as blue bars with dots.
Among all crossovers applied with AMC f G A, the values for the average P x i were relatively similar with a slight increase on FCX.This result further proves the observations stated in Algethami et al. (2016), that identified the FCX as one of the best crossovers for WSRP problem sets.Thus, this operator has been utilised more when a fitness-based performance measurement was utilised.
Nevertheless, the highest average application rate P x i among all crossovers was obtained by PMFCX, while using AMC d G A. Hence, PMFCX was utilised more often when distance-based performance measurement was applied.
Faster calculations were required when using AMC d G A and AMC h G A. Thus, FCX was used less under these methods.On the other hand, the PMFCX crossover was utilised more, due to its fast calculations as a result of mixing heuristic approach with the traditional PMX crossover method.
Traditional crossovers have provided poor quality solutions that were comparatively similar to their parents.Still, using them alongside the problem-specific method can ensure various performances in the GA.
Figure 9 illustrates a group of box plots comparisons to show overall patterns of response of change in application rates (P x i ), for each crossover operator.Crossovers distribution corresponds to the current environment state if there is change in the P x i values, i.e the bigger the box plot the more diverse is the P x i values to the correspondent crossover.
Each sub-figure corresponds to the method applied, AMC f G A, AMC d G A and AMC h G A. Each box-plot illustrates the average application rates (P x i ) range, for each crossover and are shown in 1PX (blue), 2PX (red), UX (black), HX (green), FCX (grey) a PMFCX (brown).
As it can be seen from the plots, PMFCX box plot was lower than all other plots in AMC f G A, i.e. less change in P x i values.The opposite occurs for AMC d G A and AMC h G A, where P x i values were more diverse.This indicates that PMFCX was updated more frequently.This result further supports the previous finding discussed in this section.

Effect of using a rewarding process on AMCGA
The third set of experiments aims to compare AMCGA methods with existing GA that was tailored for WSRP.The detailed results shown in Table 7 compares the indirect GA described in Algethami et al. (2016) (noted as GA) with the AMCGA variations (AMC f G A, AMC d G A and AMC h G A). Crossovers used with indirect GA were specified for each problem set as follows.for A and B GCX, for C PMGreedyX, for D, E and F FCX.The FCF mutation was also used.
In addition to the indirect GA, the AMCGAs were compared against a variation of AMCGA that used a uniform choice of operators, noted as AMC r G A. The gal for this comparison is to investigate the performance of AMCGAs even when excluding the reinforcement learning process.For each problem instance, the table shows the solution quality, noted as f (S), and the computational time in seconds, noted as Cpt, in which the best solution was found.
Best values are highlighted in bold.If more than one method achieved the same result, among cost-best equals, the time-best is highlighted in bold.
As it can be seen from Table 7, the proposed methods AMC f G A, AMC d G A and AMC h G A outperform the GA and AMC r G A in terms of computational time, in particular AMC f G A and AMC d G A. This indicates that the rewarding process in the adaptive methods were more efficient time-wise.The percentage of best cost-values overall solutions are as follows.The GA 14.29%, AMC r G A 30.95%, AMC f G A   Using a variety of operators, in which they have different performances, have provided good quality solutions for each problem set, by AMCGAs methods, especially AMC f G A and AMC d G A. Improvements have occurred due to the combined work of the operators, which allowed an extensive search with a larger variance than using one crossover.However, using two performance measurements, as the case in AMC h G A was proven to be inefficient in comparison to the other AMCGAs.This is due to the large jumps in the score values in comparison to using one performance measurement, i.e.AMC f G A and AMC d G A, where small movements into various directions in the solution space has resulted in finding undiscovered regions easier.
By contrast, in the indirect GA, the search was completed by one crossover involvement in a large proportion of the solution space.Using one crossover throughout the search has less ability to extend the search in those regions, which were most promising.This was the reason for the GA providing worse results than the AMCGAs.
On the other hand, random crossover selection has utilised different crossovers at an arbitrary level, without prior knowledge of the current search space.Still, the AMC r G A has obtained the best results on E 5 , E 7 , F 3 and F 5 .This was because of the lack of performance measurements calculations.Hence, results were computed faster under those instances.
In regard to computational time, AMCGAs proved to improve the efficiency of the algorithm by computing the results in less time.For GA, AMC r G A , AMC f G A, AMC d G A and AMC h G A, the average C pt in seconds were 6069.04, 5756.7, 5797.7, 5873.8 and 5825.6 respectively.The reason for the rapid decrease in the computational times was as follows.
In the traditional GA, one crossover was applied at each iteration.As a result, the algorithm required more time to improve the solution.On the other hand, the AMCGAs enforces the performances of different crossovers onto one population, which was made in a minimum time.Therefore, the proposed AMCGAs obtain solutions in fewer iterations than the indirect GA.In AMC r G A , the average C pt was the lowest among all compared methods, this was due to the exclusion of the feedback process.The average C pt for AMC f G A was less than the average C pt for AMC d G A and AMC h G A. This result was the outcome of computing the dissimilarities between the parent and the offspring, which resulted in a massive amount of calculations.Despite the fact that in AMC h G A the two performance measurements were applied, it obtained the results in less C pt-time than AMC d G A. This outcome was due to the use of the rewarding process that resulted in better and faster convergence.

Using an AMC d GA as a component of AMCAGA
A major advantage of AMCGAs was that they have a considerable influence on improving the efficiency of the baseline GA, especially when the distance-based performance measurement was applied.Thus, the AMC d G A method was selected to be integrated with the adaptive operator rates control method AMCAGA (noted as AMC d AG A), that perform as a full adaptive GA.The following experiments aimed to evaluate the validity of AMC d AG A by providing further insights into its performance.

Effect of using different cycle sizes in AMCAGA
In this section, the aim was to investigate the effect of the using different cycle sizes in AMC d AG A on solutions quality and computational times.
Figure 10 shows the total number of the best solutions generated by different cycle sizes.Each bar at the x-axis represents an AMC d AG A method applied with a cycle size υ.The higher the bar the better, i.e. the more "best" values.Using different cycle sizes has resulted in relatively similar performance, with a high number of best solutions in more than one cycle size.Thus, it can be claimed that using adaptive operator rates (P c and P m ) makes the performance of the AMCGA more stable than using the AMC d G A separately, as seen on Figs. 6 and 10.
The number of best solutions was the highest when υ = 15, 40 and 50 with a value of 17 best solutions.Followed by υ = 10 with 16 solutions and υ = 25 with 15 solutions.The values for the number of best solutions when υ = 5 was the lowest.The time was insufficient to score all operators.Thus, a large cycle size, i.e. 50, provide better learning outcomes with the time to retrieve the information needed in order to improve the results.
Another observation was made when recording the computational times for AMC d AG A while using different cycle sizes.Each sub-figure in On average, the lowest computational times for problem set A was when υ = 25, however, all cycle sizes have obtained the same result in less than 50 s which was relatively acceptable.Followed by the υ = 50 that also obtained the lowest average computational times for problem sets C, D and F. The results, as shown in Figs.11 and 10, indicate that combining the adaptive operator rates (Pc and Pm) approach with AMC d G A of a cycle of size 50 provides the best results in less time, especially for difficult problem sets.The larger the cycle size, the more scores were calculated for an operator effectiveness in addition to updating the operator rates.As a result, the population evolved over time into better, fitter solutions.The interchange between the two adaptive aspects was vastly exploited in AMC d AG A. Note that the lowest average computational times for problem sets B and E were when υ = 10 with a difference between 100 and 200 s to the other sizes.With reference to Fig. 10, this cycle size also obtained the second highest number of best solutions.Therefore, the results for a υ = 10 and υ = 50 was investigated further in the next section.

Overall results of AMCAGA versus AMCGA and AGA
The aim of this experiments is to compare the performance of AMC d AG A against AMC d G A and AG A to understand the effect of combining the adaptive parameter aspects into the AMCs.
Results for AMC d AG A are shown in Table 8 with the diversity-based adaptive operators rate control GA (noted as AG A) and the AMCGA variations that utilised diversity-based measurement (noted as AMC d G A).
The AMC d AG A method used in this comparison was when υ = 10 (noted as AMC d AG A 10 ) and a υ = 50 (noted as AMC d AG A 50 ).These methods are selected based on the observations taken form Figs. 6 and 7.
Note that there are two types of adaptability in Table 8.The first adaptive method used in AGA and AMC d AG A is the GA parameter control, P c and P m , The second adaptive method used in AMC d G A and AMC d AG A is multiple crossover adaptability.
The baseline GA was excluded in this comparison because it was proven that AMC d G A and AG A have provided better results as discussed earlier.Another advantage of this method was the rapid decrease in computational time that was previously reduced by using the adaptive methods separately in comparison to the baseline GA.On average, the computational times in seconds were 5653.11, 5873.79, 5069.67 and 4962.10 for AG A, AMC d G A, AMC d AG A 10 and AMC d AG A 50 respectively.Thus, using the combined method was more cost-effective than using one method individually for WSRP.

Results of AMCAGA versus WSRPs solution methods
This section provides a comparison of the best-performing methods proposed in this study, i.e., AMC d AG A 10 and AMC d AG A 50 , against three existing WSRP applications: MIP solver (Laesanklang and Landa-Silva 2016), MIP with decomposition (Laesanklang et al. 2015) and VNS algorithm (Pinheiro et al. 2016).Table 9 shows the solution quality, noted as f (S), and the computational time in seconds, noted as Cpt, in which the best solution was found.The best values are highlighted in bold.
For smaller problem sets, the proposed methods were quite competitive, matching the best-known results for many of those instances.The VNS seems to provide better overall results with 42.86% of all best solutions.However, the AMC d AG A 10 and AMC d AG A 50 outperformed the MIP with decomposition method with 21.43% and 42.86% of all best solutions respectively.
The average computational times were calculated for methods that provided solutions for all instances.i.e., MIP with decomposition AMC d AG A 1 0 and AMC d AG A 5 0. The recorded times were 4964.26,5069.67 and 4962.10 respectively.
The AMC d AG A 10 has the largest computational time among all methods.On the contrary, the AMC d AG A 50 has the lowest computation time among all methods.These findings provide implications that the use of a large cycle size provided better GA performance, especially cost-wise.A noticeable effect of using adaptive aspects was the reduction of the computation time, which was also enhanced when using AMCAGA.On the other hand, the MIP with decomposition method has obtained a low average computation time.However, the results obtained were poor in comparison to the other methods.
So far, the proposed methods were the only algorithms that provided results for all instances.Hence, adaptive GAs were able to solve real-world and highly constrained optimisation problems.Interestingly, the best cost values were obtained by diversity-based methods.This further proves the significance of maintaining a diverse population in enhancing the GA performance when tackling WSRP.Even though VNS had obtained better results than the proposed methods in this study, when results were available, this study aided in better understanding of the usability of GAs under the combined settings.This was the first time that a study on the design components of a GA, that has been used to solve WSRP.

Conclusion
Using synergies among the operators an provide better results than using one operator at a certain stage of the search (Li et al. 2014).This concept is used in this study by proposing an adaptive multiple crossover framework (AMCGA) to tackle the WSRP scenarios to improve the GA efficiency.
In this investigation, six different crossover operators, well-known and cost-based, are selected as part of the AMCGA method.The crossovers are rewarded with scores during the learning process, which affected their application rates during the run, and are later adjusted with values that reflect the search environment.Rewards are given to a crossover according on its effectiveness at the current stage.To evaluate a crossover, three performance measurements are utilised to provide a score for each operator.The crossover ability to provide an offspring with fitness improvements, to generate an offspring that is different than their parents or to provide fitness improvements and more diverse solutions together.Based on the performance measurements, three variations of the AMCGA are implemented ( AMC f G A, AMC d G A, AMC h G A) with six different cycle sizes on a total of 42 different instances of WSRP problems.The cycle-size, which provided the highest number of the best solutions, for one method is selected in a more in-depth analysis in this study.
Further analysis has identified that AMC d G A has the best score values, the highest number of best solutions with the highest computational time among all methods.This is due to diversity-based measurement requires more computation time than fitness-based measurement.However, it can provide better results.The combination of the diversity and the fitness measurements have obtained less computational times than using separate methods.This is due to the scoring process used.If only one measurement is used; a small increase in the score occurs.On the other hand, if both performance measurements are used a large increase in the score occurs.Thus, scores in AMC h G A are higher than AMC f G A and AMC d G A, henceforth, less time is required.
Another key finding is that the highest average application rate among all methods is obtained by FCX and PMFCX crossovers.Results showed that PMFCX crossover is the most utilised crossover in all methods.
Results are compared with MIP with decomposition (Laesanklang et al. 2015), VNS algorithm (Pinheiro et al. 2016), indirect GA (non-adaptive) (Algethami et al. 2016), randomly uniform variation of the AMCG A (no-learning), noted as AMC r G A and adaptive parameter control GA (AGA) (Algethami and Landa-Silva 2017).In general, the AMCGAs have outperformed the stated methods, especially when diversity measurement is applied.Using a learning scheme has effectively improved the GA performance.
The key strength of this study is the combinations of adaptive GA settings (P c and P m rates) in addition to the adaptive multiple crossover method (AMCGA).This com-bination has resulted in better quality solutions, in less computational times.Hence, using a full adaptive algorithm is the best variation of a GA to be used for the combined scheduling and routing problem settings.Further work needs to be carried out to validate the features of the set of operators used in this study, such as the number of operators and the type of operators.

Fig. 1 A
Fig. 1 A WSRP example, with 3 workers and 15 visit assignments (Colour figure online)

FFig. 2
Fig. 2 Box plots comparisons of WSRP scenarios dimensions, shown as the number of visit assignments V and the number of workers W (Colour figure online)

Fig. 4
Fig.4An illustration of the dissimilarities percentage (presented in white cells) between two individuals (parent, offspring)

Fig. 5
Fig. 5 An example illustration of the application rates calculations Fig. 6 Total number of best solutions produced by different rewarding mechanisms (AMC f G A, AMC d G A and AMC h G A) (Colour figure online) Fig. 7 Average computation-time (in seconds) produced by AMC variations with different cycle sizes υ (Colour figure online)

Fig. 8
Fig. 8 Average application rates values used under different rewarding mechanisms ( AMC f G A, AMC d G A, and AMC h G A) (Colour figure online)

Fig. 10
Fig. 10Total number of best solutions while using AMC d AG A Fig. 11 corresponds to a problem sets from A to F, and each bar illustrates the overall average of the computational time in seconds used by a cycle-size.The bars colour and pattern indicate each cycle as follows: black solid bars when cycle size = 5, grey bars with stripped lines when υ = 10, blue bars with dots when υ = 15, red bars with right inclined lines when υ = 25, green bars with left inclined lines when υ = 40 and yellow bars with a grid when υ = 50.The lower the bar the better, i.e. the less computational time.

Fig. 11
Fig. 11 Average computational-time (in seconds) produced by AMC d AG A used with different cycle sizes υ (Colour figure online)

Table 1
Objectives and constraints in the WSRP

Table 2
Weights values for WSRP

Table 3
Main features of the 42 home health care problem instances

Table 5
Parameter settings for AMCGAs and AMCAGA

Table 6
Non-parametric Friedman's Statistical results combined with performances metrics

Table 7
Results of the cost values f (S) and computational time Cpt (in seconds) produced by the proposed methods GA, AMC r

Table 7 continued
AMC d G A 47.48% and AMC h G A 21.43%.Closer inspection of the results is explained next.

Table 8
Results comparison between cost values f Bold text refers to the best result From the table, it can be seen that AMC d AG A has outperformed AG A and AMC d G A, especially AMC d AG A 50 .Closer inspection of the results shows the percentage of best cost-values overall solutions are as follows.The AGA 21.43%, AMC d G A 35.71%, AMC d AG A 10 42.86% and AMC d AG A 50 59.52%.Hence, AMC d AG A 50 has provided the best results on 59.52% of all solutions.Followed by AMC d AG A 10 that provided the best results on 42.86% of all solutions.These results indicate the power of combining two adaptive elements that work together in order to improve the GA offspring productivity.

Table 9
Cost-values f (S) and computational-time Cpt (in seconds) for WSRP instances, produced by all WSRP solution methods on 42 instances