knitr::include_graphics("fe0026434-image.png")

1 Faulty actuator of CO2 absorber

A main actuator of CO2 absorber found faulty due to broken stem and could not get fixed during online. Based on current history, the actuator still able to control flow of fluid, but opening of actuator sometime got higher and sometime got lower.

On the other hand, the plant have schedule for Turn Around (Shutdown and do some maintenance) of several equipment. The dilemma is, if the actuator fail before planned Turn Around, the plant will shutdown for long time and costly.

Hence, the prediction when this valve will fail became important to “decide” when the best time to schedule Turn Around. Hopefully, the plant could operate as maximum as possible to reach production target, delivery time for new set of actuator, manpower for maintenance but still in planned shutdown timeline.

The objective of this analysis is to predict when the actuator will reach near 100% opening. At 100% the valve can not hold the fluid and all fluid drained from the absorer tower and plant control will trigger to trip/shutdown.

2 Data Import

First, we have to import required library.

library(lubridate)
library(dplyr)
library(forecast)
library(plotly)
library(scales)
library(tseries)

Then, import the data. the actuator data consist of “Date” Column which contain date of previous date. “Percent opening” was a PV (Present value of Valve Opening). SV (Set value of valve opening) was a target of valve opening. and “Error” was the percent gap of PV and SV.

Fault_Actuator <- read.csv("Fault_Actuator3.csv")
str(Fault_Actuator)
## 'data.frame':    2190 obs. of  4 variables:
##  $ Date           : chr  "2022/01/01 00:00:00" "2022/01/01 01:00:00" "2022/01/01 02:00:00" "2022/01/01 03:00:00" ...
##  $ Percent.Opening: num  57 61 56.5 56.5 60.5 ...
##  $ SV             : int  50 50 50 50 50 50 50 50 50 50 ...
##  $ Error          : num  14 22 13 13.1 21.1 ...

2.1 Data Type Conversion

The date cariable shall be converted to datetime column, “ymd_hms” POSIXct format selected to as datetime format.

Fault_Actuator <- Fault_Actuator %>% 
  mutate(Date = ymd_hms(Date))

str(Fault_Actuator)
## 'data.frame':    2190 obs. of  4 variables:
##  $ Date           : POSIXct, format: "2022-01-01 00:00:00" "2022-01-01 01:00:00" ...
##  $ Percent.Opening: num  57 61 56.5 56.5 60.5 ...
##  $ SV             : int  50 50 50 50 50 50 50 50 50 50 ...
##  $ Error          : num  14 22 13 13.1 21.1 ...
tail(Fault_Actuator)
##                     Date Percent.Opening SV  Error
## 2185 2022-04-02 00:00:00          85.466 50 70.932
## 2186 2022-04-02 01:00:00          85.473 50 70.945
## 2187 2022-04-02 02:00:00          83.479 50 66.959
## 2188 2022-04-02 03:00:00          76.486 50 52.973
## 2189 2022-04-02 04:00:00          86.493 50 72.986
## 2190 2022-04-02 05:00:00          81.000 50 62.000

3 Frequency selection

Selected frequency used was 24, to see the seasonal data in daily interval.

Freq_1 = 24

The target of prediction is the percent opening of valve. Then, the percent opening converted into ts object.

Fault_Actuator_ts_1 <- ts(data = Fault_Actuator$Percent.Opening,
                        frequency = Freq_1)

head(Fault_Actuator_ts_1)
## Time Series:
## Start = c(1, 1) 
## End = c(1, 6) 
## Frequency = 24 
## [1] 57.007 61.014 56.521 56.527 60.534 58.541

4 Decompose Pattern with frequency with daily interval

To see the pattern of trend and seasonal data, the ts object then converted into decompose plot.

Fault_Actuator_ts_1 %>% 
  decompose () %>% 
  autoplot()

# Decompose Pattern with frequency with weekly interval

As the pattern of daily data not show straight pattern trend and contain seasonal event. Then, the frequency changed into weekly.

Freq_2 = 24*7
Fault_Actuator_ts_2 <- ts(data = Fault_Actuator$Percent.Opening,
                        frequency = Freq_2)
#,
head(Fault_Actuator_ts_2)
## Time Series:
## Start = c(1, 1) 
## End = c(1, 6) 
## Frequency = 168 
## [1] 57.007 61.014 56.521 56.527 60.534 58.541
Fault_Actuator_ts_2 %>% 
  decompose () %>% 
  autoplot()

Also as shown above, for frequency changed to weekly the trend still contain a seasonal pattern, we should increase the frequency one level to monthly interval.

5 Decompose Pattern with frequency with Monthly interval

Freq_3 = 24*7*4
Fault_Actuator_ts_3 <- ts(data = Fault_Actuator$Percent.Opening,
                        frequency = Freq_3)
#,
head(Fault_Actuator_ts_3)
## Time Series:
## Start = c(1, 1) 
## End = c(1, 6) 
## Frequency = 672 
## [1] 57.007 61.014 56.521 56.527 60.534 58.541
Fault_Actuator_ts_3 %>% 
  decompose () %>% 
  autoplot()

6 Decompose pattern of all frequency

As shown with frequency changed to monthly (2474), the trend almost close to liner and not contain seasonal factor anymore. Also, the seasonal already capture a lot of seasonal pattern

Fault_Actuator_msts <- msts(data = Fault_Actuator$Percent.Opening,
                            seasonal.periods = c(24, 24*7, 24*7*4))
Fault_Actuator_msts %>% 
  mstl() %>% 
  autoplot()

Figure above shows the different seasonal effect captured by for daily, weekly, monthly interval object ts. All seasonal pattern already captured in monthly interval.

