Crossover

Crossover is the process of mating between the two selected individuals, the process represents how genes are transferred to the offspring. There are several methods of crossover:

Binary

• Single-point Crossover

Split the parent based on a randomly generated single point. The first child would have gene from the first parent for the position 1 up to the point, while the rest is filled with genes from the second parent. The second child is opposite, position 1 up to the point is filled with genes from the second parent.

• Uniform Crossover

Every gene is exchanged between the a pair of randomly selected chromosomes with a certain probability, called as the swapping probability. Usually, the swapping probability value is taken to be 0.5.

Real-Valued

• Whole Arithmetic Crossover

Takes weighted sum of the two parental genes for each position. The weight is applied for all position. Let x = the value of the first parent and y = the value of the second parent for position/gene 1:

\[Child \ 1 = weight \ * x + (1-weight) \ * y\] \[Child \ 2 = weight \ * y + (2-weight) \ * x\]

If the weight = 0.5 two offspring are identical

• Local Arithmetic Crossover

The local crossover is similar to the whole artihmetic crossover, but uses a random alpha for each position.

Permutation

• One-point crossover

One-point crossover use only one random point as the marker where crossover should happen.

• Position-Based crossover

A number of positions from the first parent is randomly selected and fill the remaining position from the second parent.

Mutation

A mutation means change in one or several genes inside the chromosome. Just like in the natural world, mutation rarely occur, thus they have a little probability to happen, usually set at 0.01. There are some example of mutations:

Binary

• Random Mutation

Genes are randomly flipped from 1 to 0 or vice versa with uniform probability for each gene, meaning that each gene has the same chance to mutate with mutation rate or mutation probability set by the user.

Permutation

• Adjacent Exchange Mutation

Two adjacent genes are swapped randomly.

• Inverse Mutation

Inverse mutation is done by inversing the sequence between two random points.

Define Encoding and Fitness Function

Since the problem is about whether bring certain item or not, then the solution should consists a logical value of each item. If the value is 1 or TRUE, then I should bring those item, and if value is 0 or FALSE, then I should not bring the item. By this logic, the encoding for our genetic algorithm should be binary encoding because it consists only 2 values (1 or 0).

The fitness function then, should be to maximize the total point.

\[Maximize \ Total\ Point = \Sigma_{i=1}^{k} \ P_k \ * \ chosen_k \]

• \(Total\ Point\): The total point generated
• \(k\): Number of item
• \(P_k\): Point of item k
• \(chosen_k\): logical value of each item (1 or 0)

Constrained to:

\[Total\ Weight <= 15\]

Since we have constraint, we need to incorporate the constraint into the fitness function. This is done to prevent an unfeasible solution (those that violate the constraint) to be chosen or selected as optimal value

\[Maximize\ Total\ Point = \Sigma_{i=1}^{k} \ P_k \ chosen_k - 2 \ (total\ weight - 15)^2\]

Initial Population

First, GA will randomly generate a population of solution. Here, we have a population of 4 solutions/indviduals. The binary value on each individuals are generated randomly.

Fitness Assignment

Now we calculate the fitness value (total point) of each individuals based on the fitness function above. We can see that solution that violate the constraint (have total weight > 15) will have worse fitness function.

Selection

Now we will choose which individuals are gonna be crossed in order to get a new solution. We will use linear rank selection to select the parent solutions. We will choose 2 solutions since crossover require 2 parent solutions.

Linear ranking selection work as follows:

• Give rank to each individuals. Bigger rank is assigned for bigger fitness function. Individual 2 with fitness function of 19 is the highest, so it will be given rank 4.
• Calculate the probability of each individuals based on their rank. The sum of the rank is 10, so the probability of each solution is \(Rank_i/10\).
• Chose solution randomly based on the probability. Solution with higher fitness function will naturally has higher chance to be chosen.

Let’s say for example, the first and second solution are chosen as the parent

Crossover

Crossover means exchanging genes between two parents two create new solution called child solution. We will use the single-point crossover here. A single point is chosen randomly.

