4.1.1 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. 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. I’ve prepared the dummy data which you can access here.

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

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, ]
    for (i in 1:(nrow(df_car) - 1)) {
        df$s[i] <- if_else(df$color[i] == df$color[i + 1], 0, 1)
    }
    S = -(1 + sum(df$s, na.rm = T))
    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 300 iterations or if in 30 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.

The parameter of the algorithm:

• type: Type of the gene encoding, here we use “permutation”
• fitness: The fitness function that will be optimized
• lower: The lowest value of the permutation encoding, we will encode the permutation from 1 to the number of car available
• upper: The highest value of the permutation encoding
• 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, here we set it to 50 solutions that will survive
• 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

The rest of the parameters, such as the crossover and mutation probability are left to the default setting.

tic()
gann <- ga(type = "permutation", fitness = min_setup, lower = 1, upper = nrow(df_car), 
    seed = 123, elitism = 50, maxiter = 300, popSize = 300, run = 30, parallel = parallel::detectCores())

summary(gann)

-- Genetic Algorithm ------------------- 

GA settings: 
Type                  =  permutation 
Population size       =  300 
Number of generations =  300 
Elitism               =  50 
Crossover probability =  0.8 
Mutation probability  =  0.1 

GA results: 
Iterations             = 96 
Fitness function value = -7 
Solutions = 
      x1 x2 x3 x4 x5 x6 x7 x8 x9 x10  ...  x27 x28
[1,]  25  4 11 18  1 15  3  8 17  22        13  27
[2,]  19  7 21 10  5 24  3 17 15   2        11  18
[3,]  25 11  4 18  1  8  3 15  2  22        13  27
[4,]  25  4 11 18  1  8  3 15  2  22        13  27
[5,]  25  4 11 18  2  1 15 22 16  17        19   5
[6,]   1  2  8 16 15 17  3 22 23   9        11  18
[7,]  25  4 11 18  1  8  3 17 15   2        13  27
[8,]  11  4 25 18  1  8  3 15  2  22        13  27
[9,]   7 19 24 10 21  5  1 15  3   8        23   9
[10,]  4 25 11 18  2  1 15 22 16  17        19   5
 ...                                              
[24,]  4 25 11 18  1  8  3 17 15   2        13  27
[25,]  2  1 15 22 16 17  3  8 23   9        11  18

toc()

87.14 sec elapsed

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

plot(gann)

The best or optimal solutions are obtained at least at the 67th iteration and didn’t improve since then, so the algorithm stop after 30 iterations at 96th iteration.

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.

min_setup(gann@solution[1, ])

[1] -7

If you want to look at all of the fitness value from the last population, use can use the fitness attribute from the GA object. We’ll summarise the fitness value to get the distribution.

summary(gann@fitness)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -25.00  -18.00  -15.00  -14.29  -11.00   -7.00 

The worst solution in the population has fitness value of 25, while the optimal solution has fitness value of 7. From all 300 chromosomes or solutions generated, the median of fitness value is 15 while the mean is 14.29.

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:

df_car[gann@solution[1, ], ]

The second solution:

df_car[gann@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.

4.1.2 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))

df_item

Each items have quantities and their respective weight.s

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(1, 3)]

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\]

weightlimit <- 10000

evalFunc <- function(x) {
    df <- df_item_long[x == 1, ]
    total_weight <- sum(df$weight)
    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 integer-valued. Each bit represent a single logical value whether the item is selected.

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

summary(gann2)

-- 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

toc()

0.7 sec elapsed

Let’s check the plot.

plot(gann2)

Let’s choose one of the solution as our optimal choice.

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

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

df_sol

Let’s check the total weight

sum(df_sol$total_weigth)

[1] 10000

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

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.

4.2.1 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

4.2.2 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")

intrain

<1177/293/1470>

4.2.3 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))

data_train

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, trees = 500, min_n = b)
    
    # 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)
}

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.

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 https://www.rdocumentation.org/packages/GA/versions/3.2/topics/gaControl. Here, 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             = 16 
Fitness function value = 0.8805461 
Solutions = 
     x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12
[1,]  0  0  0  0  0  0  1  0  0   0   0   0
[2,]  0  0  0  0  1  0  1  0  0   0   0   0

toc()

633.28 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
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))

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()

660.06 sec elapsed

attrition_forest

Random Forest 

1974 samples
  17 predictor
   2 classes: 'no', 'yes' 

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.9489451  0.8978902
  22    0.9596744  0.9193495
  42    0.9476663  0.8953330

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.