Enhancing Digital Twins Part 3: Predictive Maintenance with Azure Databricks

 

%python //or r, scala, sql
 

 

# load Apache Spark's R interface through the ‘sparklyr’ package 

library(sparklyr, quietly = TRUE) 

# then establish a Databricks connection to Apache Spark using 'sparklyr' 

sc <- spark_connect(method = "databricks") 

# The connection is assigned to var 'sc' 
 

 

# Read our csv file into a Spark dataframe using spark_read_csv
# The first parameter must be a reference to our spark connection variable, 'sc'
# https://www.rdocumentation.org/packages/sparklyr/versions/1.0.2/topics/spark_read_csv 

sdf <- spark_read_csv(sc, name = "test_maintenance", path = "/FileStore/tables/maintenancedata.csv", memory = FALSE, infer_schema = TRUE) 
 

 

# Create a subset ('broken_set') of your Spark dataframe where broken == 1 for each row 
# The subset is created using the dplyr package pipe '%>%', which can be read as 'Then';
#(ensure the dplyr package is loaded in Cmd 1)
# The pipe can be thought of like: you take the sdf dataframe, 'Then' you filter the data by broken == 1. 

broken_set <- sdf %>% filter(broken == 1)

# https://www.rdocumentation.org/packages/graphics/versions/3.6.1/topics/layout 
# The following 3 boxplots will be positioned according to layout() 

layout(matrix(c(1,1,2,2, 

              1,1,2,2, 

              1,1,2,2, 

              3,3,3,3, 

              3,3,3,3, 

              3,3,3,3), ncol  = 4, byrow = TRUE))  

# https://www.rdocumentation.org/packages/graphics/versions/3.6.1/topics/boxplot
# For the boxplot() command, a quantitative variable is connected to a qualitative variable using the tilde ‘~’;
#You must then specify the dataset as an additional argument to the function. 

boxplot(lifetime~broken, data=broken_set, main="Broken machines", xlab="", ylab="Lifetime", col="#357EC7") 

boxplot(lifetime~branch, data=broken_set, main="Per Branch", xlab="branch", ylab="Lifetime", col="#357EC7") 

boxplot(lifetime~machine, data=broken_set, main="Per Machine", xlab="machine", ylab="Lifetime", col="#357EC7") 
 

 

# To calculate a survival analysis, we will need to use the survival package's survfit() function;
#(ensure the 'survival' package is loaded in Cmd 1)
# https://www.rdocumentation.org/packages/survival/versions/2.11-4/topics/survfit
# survfit(survivalObject, dataset), to add a grouping to the survivalObject, add it to the right of the tilde ‘~’, use ~1 if there is none. 

maintenance_graph <- survfit(Surv(lifetime,broken) ~ 1, data = sdf) 

# https://www.rdocumentation.org/packages/GGally/versions/1.4.0/topics/ggsurv 

ggsurv(maintenance_graph, cens.size = 1, cens.shape = 4) 
 

 

maintenance_graph2 <- survfit(Surv(lifetime,broken) ~ branch, data = sdf) 

ggsurv(maintenance_graph2, cens.size = 1, cens.shape = 4) 
 

 

maintenance_graph3 <- survfit(Surv(lifetime,broken) ~ machine, data = sdf) 

ggsurv(maintenance_graph3, cens.size = 1, cens.shape = 4) 
 

 

# Perform a Survival regression analysis (SRA) on only Branch A
# To do an SRA, we will need to use the survival package's survreg() function;
#(ensure the 'survival' package is loaded in Cmd 1)  

# Filter the data for only rows where branch = A  

branchA_set <- sdf %>% filter(branch == "Branch A") 

# Extract the columns 'pressure' and 'data' as vectors to force column values to be of 1 type, else Surv() will fail 

pressure = pull(branchA_set, pressure) 

broken = pull(branchA_set, broken) 

# Choose lifetime and broken as dependant variables to be used in the survival regression model. 

dependantvars <- Surv(pressure, broken) 

# Create model (use the gaussian method) 

survreg <- survreg(dependantvars~pressure+moisture+temperature+machine, dist="gaussian",data=branchA_set) 

summary(survreg) 
 

 

# p = percentile = 0.5 = expected median
# Analyse Branch A only 
# Ebreak : Estimate of how long a machine will live 

Ebreak<-predict(survreg, newdata=branchA_set, type="quantile", p=.5) 

# Extract the column 'lifetime' as a vector to force column values to be of 1 type 

lifetime <- pull(branchA_set, lifetime)

# Make forecast dataset, coerce columns from branchA_set to all be 1 type using pull() 

Forecast<-data.frame(Ebreak) 

# Lifetime: How long a machine has lived 

Forecast$lifetime=lifetime 

Forecast$broken=pull(branchA_set, broken) 

Forecast$id=pull(branchA_set, id) #Also, we attach the ID of the machines. 

# Computed Expected Remaining Lifetime of mmachine (remainingLT) 

Forecast$RemainingLT<-Forecast$Ebreak-lifetime 

# Order the elements by Expected Remaining Lifetime 

Forecast<-Forecast[order(Forecast$RemainingLT),] 

# Keep only those who are not broken yet 

ActionsPriority <- Forecast %>% filter(broken ==0) 

print(ActionsPriority) 
 

 

ActionsPriority$class <- cut(ActionsPriority$RemainingLT, c(-10,1,4,1000)) 

levels(ActionsPriority$class) <- c('urgent', 'medium', 'good') 

summary(select(ActionsPriority,-c(id)))  #leave out the id column 
 

 

# Outputs for the digital twin (Warehouse scenario = Branch A) 
# Get the information for the asset maintenance overview board 
# Total asset count 

totalBranchA <- branchA_set %>% tally()  

# Total broken asset count 

total_Broken_branchA <- branchA_set %>% filter(broken == 1) %>% tally() 

# Total working asset count 

total_Working_branchA <- branchA_set %>% filter(broken == 0) %>% tally() 

# This extracts values from the tally result table by column name 'n' 

totalBranchA <- pull(totalBranchA, n)
total_Broken_branchA <- pull(total_Broken_branchA, n)  
total_Working_branchA <-pull(total_Working_branchA, n)  

# Maintenance_class_counts = count(ActionsPriority, class) 
# To get each count separately 

urgent <- length(which(ActionsPriority$class=='urgent')) 
medium <- length(which(ActionsPriority$class=='medium'))
good <- length(which(ActionsPriority$class=='good'))