Let’s say the point is chosen to be between position 2 and 3. For the first child, it will have the first and second gene from the first parent, while the rest will be filled from the second parent. For the second child, it is the opposite, it will have the first and second gene from the second parent, while the rest will be filled from the first parent.

Mutation

Mutation is the change of value in the gene. Here we use random mutation. It will randomly swap the value of gene. If the gene in position 3 have value of 1, it would be switched into 0 and vice versa. Mutation is rarely occur, its probability can be set manually, but generally people use mutation probability of 0.1.

The selection, crossover, and mutation process will be repeated by choosing different parent until the number of children are the same with the population size of the parent. If the parent have population size of 4, then we would have a children population of 4.

Check Stopping Criteria

After mutation is done, we would check if it is time to stop the algorithm. Have we reached iteration 100? Do in the last 10 iterations the optimal fitness value didn’t change? If not, then the iteration would continue. After crossover and mutation, now we have parent population of 4 and children population of 4, so in total we now have 8 solutions. For the next iteration, we will only choose 4 individuals (in accordance with the initial population size) from the parent and the children. Generally, the individual with high fitness value will be chosen, so the next population will be better than the last one.

Knapsack Problem

Suppose we have several items to deliver into the distribution center. However, our truck can only load up to total of 10 tons or 10000 kg, so we need to choose which items to deliver and maximize the weight loaded. You can directly solve this problem using the Knapsack Method, but you can also use Genetic Algorithm. Here we will see how GA deal with a constrained problem (weight should not exceed 10 tons or 10000 kg).

df_item <- data.frame(item = c("Tires", "Bumper", "Engine", "Chasis", "Seat"),
                      freq = c(80,50,70,50,70),
                      weight = c(7, 17, 158, 100, 30))

rmarkdown::paged_table(df_item)

Each items have quantities and their respective weight.

Let’s create a new dataset with one observation correspond to only 1 item.

df_item_long <- df_item[rep(rownames(df_item), df_item$freq), c("item","weight")]
  
rmarkdown::paged_table(head(df_item_long))

Let’s set the objective function.

\[Total\ Weight = \Sigma_{i=1}^{k} \ W_k \ * \ n_k \]

• \(Total\ Weight\): The total weight (kg)
• \(k\): Number of unique item
• \(W_k\): Weight of item k
• \(n_k\): Number of item k that is chosen

Constrained to:

\[Total\ Weight <= 10000\]

We will translate the objective function into R function.

weightlimit <- 10000

evalFunc <- function(x) {
  # Subset only selected item based on the solution
  df <- df_item_long[x == 1, ]
  
  # calculate the total weight
  total_weight <- sum(df$weight)
  
  # Penalty if the total weight surpass the limit
  total_weight <- if_else(total_weight > weightlimit, 0, total_weight)
  
  return(total_weight)
}

Let’s run the algorithm. We will use binary encoding since the decision variables are logical (TRUE or FALSE), with each bit represent a single logical value whether the item is selected.

The parameter of the algorithm:

• type: Type of the gene encoding, here we use “binary”
• fitness: The fitness function that will be optimized
• nBits: Number of bits generated for the decision variables
• pmutation: probability of mutation, default 0.1
• pcrossover: probability of crossover, default 0.8
• seed: The value of random seed to get reproducible result
• elitism: Number of best solution that will survive in each iteration, default is the top 5% of the population.
• maxiter: Number of max iterations for the algorithm to run
• popSize: Number of solutions in a population
• run: Number of iteration before the algorithm stop if there is no improvement on the optimal fitness value
• parallel: Will we use parallel computing? Can be logical value or the number of core used
• lower: The lowest value of the permutation encoding
• upper: The highest value of the permutation encoding
• keepBest: A logical whether we save best solution in each generation in a separate slot

There are several options that can be chosen for each of the GA process. For example, the default parent selection of GA for binary is linear rank selection. The detailed default methods used on each encoding can be found at here⁷.

The tic() and toc() function is used to measure the computation time.

Here we will run a genetic algorithm with binary encoding, the fitness function would be the evalFunc() we created earlier. The population size (popSize) is 100 with maximum iteration (maxiter) is also 100. The maximum run (run) with no improvement is 20, after that the algorithm would be stopped. The seed for reproducibility is set at 123. Since we use binary encoding, the parameter nBits control the number of bits generated. Since each bit would represent a decision for each item, the number of bits is equal to the number of row from the df_item_long.

