The four pillars:
Prescriptive analytics seeks to provide a data-driven course of action, that is, recommendations that are likely to result in the best outcome as determined from the data. These recommendations may take the form of decision rules, prioritized lists of actions, or specific prescriptions for how to allocate resources or make trade-offs.
Let’s make use of our prediction models to identify days with unusual downtimes.
tb_subset = tb %>%
filter(day >= as.Date("2020/11/01")) %>%
mutate(uid = paste0(plant, "_", machine_id))
tb_formatted_data = tb_subset %>%
group_by(day, plant) %>%
summarize(tot_dt = sum(dt), dt_instances = n()) %>%
ungroup() %>%
group_by(plant) %>%
complete(day = seq(min(tb_subset$day), max(tb_subset$day), by = 'day')) %>%
mutate(wk = epiweek(day)) %>%
mutate(yr = year(day)) %>%
mutate(wid = paste0(yr, "_", wk)) %>%
replace_na(list(tot_dt = 0, dt_instances = 0)) %>%
mutate(wday = weekdays(day))
tb_reg_data = tb_formatted_data %>%
dummy_cols(select_columns = c("plant", "wday")) %>%
select(dt_instances, tot_dt, day, wid, plant, wday, 9:18) %>%
mutate(l1_dt_instances = lag(dt_instances, 1)) %>%
mutate(l7_dt_instances = lag(dt_instances, 7)) %>%
drop_na() %>%
filter(dt_instances > 5)
cat_nms = names(tb_reg_data)[7:16]
cat_vars = expand_grid(plant = cat_nms[1:3], wday = cat_nms[4:10]) %>%
mutate(vars = paste0(plant, ":", wday)) %>%
pull(vars) %>%
paste0(collapse = " + ")
tot_inst_formula = paste0(
"dt_instances ~ ", cat_vars, "+ l1_dt_instances + l7_dt_instances - 1"
) %>% as.formula()
mod_inst = lm(tot_inst_formula, tb_reg_data)
tb_inst_results = tb_reg_data %>%
mutate(residuals = mod_inst$residuals)
tb_inst_results %>%
mutate(formatted_date = paste0(wday, ", ", format.Date(day))) %>%
plot_ly(
color = ~plant, x = ~residuals, text = ~formatted_date, type = "box",
boxpoints = "all", quartilemethod = "inclusive",
pointpos = 0, jitter = 0.5, alpha = 0.5
)
tb_reg_data2 = tb_formatted_data %>%
dummy_cols(select_columns = c("plant", "wday")) %>%
select(dt_instances, tot_dt, day, wid, plant, wday, 9:18) %>%
mutate(l1_tot_dt = lag(tot_dt, 1)) %>%
mutate(l7_tot_dt = lag(tot_dt, 7)) %>%
drop_na() %>%
filter(dt_instances > 5)
tot_dt_formula = paste0(
"tot_dt ~ ", cat_vars, "+ l1_tot_dt + l7_tot_dt - 1"
) %>% as.formula()
mod_tot_dt = lm(tot_dt_formula, tb_reg_data2)
tb_tot_dt_results = tb_reg_data2 %>%
mutate(residuals = mod_tot_dt$residuals)
tb_tot_dt_results %>%
mutate(formatted_date = paste0(wday, ", ", format.Date(day))) %>%
plot_ly(
color = ~plant, x = ~residuals, text = ~formatted_date, type = "box",
boxpoints = "all", quartilemethod = "inclusive",
pointpos = 0, jitter = 0.5, alpha = 0.5
)
# Grab all outliers above the mean
plants = tb_inst_results$plant %>% unique()
tb_date_outliers = lapply(plants, function(p) {
inst_vec = tb_inst_results %>%
filter(plant == p) %>%
pull(residuals)
totd_vec = tb_tot_dt_results %>%
filter(plant == p) %>%
pull(residuals)
inst_out = boxplot.stats(inst_vec)$out[boxplot.stats(inst_vec)$out > mean(inst_vec)]
totd_out = boxplot.stats(totd_vec)$out[boxplot.stats(totd_vec)$out > mean(totd_vec)]
tmp1 = tb_inst_results %>%
filter(plant == p) %>%
filter(residuals %in% inst_out) %>%
select(day, wday, plant, dt_instances:tot_dt)
tmp2 = tb_tot_dt_results %>%
filter(plant == p) %>%
filter(residuals %in% totd_out) %>%
select(day, wday, plant, dt_instances:tot_dt)
bind_rows(tmp1, tmp2) %>%
mutate(day = format.Date(day))
}) %>% bind_rows() %>%
distinct()
tb_date_outliers %>%
DT::datatable()
We’ll use our behavior models from the previous section to determine outliers.
First let’s look at downtime instances.
unique_states = tb_subset$reason_code_id %>% unique()
unique_machines = tb_subset %>%
group_by(uid) %>%
summarize(n = n()) %>%
filter(n > 100) %>%
pull(uid)
tm_list = lapply(unique_machines, function(u_id) {
state_seq = tb_subset %>%
filter(uid == u_id) %>%
pull(reason_code_id)
mc_fit = markovchainFit(data = state_seq, possibleStates = unique_states)
return(mc_fit$estimate@transitionMatrix)
}) %>% setNames(unique_machines)
# Calculate the distance matrix
tb_machine_dist = expand_grid(row = unique_machines, col = unique_machines) %>%
rowwise() %>%
mutate(d = norm(tm_list[[row]] - tm_list[[col]], 'F'))
# Create the distance matrix
M = tb_machine_dist %>%
pivot_wider(names_from = col, values_from = d) %>%
select(-row) %>%
as.matrix() %>%
magrittr::set_rownames(colnames(.))
# Get the hierarchical tree
M_dist = as.dist(M)
M_tree = hclust(M_dist)
M_dend = as.dendrogram(M_tree)
K = 4
clusters = cutree(M_dend, k=K, order_clusters_as_data = FALSE)
tb_clusters = tibble(uid = names(clusters), cluster = clusters)
uid_order = tb_clusters$uid
tb_per_machine = tb_subset %>%
group_by(uid, plant, category, description) %>%
summarize(dt_instances = n(), avg_dt = mean(dt), med_dt = median(dt), tot_dt = sum(dt)) %>%
filter(dt_instances > 100) %>%
ungroup() %>%
left_join(tb_clusters, by = 'uid') %>%
mutate(cluster = factor(cluster))
tb_per_machine %>%
mutate(detail_description = paste0(description, " (", category, ")")) %>%
mutate(cluster = fct_reorder(cluster, tot_dt)) %>%
plot_ly(
color = ~cluster, x = ~dt_instances, type = "box",
boxpoints = "all", text = ~detail_description, quartilemethod = "inclusive",
pointpos = 0, jitter = 0.5, alpha = 0.5
) %>%
layout(
xaxis = list(title = "downtime instances"),
yaxis = list(title = "behavior cluster"),
title = "downtime instances per machine"
)
Next, let’s look at total downtime.
tb_per_machine %>%
mutate(detail_description = paste0(description, " (", category, ")")) %>%
mutate(cluster = fct_reorder(cluster, tot_dt)) %>%
mutate(tot_dt = tot_dt/(60*24)) %>%
plot_ly(
color = ~cluster, x = ~tot_dt, type = "box",
boxpoints = "all", text = ~detail_description, quartilemethod = "inclusive",
pointpos = 0, jitter = 0.5, alpha = 0.5
) %>%
layout(
xaxis = list(title = "total downtime (days)"),
yaxis = list(title = "behavior cluster"),
title = "total downtime per machine"
)
Let’s put this information together.
# Get all anomalies (outliers above the mean, cluster 1)
tb_anomalies = lapply(1:4, function(c) {
tmp = tb_per_machine %>%
filter(cluster == c)
out_inst = boxplot.stats(tmp$dt_instances)$out
out_inst = out_inst[out_inst > mean(tmp$dt_instances)]
out_totd = boxplot.stats(tmp$tot_dt)$out
out_totd = out_totd[out_totd > mean(tmp$tot_dt)]
if(c == 1) {
tmp
} else {
tmp %>%
filter(dt_instances %in% out_inst | tot_dt %in% out_totd)
}
}) %>% bind_rows() %>%
mutate(anomaly = T)
plants = c("all", tb_subset$plant %>% unique())
tb_tmp = lapply(plants, function(p) {
tmp = tb_subset %>%
mutate(uid = paste0(plant, "_", machine_id)) %>%
group_by(uid, category, plant) %>%
summarize(tot_dt = sum(dt)/(60*24), dt_inst = n()) %>%
mutate(group = p)
# To help in redraw of visualization
dummy_data = tibble(category = tb$category %>% unique()) %>%
mutate(tot_dt = 0) %>%
mutate(dt_inst = 0) %>%
mutate(plant = p)
if (p != "all") {
tmp = tmp %>% filter(plant == p)
}
tmp %>%
select(-plant) %>%
rename(plant = group) %>%
bind_rows(dummy_data) %>%
return()
}) %>% bind_rows() %>%
left_join(tb_anomalies %>% select(uid, anomaly), by = "uid") %>%
replace_na(list(anomaly = F)) %>%
left_join(tb_per_machine %>% select(uid, description, cluster), by = "uid") %>%
mutate(
detail_description = paste0(
description, " (", category, ")(cluster = ", cluster, ")"
)
)
tb_tmp %>%
plot_ly(x = ~tot_dt, y = ~dt_inst, text = ~detail_description) %>%
add_markers(data = tb_tmp, frame = ~plant, color = ~anomaly, alpha = 0.7) %>%
layout(
xaxis = list(
title = "total downtime (days)"
),
yaxis = list(
title = "total downtime instances"
),
title = "identifying anomalous behavior in machines"
) %>%
animation_button(visible = F) %>%
animation_opts(redraw = F, transition = 0) %>%
animation_slider(currentvalue = list(font = list(size = 0, color = "white")))
tb_anomalies %>%
select(-anomaly) %>%
DT::datatable()