Trend from montly interval is smooth and there is no seasonal left in trend. Hence, the ts object ready for future process. Then, we shall divide our data into train and test data. About 1518 data used as train data and 672 data used as test data.

Fault_Actuator_train <- head(Fault_Actuator_msts, -Freq_3)
Fault_Actuator_test <- tail(Fault_Actuator_msts, Freq_3)
length(Fault_Actuator_train) 
## [1] 1518
length(Fault_Actuator_test)
## [1] 672
Fault_Actuator_train %>% 
  autoplot()

Based on train data graph, the type of time series data more likely “additive” due to constant seasonal pattern. Hence, additive model used as parameter to specify the model.

7 Stationary Check

Then we should check our ts object already stationary or not with “adf.test”.

adf.test(Fault_Actuator_train)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  Fault_Actuator_train
## Dickey-Fuller = -5.5148, Lag order = 11, p-value = 0.01
## alternative hypothesis: stationary

the p-value from adf.test already 0.01 and below alpha 0.05, it means that object already stationer.

8 Holt Winter Exponential Smoothing

Holt’s Winters Exponential (Triple Exponential Smoothing) here used as initial trial and due to its simple use. Then create the model named “Fault_Actuator_hw_add”.

Fault_Actuator_hw_add <- HoltWinters(Fault_Actuator_train)

Then use the model to predict/forecast the next three months data.

Fault_Actuator_forecast_add <- forecast(Fault_Actuator_hw_add, h = 24*7*4*3)

The accuracy of HW model to predict the next data use data test to calculate the error produced from gap of data train and calculated data.

accuracy(Fault_Actuator_forecast_add$mean, Fault_Actuator_test)
##                ME    RMSE     MAE      MPE     MAPE      ACF1 Theil's U
## Test set 3.087498 8.25143 6.72383 3.194145 7.784647 0.4604425  1.430185

The MAPE of HW model 7.7% which quite good to fit the data.

Fault_Actuator_train %>% 
  autoplot() +
  autolayer(Fault_Actuator_forecast_add) 

max(Fault_Actuator_forecast_add$mean)
## [1] 112.7905

Maximum opening reached in the next three months maximum will reach 112.79% opening which beyond target 100%. The opening means will reach 100% before three months.

9 STLM Model

In order to reduce the MAPE, we can try alternative and then compare the the MAPE value. Here, the STLM used as advanced model which can capture data which can not captured by ordinary ETS and ARIMA, and preferably used as to model the additive type data.

Fault_Actuator_arima <- stlm( Fault_Actuator_train, s.window = Freq_3, method = "arima")

Same as HW model above, monthly interval used in STLM model and forecast date for next three months.

Fault_Actuator_forecast_arima <- forecast(Fault_Actuator_arima, h = 24*7*4*3)
accuracy(Fault_Actuator_forecast_arima$mean, Fault_Actuator_test)
##                ME     RMSE      MAE     MPE     MAPE      ACF1 Theil's U
## Test set 7.211447 9.752791 8.222783 8.06121 9.369746 0.4686135  1.670626

Unfortunately, the MAPE of STLM model is 9.36% which larger if compared with HW model.

Fault_Actuator_train %>% 
  autoplot() +
  autolayer(Fault_Actuator_forecast_arima)

STLM model trend not increased but stable as shown from graph above. The STLM model can not capture the trend and

So, here to predict the Holt’s Winters is better than than STLM to forecast based on MAPE value.

Then, Holt’s Winters used as basis to predict the percent opening actuator.

Forcasted_Date <- as.data.frame(seq(as.POSIXct("2022-04-02 05:00:00"), as.POSIXct("2022-06-25 04:00:00"), by = "hour"))
Forcasted_HW <- as.data.frame(Fault_Actuator_forecast_add)
Forcasted_Opening <- cbind(Forcasted_Date,Forcasted_HW)

colnames(Forcasted_Opening) <- c("Forcasted_Dates","Forcasted_Opening", "Lo_80", "Hi_80", "Lo_95", "Hi_95" )

head(Forcasted_Opening)
##              Forcasted_Dates Forcasted_Opening    Lo_80    Hi_80    Lo_95
## 3.258929 2022-04-02 05:00:00          86.94363 81.21120 92.67605 78.17664
## 3.260417 2022-04-02 06:00:00          77.94563 72.13898 83.75228 69.06513
## 3.261905 2022-04-02 07:00:00          81.45693 75.57699 87.33687 72.46434
## 3.263393 2022-04-02 08:00:00          79.47121 73.51889 85.42354 70.36792
## 3.264881 2022-04-02 09:00:00          85.98512 79.96128 92.00896 76.77245
## 3.266369 2022-04-02 10:00:00          75.99642 69.90190 82.09094 66.67566
##             Hi_95
## 3.258929 95.71061
## 3.260417 86.82613
## 3.261905 90.44952
## 3.263393 88.57451
## 3.264881 95.19779
## 3.266369 85.31718
plot_ly(Forcasted_Opening, x=~Forcasted_Dates,y=~Forcasted_Opening ,type = 'scatter')

Based on graph above, high potential for the actuator will Fail or will repetitively hit 100% opening at second week of May 2022. So, the management shall arrange the Turn Around maximum in early May 2022 to avoid unplanned trip.

10 Conclusion

  1. The timeseries data is additive type.
  2. Holt’s Winters is better than than STLM to forecast based on MAPE value with monthly as interval.
  3. High potential for the actuator will Fail or will repetitively hit 100% opening at second week of May 2022.
  4. The management shall arrange the Turn Around maximum in early May 2022 to avoid unplanned trip.