tic()
library(GA)
ga_knapsack <- ga(type = "binary", fitness = evalFunc, popSize = 100, 
                  maxiter = 100, 
                  run = 20, nBits = nrow(df_item_long), seed = 123)
toc()

## 0.75 sec elapsed

The iteration stop at iterations 22 since after 20 iterations no improvement in the total weight. The optimal solution can be seen at Fitness function value which is 10000 (10,000 kg).

The Solution output is the solution that result in the optimal fitness value. We can choose one of them since they will result in the same fitness value. Let’s check one of them.

summary(ga_knapsack)

## -- Genetic Algorithm ------------------- 
## 
## GA settings: 
## Type                  =  binary 
## Population size       =  100 
## Number of generations =  100 
## Elitism               =  5 
## Crossover probability =  0.8 
## Mutation probability  =  0.1 
## 
## GA results: 
## Iterations             = 22 
## Fitness function value = 10000 
## Solutions = 
##      x1 x2 x3 x4 x5 x6 x7 x8 x9 x10  ...  x319 x320
## [1,]  1  1  0  1  1  1  1  0  1   1          0    1
## [2,]  1  1  0  1  1  1  1  0  1   1          1    1
## [3,]  0  1  1  1  1  0  0  1  1   0          0    0
## [4,]  1  1  0  1  1  1  1  0  1   1          1    1
## [5,]  1  1  0  1  1  1  1  0  1   1          1    0

The worst solution in the population has fitness value of 0 (those who has penalized because it passed the weight limit constrain), while the optimal solution has fitness value of 10000. From all of 100 chromosomes or solutions generated, the median of fitness value is 9572 while the mean is 6733.

To check the summary, we use summary() function on GA object

summary(ga_knapsack)

## -- Genetic Algorithm ------------------- 
## 
## GA settings: 
## Type                  =  binary 
## Population size       =  100 
## Number of generations =  100 
## Elitism               =  5 
## Crossover probability =  0.8 
## Mutation probability  =  0.1 
## 
## GA results: 
## Iterations             = 22 
## Fitness function value = 10000 
## Solutions = 
##      x1 x2 x3 x4 x5 x6 x7 x8 x9 x10  ...  x319 x320
## [1,]  1  1  0  1  1  1  1  0  1   1          0    1
## [2,]  1  1  0  1  1  1  1  0  1   1          1    1
## [3,]  0  1  1  1  1  0  0  1  1   0          0    0
## [4,]  1  1  0  1  1  1  1  0  1   1          1    1
## [5,]  1  1  0  1  1  1  1  0  1   1          1    0

As we can see, there are several solutions for the problem, all of them have the same fitness value: 10000. Therefore, we have the flexibility to choose which solution we will pick. All optimal solution has no value greater than its constraint, which is also 10000.

Let’s check the plot of fitness value progression using plot().

plot(ga_knapsack)

We can see that the optimal solution (Best) is already achieved at iteration 3.

Let’s choose one of the solution as our decision.

df_sol <- df_item_long[ga_knapsack@solution[1, ] == 1, ]

df_sol <- df_sol %>%  
  group_by(item, weight) %>% 
  summarise(freq = n()) %>% 
  mutate(total_weigth = freq * weight)

rmarkdown::paged_table(df_sol)

Let’s check the total weight

sum(df_sol$total_weigth)

## [1] 10000

The optimal fitness value is same with the constraint, so we can fully load the vehicle based on it’s maximum capacity.

Let’s check the third solution

df_sol <- df_item_long[ga_knapsack@solution[3, ] == 1, ]

df_sol <- df_sol %>%  
  group_by(item, weight) %>% 
  summarise(freq = n()) %>% 
  mutate(total_weigth = freq * weight)

rmarkdown::paged_table(df_sol)

sum(df_sol$total_weigth)

## [1] 10000

Even though the first solution and the second solution has different decisions, the fitness function or the total weights are the same, therefore both are optimal solutions and the decision to choose which one should be put into actual use is being handed to the decision maker.

Traveling Salesman Problem

The traveling salesperson problem (TSP) is one of the most widely discussed problems in combinatorial optimization. In its simplest form, consider a set of n cities with known symmetric intra-distances, the TSP involves finding an optimal route for visiting all the cities and return to the starting point such that the distance traveled is minimized. The set of feasible solutions is given by the total number of possible routes, which is equal to (n − 1)!/2.

Consider this eurodist dataset. It consists of 21 x 21 matrix of distance between cities. The total number of possible routes would be (21 - 1)!/2 = 1216451004088320000. Trying each one of them would take a great amount of time. Simple random search may fail to find the global optima.

cities <- as.matrix(eurodist)

head(cities)

##           Athens Barcelona Brussels Calais Cherbourg Cologne Copenhagen Geneva
## Athens         0      3313     2963   3175      3339    2762       3276   2610
## Barcelona   3313         0     1318   1326      1294    1498       2218    803
## Brussels    2963      1318        0    204       583     206        966    677
## Calais      3175      1326      204      0       460     409       1136    747
## Cherbourg   3339      1294      583    460         0     785       1545    853
## Cologne     2762      1498      206    409       785       0        760   1662
##           Gibraltar Hamburg Hook of Holland Lisbon Lyons Madrid Marseilles
## Athens         4485    2977            3030   4532  2753   3949       2865
## Barcelona      1172    2018            1490   1305   645    636        521
## Brussels       2256     597             172   2084   690   1558       1011
## Calais         2224     714             330   2052   739   1550       1059
## Cherbourg      2047    1115             731   1827   789   1347       1101
## Cologne        2436     460             269   2290   714   1764       1035
##           Milan Munich Paris Rome Stockholm Vienna
## Athens     2282   2179  3000  817      3927   1991
## Barcelona  1014   1365  1033 1460      2868   1802
## Brussels    925    747   285 1511      1616   1175
## Calais     1077    977   280 1662      1786   1381
## Cherbourg  1209   1160   340 1794      2196   1588
## Cologne     911    583   465 1497      1403    937

Since the objective is to minimize the distance, the fitness value would be in negative.

tsp_fitness <- function(x) {
  x <- c(x, x[1])
  
  # create from-to matrix 
  route <- embed(x, 2)[, 2:1]
 
 # Calculate distance
 distance <- -sum(cities[route])
 return(distance)
 }

Let’s run the algorithm.

tic()
ga_cities <- ga(type = "permutation", fitness = tsp_fitness, monitor = F,
         lower = 1, upper = nrow(cities), popSize = 100, maxiter = 5000, 
         run = 500, seed = 123)
summary(ga_cities)

## -- Genetic Algorithm ------------------- 
## 
## GA settings: 
## Type                  =  permutation 
## Population size       =  100 
## Number of generations =  5000 
## Elitism               =  5 
## Crossover probability =  0.8 
## Mutation probability  =  0.1 
## 
## GA results: 
## Iterations             = 1739 
## Fitness function value = -12919 
## Solutions = 
##      x1 x2 x3 x4 x5 x6 x7 x8 x9 x10  ...  x20 x21
## [1,] 15  2  9 12 14  5 18  4  3  11         8  13
## [2,]  5 18  4  3 11  7 20 10  6  17        12  14
## [3,]  1 21 17  6 10 20  7 11  3   4        16  19
## [4,]  6 10 20  7 11  3  4 18  5  14        21  17
## [5,]  1 19 16  8 13 15  2  9 12  14        17  21

toc()

## 14.99 sec elapsed

We have 5 solutions/routes that has the minimal distance. Transform the information into readable route

route_optimal <- ga_cities@solution

cities_name <- names(cities[1 , ])

cities_name <- apply(route_optimal, 1, function(x) {cities_name[x]})

apply(cities_name, 2, paste, collapse = " > ")

## [1] "Marseilles > Barcelona > Gibraltar > Lisbon > Madrid > Cherbourg > Paris > Calais > Brussels > Hook of Holland > Copenhagen > Stockholm > Hamburg > Cologne > Munich > Vienna > Athens > Rome > Milan > Geneva > Lyons"
## [2] "Cherbourg > Paris > Calais > Brussels > Hook of Holland > Copenhagen > Stockholm > Hamburg > Cologne > Munich > Vienna > Athens > Rome > Milan > Geneva > Lyons > Marseilles > Barcelona > Gibraltar > Lisbon > Madrid"
## [3] "Athens > Vienna > Munich > Cologne > Hamburg > Stockholm > Copenhagen > Hook of Holland > Brussels > Calais > Paris > Cherbourg > Madrid > Lisbon > Gibraltar > Barcelona > Marseilles > Lyons > Geneva > Milan > Rome"
## [4] "Cologne > Hamburg > Stockholm > Copenhagen > Hook of Holland > Brussels > Calais > Paris > Cherbourg > Madrid > Lisbon > Gibraltar > Barcelona > Marseilles > Lyons > Geneva > Milan > Rome > Athens > Vienna > Munich"
## [5] "Athens > Rome > Milan > Geneva > Lyons > Marseilles > Barcelona > Gibraltar > Lisbon > Madrid > Cherbourg > Paris > Calais > Brussels > Hook of Holland > Copenhagen > Stockholm > Hamburg > Cologne > Munich > Vienna"

Better representation for the route can be visualized through network or graph.

Production Scheduling

Supposed that we have various units of car to be produced at our manufacturing plant. There are 7 different types of product colors. If the following unit has different color than the previous car, there will changeover cost (cost incurred due to change in product types), such as the cost of material required to clean the equipment from the previous paint color. There are many way to optimize the color sequence, but we will use the simplest one in here for illustration. We want to minimize the number of setup, making the number of color switch or changover occur less, therefore the changeover cost can also be minimized.

First we import the data.

df_car <- read.csv("data_input/car_data.csv")

rmarkdown::paged_table(df_car)

The objective function will be: \[ S = 1 + \Sigma_{k=2}^{DT} \ S_k\]

• \(S\) = number of setup
• \(D_T\) = Number of car/unit to be produced
• \(k\) = The position of car in the sequence
• \(S_k\) = Is setup required? 1 if the next car (\(K+1\)) doesn’t have the same color with the current \(k^{th}\) car, 0 if the next car does

We will translate the objective function into an R function of min_setup. The algorithm in GA package run in context of maximization, so in order to do a minimization problem, we will make our fitness value negative.

min_setup <- function(x){
  df <- df_car[x, ]
  print(df)
  S <- df %>% 
  mutate(color_next = lead(color),
         setup = if_else(color == color_next, 0, 1)) %>% 
  head(nrow(df)-1) %>% 
  pull("setup") %>% 
  sum() + 1
  S <- -S
  return(S)
}

Now we run the algorithm. Since the problem is a scheduling problem, we will encode our solutions with permutation encoding. We will stop the algorithm if the number of iteration has reached 1000 iterations or if in 100 iterations there is no improvement in the optimal fitness value.

IMPORTANT The choice of population size and the number of iteration will greatly affect the performance of GA. If the population size is too little, there will be limited gene pool in our population so the the new solutions will not differ too much from the previous generations. If the number of iterations is too little, there will not be enough time for the algorithm to find new optimal solutions. GA package can be run with parallel computing.

tic()
ga_schedule <- ga(type = "permutation", fitness = min_setup, lower = 1, upper = nrow(df_car), 
                  seed = 123, maxiter = 1000, popSize = 100, elitism = 50, monitor = F, 
                  run = 100, parallel = T)

summary(ga_schedule)

## -- Genetic Algorithm ------------------- 
## 
## GA settings: 
## Type                  =  permutation 
## Population size       =  100 
## Number of generations =  1000 
## Elitism               =  50 
## Crossover probability =  0.8 
## Mutation probability  =  0.1 
## 
## GA results: 
## Iterations             = 223 
## Fitness function value = -6 
## Solutions = 
##       x1 x2 x3 x4 x5 x6 x7 x8 x9 x10  ...  x27 x28
## [1,]   1 27  8  2 23  3 16 20 12  15        14  22
## [2,]  14 22 21  6 12 15 16 27  8   2        19   9
## [3,]  14 22 21  6 12 15 16  1 27   3        18  11
## [4,]  21  6 12 15 27  2 16  8 23   1        14  22
## [5,]  21  6 12 15  1  3 23  8 27  16        18  11
## [6,]  14 22 21 12  6 15 27  2 20   1        19   9
## [7,]   9 19 14 22 21  6 12 15  7  28        11  18
## [8,]  21  6 12 15  1 27  2 23  3  20        14  22
## [9,]  21  6 12 15 16  1 27  3 20   8        14  22
## [10,] 21  6 12 15  2 23 16 20 27   1        14  22
##  ...                                              
## [55,] 21 12  6 15  1  3 27 16 20   8        14  22
## [56,] 16  2 27 23  8  1  3 20 21   6        14  22

toc()

## 54.41 sec elapsed

We can plot the progression of the algorithm by using plot() function on the GA object.

plot(ga_schedule)

The best or optimal solutions are obtained at least at the 102nd iteration and didn’t improve since then, so the algorithm stop after 100 iterations at 233nd iteration.

Let’s put the optimal sequence to our data. The resulting data frame will be our production schedule. Since there are multiple solutions, let’s just look at the first two optimal solutions.

The first solution:

rmarkdown::paged_table(df_car[ga_schedule@solution[1, ], ])

The second solution:

rmarkdown::paged_table(df_car[ga_schedule@solution[2, ], ])

Even though the first solution and the second solution has different sequence, the fitness function or the number of setup are the same, therefore both are optimal solutions and the decision to choose which one should be put into actual use is being handed to the decision maker, in these case, the production control staff.

There are still many application of GA in various field, especially in manufacturing and engineering process. As long as there is a fitness function and decision variables, GA can be appllied toward the problem.

Import Data

attrition <- read.csv("data_input/attrition.csv")

attrition <- attrition %>% 
  mutate(job_involvement = as.factor(job_involvement),
         education = as.factor(education),
         environment_satisfaction = as.factor(environment_satisfaction),
         performance_rating = as.factor(performance_rating),
         relationship_satisfaction = as.factor(relationship_satisfaction),
         work_life_balance = as.factor(work_life_balance),
         attrition = factor(attrition, levels = c("yes", "no")))

rmarkdown::paged_table(attrition)

Cross-Validation

We split the data into training set and testing dataset, with 80% of the data will be used as the training set.

set.seed(123)
intrain <- initial_split(attrition, prop = 0.8, strata = "attrition")

Data Preprocessing

We will preprocess the data with the following step:

• Upsample to increase the number of the minority class and prevent class imbalance
• Scaling all of the numeric variables
• Remove numeric variable with near zero variance
• Use PCA to reduce dimensions with threshold of 90% information kept
• Remove over_18 variable since it only has 1 levels of factor

rec <- recipe(attrition ~ ., data = training(intrain)) %>% 
  step_upsample(attrition) %>% 
  step_scale(all_numeric()) %>% 
  step_nzv(all_numeric()) %>% 
  step_pca(all_numeric(), threshold = 0.9) %>% 
  step_rm(over_18) %>% 
  prep()

data_train <- juice(rec)
data_test <- bake(rec, testing(intrain))

Create Fitness Function

We will train the model using random forest. We will optimize the mtry (number of nodes in each tree) and the min_n (minimum number of member on each tree in order to be splitted) value by maximizing the model accuracy on the testing dataset. Since the number of mtry and min_n is integer, we better use binary encoding instead of real-valued encoding.

We will set the range of mtry from 1 to 32. We will limit the range of min_n from 1 to 128. Let’s check how many bits we need to cover those numbers.

# max value of mtry = 32
2^5

## [1] 32

# max value of min_n = 128
2^7

## [1] 128

So, we will need at least 12 bits of binary. If the converted value of bits for mtry or min_n is 0, we change it to 1.

We will write the model as fitness function. We will maximize the accuracy of our model on the testing dataset.

fit_function <- function(x){
  a <- binary2decimal(x[1:5])
  b <- binary2decimal(x[6:12])
  
  if (a == 0) {
    a <- 1
  }
  if (b == 0) {
    b <- 1
  }
  if (a > 17) {
    a <- 17
  }
  
#define model spec
model_spec <- rand_forest(
  mode = "classification",
  mtry = a,
  min_n = b,
  trees = 500)

#define model engine
model_spec <- set_engine(model_spec,
                         engine = "ranger",
                         seed = 123,
                         num.threads = parallel::detectCores(),
                         importance = "impurity")

#model fitting
set.seed(123)
model <- fit_xy(
  object = model_spec,
  x = select(data_train, -attrition),
  y = select(data_train, attrition)
)

pred_test <- predict(model, new_data = data_test %>% select(-attrition))
acc <- accuracy_vec(truth = data_test$attrition, estimate = pred_test$.pred_class)
return(acc)
}

Run the Algorithm

IMPORTANT Since some model require a lot of time to train, you may need to be extra careful on choosing the population size and the number of iteration. If you don’t have a problem to run the GA on a long period of time, than you may choose bigger population size and number of iterations. Here I will set the population size to be only 100 and the maximum number of iteration is 100. You can use parallel = TRUE if you want faster computation with parallel computing.

We will try to use tournament selection instead of the default linear rank selection.

tic()
ga_rf <- ga(type = "binary", fitness = fit_function, nBits = 12, seed = 123, 
            popSize = 100, maxiter = 100, run = 10, parallel = T, 
            selection = gabin_tourSelection)
summary(ga_rf)

## -- Genetic Algorithm ------------------- 
## 
## GA settings: 
## Type                  =  binary 
## Population size       =  100 
## Number of generations =  100 
## Elitism               =  5 
## Crossover probability =  0.8 
## Mutation probability  =  0.1 
## 
## GA results: 
## Iterations             = 18 
## Fitness function value = 0.8805461 
## Solutions = 
##      x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12
## [1,]  0  0  0  0  1  0  0  1  0   0   1   0
## [2,]  0  0  0  0  0  0  0  1  0   0   1   0

toc()

## 478.7 sec elapsed

Let’s plot the progression of GA

plot(ga_rf)

Let’s check each decision variables and the accuracy that is used as the fitness value.

# Transform decision variables into data frame
ga_rf_pop <- as.data.frame(ga_rf@population)

# Convert decision variables into mtry and min_n
ga_mtry <- apply(ga_rf_pop[ , 1:5], 1, binary2decimal)
ga_min_n <- apply(ga_rf_pop[ , 6:12], 1, binary2decimal)

# Show the list of solutions in the final population
rmarkdown::paged_table(
  data.frame(mtry = ga_mtry,
           min_n = ga_min_n,
           accuracy = ga_rf@fitness) %>% 
  mutate(mtry = ifelse(mtry == 0, 1, mtry),
         min_n = ifelse(min_n == 0, 1, min_n),
         mtry = ifelse(mtry > 17, 17, mtry)) %>% 
  distinct() %>%
  arrange(desc(accuracy))
  )

Compare with K-fold Cross Validation

You may want to compare the computation time of GA with the required time for repeated K-fold cross-validation using the conventional caret package. We will use k = 5 and 10 repeats.

tic()
set.seed(123)
ctrl <- trainControl(method="repeatedcv", number=5, repeats=10)
attrition_forest <- train(attrition ~ ., data=data_train, method="rf", trControl = ctrl)
toc()

## 702.37 sec elapsed

attrition_forest

## Random Forest 
## 
## 1974 samples
##   17 predictor
##    2 classes: 'yes', 'no' 
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 10 times) 
## Summary of sample sizes: 1580, 1579, 1578, 1580, 1579, 1578, ... 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa    
##    2    0.9462064  0.8924131
##   22    0.9587116  0.9174242
##   42    0.9475654  0.8951309
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 22.

Let’s check the accuracy of the random forest model with caret

pred_rf <- predict(attrition_forest, newdata = data_test)
accuracy_vec(truth = data_test$attrition, estimate = pred_rf)

## [1] 0.8634812

The model performance by GA is better than those optimized with simple K-fold cross validation, even though GA require more time to compute. This trade-off of computation time vs performance should be considered.