Ciência dos Dados - Trabalho de Predição

ENGF08 - TOPICOS ESPECIAIS EM ENGENHARIA I

Universidade Federal da Bahia - UFBA
Escola Politécnica
Programa de Pós-Graduação em Engenharia Industrial - PEI/UFBA
Departamento de Engenharia Química - DEQ
Ciência dos Dados na Engenharia
ENG438 - Tópicos Especiais em Engenharia Química

---

Modelagem - Sistema de Lagoas Aeradas de uma ETE em Indústria de Papel e Celulose

Introdução ao problema: A P & C tem supervisionado e controlado seu sistema biológico de tratamento (sistema de lagoas aeradas) através da demanda bioquímica de oxigênio (DBO), demanda química de oxigênio (DQO), pH, sólidos em suspensão (SS), nitrogênio nitrato (NN), nitrogênio amoniacal (NAm), fósforo (P), cor (Col), temperatura (T), condutividade (Cond) e vazão do efluente (FR) na entrada (in) e saída (out) do sistema de lagoas aeradas. Eles também dispõem de dados diários de precipitação (RF), produção de papel (Pap) e produção de celulose (Pulp). A DBO de 5 dias é um dos parâmetros utilizados na monitorização da qualidade do efluente. Medidas de controle da qualidade do efluente tratado são dificultadas tendo em vista o tempo de 5 dias requerido na análise laboratorial.

* Perguntas a serem respondidas no trabalho:

  + É possível propor um modelo para predição da DBO, DQO ou outro parâmetro de qualidade que justifique a modelagem?
  + A análise exploratória, temporal, agrupamentos, correlações, distribuições, outliers, etc. trazem informações relevantes para o problema?
  + Que outras perguntas poderiam nortear novos estudos?

---

Modelos propostos

Para a concepção do trabalho foram criados modelos de aprendizado de máquina no software R a partir do dataset disponibilizado.
No primeiro momento para o efluente de entrada no sistema de lagoas aeradas, foi realizado análise de correlação e análise de clusters tendo em vista melhor compreensão do fenômenos em questão e avaliação de possíveis outliers.
Buscando um possível melhor modelo e com menor quantidade de variáveis independentes, foi analisado qual o principal parâmetro de qualidade da água a ser analisado e predito na entrada do sistemas de lagoas. E, o parâmetro de maior destaque é a DBO, pois para analisar este parâmetro é necessário pelo menos 5 dias de “espera” para obtenção de resultados laboratoriais do valor do parâmetro com baixa incerteza, que é a DBO5.
Quanto a possível análise de para os parâmetros de qualidade da água do efluente a ser lançado no rio Mogi-Guaçu, presou-se também pela estimativa da DBO em função do lag para análise e de ser o único parâmetro com informação disponível que é necessário para avaliação do lançamento do efluente industrial no corpo receptor, enquadrado como águas doces de classe 2, segundo a resolução CONAMA 430, tendo como exigência a remoção mínima de 60% de DBO.
Os algoritmos de aprendizado de máquina testados para o efluente a ser tratado foram: Random Forest (rf), Ridge/Lasso/Elastic, Linear Model (lm), Modelos Lineares Generalizados (glm), Support Vector Machines (SVM) e eXtreme Gradient Boosting (xgboost). Entretanto, apenas o código do principal algorítimo/modelo escolhido será apresentado neste documento, e ressaltamos que xgboost e svm não foram totalmente concluídos pois ao longo das análises percebemos que há um elevado custo computacional associado a esses algoritmos (resultando em um problema de escalabilidade) e nos teste com dados re-amostrados não percebemos um ganho satisfatório no possível resultado final.
Para o modelo de DBO de saída, foi considerado como input dados preditos do modelo escolhido no modelo de DBO de entrada, tendo em vista a possibilidade de ajuste do efluente durante o processo de tratamento dado um cenário totalmente desfavorável, com lançamento de efluentes fora do range de adequabilidade.
Para o modelo preditor, tendo em vista a escolha do modelo de DBO de saída, estabeleceu-se um critério de classificação, se atende ou não atende a CONAMA 430. Esta proposta também orienta para possíveis problemas relacionados ao número de amostras, n de cada classe (lançamento de efluente adequado ou não).

Biblioteca Utilizada

Para o desenvolvimento dos modelos, foram utilizadas as bibliotecas (packages) descritas no chunk abaixo. Vale o destaque para as bibliotecas de manipulação de dados (tidyverse) e de modelagem (caret).

if(!require("tidyverse")) install.packages("tidyverse") ; library(tidyverse)
if(!require("GGally")) install.packages("GGally") ; library(GGally) 
if(!require("rgl")) install.packages("rgl") ; library(rgl) 
if(!require("car")) install.packages("car") ; library(car)
if(!require("knitr")) install.packages("knitr") ; library(knitr) 
if(!require("reshape2")) install.packages("reshape2") ; library(reshape2)
if(!require("caret")) install.packages("caret") ; library(caret)
if(!require("randomForest")) install.packages("randomForest") ; library(randomForest)
if(!require("readxl")) install.packages("readxl") ; library(readxl)
if(!require("factoextra")) install.packages("factoextra") ; library(factoextra)
if(!require("DMwR")) install.packages("DMwR") ; library(DMwR) 
if(!require("DT")) install.packages("DT") ; library(DT) 
if(!require("ggpmisc")) install.packages("ggpmisc") ; library(ggpmisc) 
if(!require("plotly")) install.packages("plotly") ; library(plotly)
if(!require("shiny")) install.packages("shiny") ; library(shiny)

Dataset e segmentação

O dataset contendo as informações do trabalho foi importado e segmentado em duas parte, entrada e saída. Tal decisão foi realizada para facilitar futuras manipulações de dados considerando as duas vertentes dos modelos preditores, DBO de entrada e da saída das lagoas.

# Caminho da pasta de trabalho
# setwd("D:/PEI/Ciencia_dos_Dados/PREDICAO_TRAB/Trabalho_DQO_DBO")

# Import dataset - Arquivo modificado para MS Excel 2013
Dados_Brutos <- readxl::read_excel("D:/PEI/Ciencia_dos_Dados/PREDICAO_TRAB/Trabalho_DQO_DBO/Dados Brutos.xlsx")

head(Dados_Brutos, 5)

## # A tibble: 5 x 24
##   Datein              `FR(m3/dia)` `BODin (ppm)` `CODin (ppm)` `SSin(ppm)`
##   <dttm>                     <dbl>         <dbl>         <dbl>       <dbl>
## 1 1996-09-01 00:00:00        79757           236           465          NA
## 2 1996-09-02 00:00:00        78818           270           483          NA
## 3 1996-09-03 00:00:00        82039           292           571          NA
## 4 1996-09-04 00:00:00        79501           264           461          NA
## 5 1996-09-05 00:00:00        82461           299           471         181
## # ... with 19 more variables: pHin <dbl>, `NAmin(ppm)` <dbl>, `NNin
## #   (ppm)` <dbl>, `Pin (ppm)` <dbl>, `Colin(ppm)` <dbl>, `Tin(°C)` <dbl>,
## #   `Condin L1` <dbl>, `RF(mm)` <dbl>, `Pulp(ton/dia)` <dbl>, `Pap
## #   (ton/dia)` <dbl>, X__1 <lgl>, Dateout <dttm>, `BODout(ppm)` <dbl>,
## #   `CODout (ppm)` <dbl>, `SSout(ppm)` <dbl>, `Colout(ppm)` <dbl>,
## #   `Tout(°C)` <dbl>, Condout <dbl>, `FRout(m3/dia)` <dbl>

#Segmentação do dataset
entrada <- Dados_Brutos[,1:15]

saida   <- Dados_Brutos[,17:24]

# Plot de ts da DBO de entrada
shiny::div(plotly::ggplotly(ggplot2::ggplot(entrada) +
                              geom_line(aes(x  = Datein,
                                            y  = `BODin (ppm)`)) +
                              labs(title = "DBO in",
                                   x = "Data",
                                   y = "ppm") +
                              theme_bw()), align = "center")

# Plot de ts da DBO de saída
shiny::div(plotly::ggplotly(ggplot2::ggplot(saida) +
                              geom_line(aes(x  = Dateout,
                                            y  = `BODout(ppm)`)) +
                              labs(title = "DBO out",
                                   x = "Data",
                                   y = "ppm") +
                              theme_bw()), align = "center")

Para uma compreensão inicial sobre a DBO de entrada na ETE, realizou-se plot da série temporal (ts) do parâmetro em questão.
A ts de entrada apresenta ser estacionário, porém a mesma tem maiores variações para os menores valor de DBO em jun/200. Aparentemente a ts da DBO de saída apresenta também ser estacionária, entretanto de set/1996 à aproximadamente fev/1997 e no fim do ano 2000 a série plotada demonstra indícios de médias e variância distintas quando visualizado o período de fev/1997 à aproximadamente jun/2000.

Manipulação de Dados

Notou-se que o presente dataset contém muitos valores faltantes (NA), e este foi o primeiro fator que foi escolhido como critério de redução de variáveis, onde variáveis com muitos NA's foram excluídas. Tal tomada de decisão foi escolhida também de acordo com a expertise sobre o fenômeno em questão. Sabe-se que a DBO pode ser parcialmente analisada e correlacionada por parâmetros de qualidade da água, como: DQO, cor total, cor aparente, pH, sólidos dissolvidos e sólidos em suspensão.
Espera-se também que com este critério de remoção de variáveis com muitos NA, os modelos aqui propostos possam ser “implementados”, pois a ausência da medição dos parâmetros de qualidade da água no dia-a-dia da empresa pode ser uma entrave para aplicação real do modelo.

# Check de NA
df <- as.data.frame(unclass(summary(entrada)),
                    check.names      = FALSE, 
                    stringsAsFactors = FALSE)

df <- df[7,1:15]

# View(df)
head(df)

##            Datein    FR(m3/dia)   BODin (ppm)   CODin (ppm)     SSin(ppm)
## X.6 NA's   :107   NA's   :107   NA's   :196   NA's   :196   NA's   :969  
##              pHin    NAmin(ppm)     NNin (ppm)      Pin (ppm)
## X.6 NA's   :160   NA's   :877   NA's   :1258   NA's   :1277  
##        Colin(ppm)       Tin(°C)     Condin L1        RF(mm) Pulp(ton/dia)
## X.6 NA's   :158   NA's   :574   NA's   :163   NA's   :362   NA's   :210  
##     Pap (ton/dia)
## X.6 NA's   :199

Para o modelo de saída da ETE, foi considerado apenas o lag de 02 dias no sistema, tempo este relacionado ao tempo de detenção hidráulica (TDH).

# Removendo variáveis com demasiada quantidade de NA's
df <- entrada %>% 
  dplyr::select(c('Datein', 'FR(m3/dia)', 'BODin (ppm)', 
                  'CODin (ppm)', 'pHin', 'Colin(ppm)', 
                  'Condin L1', 'Pulp(ton/dia)', 'Pap (ton/dia)')) %>% 
  as.data.frame() %>% 
  dplyr::mutate('Pulp(ton/dia)' = ifelse(`Pulp(ton/dia)` < 0, NA, `Pulp(ton/dia)`)) %>% 
  na.omit()
  # Modificação de dados aberrantes de produção de polpa (dados negativos)

# View do dataframe
DT::datatable(df, options = list(
  initComplete = JS(
    "function(settings, json) {",
    "$(this.api().table().header()).css({'font-size': '9px', 'background-color': '#c2d1f0', 'color': '#000'});",
    "}"))) %>%  
  formatStyle(columns = colnames(.$x$data), 
              `font-size` = '9px')

Análise descritiva

Para verificar visualmente possível correlação e direcionar a interpretação inicial do possível comportamento do fenômenos, foi gerado um quadro com informações descritivas básicas e plot da matriz de correlação, sendo apresentada na figura abaixo.

# Para uso futuro de DBO de saída
df_data <- df

df <- df %>% 
  dplyr::select(-Datein)

summary(df)

##    FR(m3/dia)     BODin (ppm)     CODin (ppm)         pHin       
##  Min.   : 5913   Min.   : 41.0   Min.   :136.0   Min.   : 0.850  
##  1st Qu.:61751   1st Qu.:222.0   1st Qu.:515.0   1st Qu.: 6.840  
##  Median :65959   Median :247.0   Median :572.0   Median : 7.110  
##  Mean   :67024   Mean   :247.9   Mean   :572.5   Mean   : 7.448  
##  3rd Qu.:75836   3rd Qu.:273.0   3rd Qu.:626.0   3rd Qu.: 7.480  
##  Max.   :97850   Max.   :449.0   Max.   :925.0   Max.   :12.530  
##    Colin(ppm)       Condin L1    Pulp(ton/dia)    Pap (ton/dia)   
##  Min.   :  41.0   Min.   : 771   Min.   :   0.0   Min.   : 382.4  
##  1st Qu.: 400.0   1st Qu.:1314   1st Qu.: 861.2   1st Qu.:1003.0  
##  Median : 479.0   Median :1467   Median : 918.6   Median :1053.6  
##  Mean   : 465.7   Mean   :1538   Mean   : 886.9   Mean   :1040.7  
##  3rd Qu.: 532.0   3rd Qu.:1685   3rd Qu.: 962.6   3rd Qu.:1098.8  
##  Max.   :1151.0   Max.   :4400   Max.   :1108.4   Max.   :1304.8

GGally::ggpairs(df) +
  labs(title = "Correlação, dispersão e distribuição de dados de entrada na ETE") + 
  theme(axis.text.x  = element_text(size = 5)) +
  theme_bw(base_size = 8)

Observa-se que os parâmetros DQO e cor são os mais bem correlacionados com a DBO de entrada. Os pontos dispersos nos plots indicam uma possível dificuldade de estabelecimento de modelos bem ajustados, com relativa variância e baixo bias.

Modelagem - Modelo geral para DBO in

Buscando um bom ajuste e para o modelo de predição de DBO na entrada da lagoa, realizou-se uso do algoritmo rf considerando o data frame df, sendo a variável dependente a DBO e como variável independente todos os outros restantes parâmetros de qualidade da água.
Nos modelos de rf aqui desenvolvidos para regressão de cada cluster ou modelo geral adotou-se o pré-processamento de dados baseado apenas em z-score (scale()) e centralização dos dados (center()), o treinamento é controlado por validação cruzada com k = 10, e o número de árvores igual a 1000 (um mil), como sugerido por Breiman (2002).
No chunck abaixo é desenvolvido o modelo de regressão para DBO de entrada utilizando os seguintes parâmetros de qualidade da água do efluente de: DQO [ppm], pH, cor [ppm], condutividade elétrica, produção de polpa [ton/dia] e produção de papel [ton/dia].

dft <- df

{
  set.seed(1)
  
  intrain  <- createDataPartition(y = dft$`BODin (ppm)`, p = 0.75, list = FALSE)
  training <- dft[intrain,]
  testing  <- dft[-intrain,]
  
  control  <- caret::trainControl(method  = "repeatedcv", 
                                  number  = 10, 
                                  repeats = 10,
                                  savePredictions = T)
  
  rf_g <- caret::train(`BODin (ppm)` ~ ., 
                       data   = training, 
                       method = "rf",  
                       preProcess = c("center","scale"),
                       ntree = 1000, 
                       importance = TRUE, 
                       verbose = TRUE,
                       metric  = "RMSE",
                       trControl  = control)
  
  print(rf_g)
  rf_g$results
  }

## Random Forest 
## 
## 922 samples
##   7 predictor
## 
## Pre-processing: centered (7), scaled (7) 
## Resampling: Cross-Validated (10 fold, repeated 10 times) 
## Summary of sample sizes: 830, 829, 830, 832, 830, 830, ... 
## Resampling results across tuning parameters:
## 
##   mtry  RMSE      Rsquared   MAE     
##   2     32.82906  0.4725390  24.92251
##   4     32.35651  0.4799343  24.36392
##   7     32.21932  0.4822038  24.07555
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 7.

##   mtry     RMSE  Rsquared      MAE   RMSESD RsquaredSD    MAESD
## 1    2 32.82906 0.4725390 24.92251 3.065406  0.1035840 1.606316
## 2    4 32.35651 0.4799343 24.36392 3.109608  0.1030105 1.606265
## 3    7 32.21932 0.4822038 24.07555 3.069469  0.1046997 1.600809

x <- varImp(rf_g)

x <- as.data.frame(x[["importance"]])

x <- data.frame(Parametro = rownames(x), Importancia = x$Overall)

# levels(x$Parametro)

levels(x$Parametro) <- c("DQO", "Cor", "Condutividade Elétrica",
                         "Vazão de Entrada", "Produçãp de Papel", " Produção de Polpa",
                         "pH")

# Plot de importância de variáveis
ggplot2::ggplot(x) +
  geom_bar(aes(x = Parametro,
               y = Importancia,
               fill = Parametro),
           col  = "black",
           stat = "identity") +
  labs(title = "Importância das variáveis - Modelo geral",
       x = "Parâmetros Análisados", 
       y = "Scores de Importância", 
       color = "Legenda") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90))

final <- data.frame(Real    = testing$`BODin (ppm)`,
                    Pred_rf = predict(rf_g, newdata = testing))

final <- final %>%
  dplyr::mutate(ord = seq(1, nrow(final),1)) %>% 
  reshape2::melt(id.vars = c("ord"))
  
# Plot de ajuste ao longo do tempo
shiny::div(plotly::ggplotly(ggplot2::ggplot(final) + 
                              geom_line(aes(x   = ord,
                                            y   = value, 
                                            col = variable),
                                        size  = 1,
                                        alpha = 0.4) +
                              scale_color_discrete(labels = c("Observado (treinamento)",
                                                              "Predito")) +
                              labs(title = "DBOin [ppm] - Modelo geral",
                                   x = "Período", 
                                   y = "DBO [ppm]", 
                                   color = "Legenda") +
                              theme_bw()), align = "center")

dft$pred <- predict(rf_g, newdata = dft)

dftx <- dft %>% 
  dplyr::select(pred, `BODin (ppm)`, `CODin (ppm)`) %>% 
  dplyr::rename("DQO" = `CODin (ppm)`) %>% 
  reshape2::melt(id.vars = c("DQO"))

ggplot2::ggplot(dftx) +
  geom_point(aes(x = DQO,
                 y = value,
                 col = variable),
             alpha = 0.4) +
  scale_color_discrete(labels = c("Observado (treinamento)", "Predito")) +
  labs(title = "DBOin vs. DQO - Modelo geral",
       x = "DQO [ppm]", 
       y = "DBO [ppm]", 
       color = "Legenda DBO") +
  theme_bw()

dftx <- dft %>% 
  dplyr::select(pred, `BODin (ppm)`, `Colin(ppm)`) %>% 
  dplyr::rename("Cor" = `Colin(ppm)`) %>% 
  reshape2::melt(id.vars = c("Cor"))

ggplot2::ggplot(dftx) +
  geom_point(aes(x = Cor,
                 y = value,
                 col = variable),
             alpha = 0.4) +
  scale_color_discrete(labels = c("Observado (treinamento)", "Predito")) +
  labs(title = "DBO vs. Cor - Modelo geral",
       x = "Cor [ppm]", 
       y = "DBO [ppm]", 
       color = "Legenda DBO") +
  theme_bw()

Este modelo de rf generalista (com todas as variáveis já filtradas) apresenta um bom ajuste, tendo como melhor resultado o modelo com R² = 0,48, RMSE = 32,22 e mtry = 7. Este modelo preditor segue o mesmo comportamento dos valores de observados para DBO.
As mais importantes variáveis independentes para este modelo são: DQO in, cor in e vazão de efluente para o sistema de lagoas.

Modelagem - Clustering e outliers

Verificado o possível valor baixo de R², O primeiro passo para desenvolvimento de uma nova proposta de regressão foi o emprego da clusterização dos objetos utilizando o algoritmo kmeans crisp. Onde foi considerado a dissimilaridade a partir da distância Euclidiana.
Do dataset em análise, cada dia é interpretado como um objeto, e para cada relevante agrupamento de objetos foi testado modelos de regressão em rf para checar uma possível melhor ajuste para o modelo.

dfts <- scale(df)

{
  set.seed(1)
  
  factoextra::fviz_nbclust(dfts, kmeans, method = "silhouette") +
  labs(title = "Gráfico de Silhueta",
       x = "Número de k clusters", 
       y = "Média interna da silhueta")
  }

{
  set.seed(1)
  
  km.res <- kmeans(scale(df), 3, nstart = 25)
  }

factoextra::fviz_cluster(km.res, data = dfts,
             ggtheme = theme_bw(),
             main = "Clustering")

O número ótimo sugerido de clusters é 3 (três), como pode ser visualizado no gráfico de silhueta. Entretanto, percebe-se que o cluster 1, visualmente, não apresenta boa coesão entre os objetos, apesar das outras dimensões. Deste modo, foi realizado análise para detecção de outliers utilizando o algorítimo LOF (Local Outlier Factor) desenvolvido por Breunig et al., (2000), que identifica outliers locais baseados em densidade.

dft <- df

outlier.scores <- DMwR::lofactor(dft, k = 5)

# plot(density(outlier.scores))

dfo <- boxplot(outlier.scores, xlab = "Fator outlier local (LOF)")

dfo <- dfo$out

dfo <- data.frame(dfo, dfo_2 = rep(1))

outlier.scores <- outlier.scores %>% 
  as.data.frame() %>% 
  dplyr::mutate(ordem = seq(1, nrow(dft), 1)) %>% 
  dplyr::left_join(dfo, by = c("." = "dfo")) %>% 
  replace(is.na(.), 0)

dft <- dft %>% 
  dplyr::mutate(outlier = outlier.scores$dfo_2) %>% 
  dplyr::filter(outlier == 0) %>% 
  dplyr::select(-outlier)

Feito a remoção dos outliers analisando a distribuição de dados referente ao LOF, foi promovida uma nova clusterização, onde foram estabelecidos apenas 2 clusters para a modelagem.

dfx <- scale(dft)

{
  set.seed(1)
  
  factoextra::fviz_nbclust(dfx, kmeans, method = "silhouette") +
  labs(title = "Gráfico de Silhueta",
       x = "Número de k clusters", 
       y = "Média interna da silhueta")
  }

{
  set.seed(1)
  
  km.res <- kmeans(dfx, 2, nstart = 25)
  }

factoextra::fviz_cluster(km.res, data = dfx,
             ggtheme = theme_bw(),
             main = "Clustering")

dft$cluster <- km.res[["cluster"]]

paste0("Quantidade de objetos por cluster:",
       table(dft$cluster))

## [1] "Quantidade de objetos por cluster:881"
## [2] "Quantidade de objetos por cluster:214"

knit_hooks$set(webgl = hook_webgl)

shiny::div(car::scatter3d(x = dft$`BODin (ppm)`, 
          y = dft$`CODin (ppm)`, 
          z = dft$`Colin(ppm)`, 
          groups = as.factor(dft$cluster),
          data = dft, 
          grid = FALSE, fit = "smooth"), align = "center")

You must enable Javascript to view this page properly.

A remoção de outliers com LOF removeu alguns dos objetos que estavam mais distantes das maiores aglomerações. Entretanto os objetos que continuam causando essa espécia de desbalanceamento espacial, aparentemente aumentam a desordem, entropia, dos clusters gerados.

Modelagem - Random Forest para clusters

Para os dois clusters gerados, foram desenvolvidos modelo de regressão em rf. Nos chunck abaixo cada cluster é modelado independente do outro, sendo utilizado o cluster 1 para o primeiro modelo (rf_c1) e o cluster 2 para o segundo modelo (rf_c2).

dft1 <- dft %>%
  dplyr::filter(cluster == 1) %>% 
  dplyr::select(-cluster)

{
  set.seed(1)
  
  intrain  <- createDataPartition(y = dft1$`BODin (ppm)`, p = 0.75, list = FALSE)
  training <- dft1[intrain,]
  testing  <- dft1[-intrain,]
  
  control <- caret::trainControl(method   = "repeatedcv", 
                                 number   = 10, 
                                 repeats  = 10,
                                 savePredictions = T)
  
  rf_c1 <- caret::train(`BODin (ppm)` ~ ., 
                        data   = training, 
                        method = "rf",  
                        preProcess = c("center","scale"),
                        ntree = 1000, 
                        importance = TRUE, 
                        verbose = TRUE,
                        metric  = "RMSE",
                        trControl  = control)
  
  print(rf_c1)
  rf_c1$results
  }

## Random Forest 
## 
## 662 samples
##   7 predictor
## 
## Pre-processing: centered (7), scaled (7) 
## Resampling: Cross-Validated (10 fold, repeated 10 times) 
## Summary of sample sizes: 596, 596, 596, 598, 595, 595, ... 
## Resampling results across tuning parameters:
## 
##   mtry  RMSE      Rsquared   MAE     
##   2     30.75889  0.3209368  24.12479
##   4     30.81573  0.3160128  24.13796
##   7     30.90288  0.3145542  24.24404
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 2.

##   mtry     RMSE  Rsquared      MAE   RMSESD RsquaredSD    MAESD
## 1    2 30.75889 0.3209368 24.12479 2.842059  0.1300745 2.086871
## 2    4 30.81573 0.3160128 24.13796 2.951500  0.1265886 2.217602
## 3    7 30.90288 0.3145542 24.24404 3.029293  0.1240928 2.281119

Apesar do razoável valor de RMSE encontrado, o máximo valor do R² foi muito baixo (0,316 com RMSE = 30,8). Para possível aprimoramento foi feito análise visual dos plots dos clusters, e sabendo da sensibilidade do algorítimo de rf para outliers, realizou-se nova clusterização objetivada em segregar os dados mais afastados, que tendenciam o modelo para uma situação com péssimo ajuste, underfitting.

dfx <- dft %>% 
  dplyr::select(-cluster) %>% 
  scale()

{
  set.seed(1)
  
  km.res <- kmeans(scale(dft), 3, nstart = 25)
  }

factoextra::fviz_cluster(km.res, data = dfx,
             ggtheme = theme_bw(),
             main = "Clustering")

dft$cluster <- km.res[["cluster"]]

rm(dfo, dfts, outlier.scores, dfx)

paste0("Quantidade de objetos por cluster:",
       table(dft$cluster))

## [1] "Quantidade de objetos por cluster:208"
## [2] "Quantidade de objetos por cluster:870"
## [3] "Quantidade de objetos por cluster:17"

knit_hooks$set(webgl = hook_webgl)

shiny::div(car::scatter3d(x = dft$`BODin (ppm)`, 
          y = dft$`CODin (ppm)`, 
          z = dft$`Colin(ppm)`, 
          groups = as.factor(dft$cluster),
          data = dft, 
          grid = FALSE, fit = "smooth"), align = "center")

You must enable Javascript to view this page properly.

Claramente percebe-se que os objetos mais afastados podem ser um clustering a parte. Entretanto, como o número de objetos no cluster 3 é pequeno, logo este cluster não será considerado na modelagem, os mesmos serão excluídos como outliers.

dft1 <- dft %>%
  dplyr::filter(cluster == 1) %>% 
  dplyr::select(-cluster)

{
  set.seed(1)
  
  intrain  <- createDataPartition(y = dft1$`BODin (ppm)`, p = 0.75, list = FALSE)
  training <- dft1[intrain,]
  testing  <- dft1[-intrain,]
  
  control <- caret::trainControl(method   = "repeatedcv",
                                 number   = 10, 
                                 repeats  = 10,
                                 savePredictions = T)
  
  rf_c1 <- caret::train(`BODin (ppm)` ~ ., 
                        data   = training, 
                        method = "rf",  
                        preProcess = c("center","scale"),
                        ntree = 1000, 
                        importance = TRUE, 
                        verbose = TRUE,
                        metric  = "RMSE",
                        trControl  = control)
  
  print(rf_c1)
  rf_c1$results
  }

## Random Forest 
## 
## 157 samples
##   7 predictor
## 
## Pre-processing: centered (7), scaled (7) 
## Resampling: Cross-Validated (10 fold, repeated 10 times) 
## Summary of sample sizes: 142, 142, 141, 141, 141, 142, ... 
## Resampling results across tuning parameters:
## 
##   mtry  RMSE      Rsquared   MAE     
##   2     26.20264  0.6159159  19.86084
##   4     23.93309  0.6576240  17.20532
##   7     22.81637  0.6761823  15.24148
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 7.

##   mtry     RMSE  Rsquared      MAE   RMSESD RsquaredSD    MAESD
## 1    2 26.20264 0.6159159 19.86084 5.638319  0.1653925 4.197645
## 2    4 23.93309 0.6576240 17.20532 6.123401  0.1670130 4.223911
## 3    7 22.81637 0.6761823 15.24148 6.683228  0.1710768 4.308869

x <- varImp(rf_c1)

x <- as.data.frame(x[["importance"]])

x <- data.frame(Parametro = rownames(x), Importancia = x$Overall)

# levels(x$Parametro)

levels(x$Parametro) <- c("DQO", "Cor", "Condutividade Elétrica",
                         "Vazão de Entrada", "Produçãp de Papel", " Produção de Polpa",
                         "pH")

# Plot de importância de variáveis
ggplot2::ggplot(x) +
  geom_bar(aes(x = Parametro,
               y = Importancia,
               fill = Parametro),
           col  = "black",
           stat = "identity") +
  labs(title = "Importância das variáveis - cluster 1",
       x = "Parâmetros Análisados", 
       y = "Scores de Importância", 
       color = "Legenda") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90))

final <- data.frame(Real    = testing$`BODin (ppm)`,
                    Pred_rf = predict(rf_c1, newdata = testing))

final <- final %>%
  dplyr::mutate(ord = seq(1, nrow(final),1)) %>% 
  reshape2::melt(id.vars = c("ord"))
  
# Plot de ajuste ao longo do tempo
shiny::div(plotly::ggplotly(ggplot2::ggplot(final) + 
  geom_line(aes(x   = ord,
                y   = value, 
                col = variable),
            size  = 1,
            alpha = 0.4) +
  scale_color_discrete(labels = c("Observado (treinamento)", "Predito")) +
  labs(title = "DBOin [ppm] - cluster 1",
       x = "Período", 
       y = "DBO [ppm]", 
       color = "Legenda") +
  theme_bw()), align = "center")

dft1$pred <- predict(rf_c1, newdata = dft1)

dft1x <- dft1 %>% 
  dplyr::select(pred, `BODin (ppm)`, `CODin (ppm)`) %>% 
  dplyr::rename("DQO" = `CODin (ppm)`) %>% 
  reshape2::melt(id.vars = c("DQO"))

ggplot2::ggplot(dft1x) +
  geom_point(aes(x = DQO,
                 y = value,
                 col = variable),
             alpha = 0.4) +
  scale_color_discrete(labels = c("Observado (treinamento)", "Predito")) +
  labs(title = "DBOin vs. DQO - cluster 1",
       x = "DQO [ppm]", 
       y = "DBO [ppm]", 
       color = "Legenda DBO") +
  theme_bw()

dft1x <- dft1 %>% 
  dplyr::select(pred, `BODin (ppm)`, `Colin(ppm)`) %>% 
  dplyr::rename("Cor" = `Colin(ppm)`) %>% 
  reshape2::melt(id.vars = c("Cor"))

ggplot2::ggplot(dft1x) +
  geom_point(aes(x = Cor,
                 y = value,
                 col = variable),
             alpha = 0.4) +
  scale_color_discrete(labels = c("Observado (treinamento)", "Predito")) +
  labs(title = "DBO vs. Cor - cluster 1",
       x = "DQO [ppm]", 
       y = "Cor [ppm]", 
       color = "Legenda DBO") +
  theme_bw()

Como pode ser visualizado nos plots de valores preditos versus valores observados, o modelo proposto para o cluster 1 de DBO de entrada apresenta um primoroso ajuste ainda que dados como do parâmetro cor apresentem dois padrões de distribuição. E os parâmetros de qualidade da água mais importantes para este modelo são a cor, DQO, pH e vazão de entrada na lagoa.

Este modelo proposto apresenta um R² = 0,676 e RMSE = 22,8 para o modelo de rf com mtry = 7. Esquerre (2004), utilizando modelo de regressão linear múltipla e aplicações de análise de componente principal e mínimos quadrados parciais, encontrou um modelo com os parâmetros de DQO, vazão de entrada, temperatura, pH e produção de papel com R² = 0,458.
Para o cluster 2, foi desenvolvido modelo de rf seguindo a mesma metodologia aqui proposta para o cluster 1.

dft2 <- dft %>%
  dplyr::filter(cluster == 2) %>% 
  dplyr::select(-cluster)

{
  set.seed(1)

  intrain  <- createDataPartition(y = dft2$`BODin (ppm)`, p = 0.75, list = FALSE)
  training <- dft2[intrain,]
  testing  <- dft2[-intrain,]

  control <- caret::trainControl(method   = "repeatedcv",
                                 number   = 10,
                                 repeats  = 10,
                                 savePredictions = T)

  rf_c2 <- caret::train(`BODin (ppm)` ~ .,
                        data   = training,
                        method = "rf",
                        preProcess = c("center", "scale"),
                        ntree = 1000,
                        importance = TRUE,
                        verbose = TRUE,
                        metric  = "RMSE",
                        trControl  = control)

  print(rf_c2)
  rf_c2$results
  }

## Random Forest 
## 
## 654 samples
##   7 predictor
## 
## Pre-processing: centered (7), scaled (7) 
## Resampling: Cross-Validated (10 fold, repeated 10 times) 
## Summary of sample sizes: 588, 589, 588, 588, 589, 589, ... 
## Resampling results across tuning parameters:
## 
##   mtry  RMSE      Rsquared   MAE     
##   2     30.42500  0.2244020  23.68731
##   4     30.28248  0.2283126  23.58117
##   7     30.38385  0.2254397  23.64639
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 4.

##   mtry     RMSE  Rsquared      MAE   RMSESD RsquaredSD    MAESD
## 1    2 30.42500 0.2244020 23.68731 2.327493 0.08258643 1.729179
## 2    4 30.28248 0.2283126 23.58117 2.327719 0.08068875 1.787745
## 3    7 30.38385 0.2254397 23.64639 2.357457 0.07974497 1.828753

x <- varImp(rf_c2)

x <- as.data.frame(x[["importance"]])

x <- data.frame(Parametro = rownames(x), Importancia = x$Overall)

# levels(x$Parametro)

levels(x$Parametro) <- c("DQO", "Cor", "Condutividade Elétrica",
                         "Vazão de Entrada", "Produçãp de Papel", " Produção de Polpa",
                         "pH")

# Plot de importância de variáveis
ggplot2::ggplot(x) +
  geom_bar(aes(x = Parametro,
               y = Importancia,
               fill = Parametro),
           col  = "black",
           stat = "identity") +
  labs(title = "Importância das variáveis - cluster 1",
       x = "Parâmetros Análisados", 
       y = "Scores de Importância", 
       color = "Legenda") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90))

Apesar de acompanhar as tendências de comportamento, o presente modelo proposto de rf para o cluster 2 não apresenta um bom ajuste. Logo, para teste de um possível melhor modelo, foram removidas da análise as variáveis independentes de menor importância e ou correlação para que possa ser gerado um novo modelo.

dft2 <- dft %>%
  dplyr::filter(cluster == 2) %>% 
  dplyr::select(-cluster,
                -`Pulp(ton/dia)`,
                -`Pap (ton/dia)`)

GGally::ggpairs(dft2) +
  labs(title = "Correlação, dispersão e distribuição de dados de entrada na ETE") + 
  theme(axis.text.x = element_text(size = 5)) +
  theme_bw(base_size = 8)

{
 set.seed(1)

  intrain  <- createDataPartition(y = dft2$`BODin (ppm)`, p = 0.75, list = FALSE)
  training <- dft2[intrain,]
  testing  <- dft2[-intrain,]

  control  <- caret::trainControl(method   = "repeatedcv",
                                  number   = 10,
                                  repeats  = 10,
                                  savePredictions = T)

   rf_c2 <- caret::train(`BODin (ppm)` ~ .,
                         data   = training,
                         method = "rf",
                         preProcess = c("center", "scale"),
                         ntree = 1000,
                         importance = TRUE,
                         verbose = TRUE,
                         metric  = "RMSE",
                         trControl  = control)

   print(rf_c2)
   rf_c2$results

    }

## Random Forest 
## 
## 654 samples
##   5 predictor
## 
## Pre-processing: centered (5), scaled (5) 
## Resampling: Cross-Validated (10 fold, repeated 10 times) 
## Summary of sample sizes: 588, 589, 588, 588, 589, 589, ... 
## Resampling results across tuning parameters:
## 
##   mtry  RMSE      Rsquared   MAE     
##   2     31.03505  0.1892929  24.40908
##   3     31.11182  0.1888816  24.48453
##   5     31.32499  0.1853248  24.63207
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 2.

##   mtry     RMSE  Rsquared      MAE   RMSESD RsquaredSD    MAESD
## 1    2 31.03505 0.1892929 24.40908 2.260341 0.07704401 1.788039
## 2    3 31.11182 0.1888816 24.48453 2.233202 0.07662852 1.773789
## 3    5 31.32499 0.1853248 24.63207 2.262131 0.07713471 1.788345

final_imp <- data.frame(Real = testing$`BODin (ppm)`,
                        Pred_rf_imp = predict(rf_c2, newdata = testing))

final_imp <- final_imp %>%
  cbind(DQO  = testing$`CODin (ppm)`) %>%
  reshape2::melt(id.vars = c("DQO"))

ggplot2::ggplot(final_imp) +
  geom_point(aes(x   = DQO,
                 y   = value,
                col  = variable),
            size  = 3,
            alpha = 0.4) +
  scale_color_discrete(labels = c("Observado (treinamento)", "Predito")) +
  labs(title = "DBO in [ppm] - cluster 2",
       x = "DQO [ppm]",
       y = "DBO [ppm]",
       color = "Legenda") +
  theme_bw()

Mesmo com a redução das variáveis relacionadas a produção, com baixos valores de importância e de menor correlação, o modelo da rf para o cluster 2 não é tão satisfatório (R² = 0.189), menor que o modelo “completo” para o cluster 2. Pode-se perceber nos gráficos de correlação que os preditores da DBO para este agrupamento não apresentam padrão, os pontos observados são totalmente dispersos. Mesmo sabendo da não-linearidade do modelo a ser predito, a nível de teste/didático, um lm também foi realizado, tornando evidente, mais uma vez, a elevada dispersão entre os dados e e ausência de tendências ou padrões para o cluster 2.

my.formula <- y ~ x

ggplot(data = dft2, aes(x = `CODin (ppm)`,
                        y = `BODin (ppm)`)) +
   geom_smooth(method   = "lm", se = T, color = "black",
               formula  = my.formula) +
   stat_poly_eq(formula = my.formula,
                eq.with.lhs = "italic(DBO)~`=`~",
                eq.x.rhs  = "~italic(DQO)",
                aes(label = paste(..eq.label.., ..rr.label.., sep = "*plain(\",\")~")), 
                parse = TRUE) +         
   geom_point() +
   labs(title = "DBO in [ppm] - cluster 2",
       x = "DQO [ppm]",
       y = "DBO [ppm]",
       color = "Legenda") +
   theme_bw()

Deste modo, o modelo escolhido para predição da DBO de entrada na lagoa é o modelo rf geral sem clusters e tratamento de possíveis outliers, sendo utilizado como variáveis independentes: DQO, pH, cor, condutividade, produção de polpa, produção de papel e vazão de entrada.
O modelo preditor rf tem R² = 0,48, RMSE = 32,2, MAE = 24,07 e mtry = 7, onde mtry é o número de variáveis disponíveis para dividir em cada nó da árvore.
Vale a pena ressaltar que o parâmetro temperatura pode ser um bom preditor para DBO, mas a ausência do mesmo no modelo proposto representa um maior ganho de aplicabilidade real, tendo em vista que há um elevado registro de dados ausentes (NA) para esse parâmetro.

Modelagem - DBO out

Sabendo que a análise laboratorial de DBO demora até 5 dias, considerando que a DBO de entrada deve ser a variável mais latente para predição de DBO de saída e que o modelo de DBO de entrada proposto representa muito bem o comportamento do fenômeno com baixo erro associado a predição, e o TDH de dois dias no sistema em análise, pressupõem-se que um modelo de classificação baseado na predição de remoção de 60% de DBO (CONAMA 430) é o ideal para orientar as operações da ETE da empresa.
Logo, será criado um atributo de classificação do sistema baseado no percentual de remoção de DBO que rotulará o modelo de classificação baseado em rf considerando os atributos de maior importância para predição de DBO de entrada, a DBO de entrada. Tal cenário proposto é idealizado para uma possível automação de correção de falhas no sistema, pois o tempo necessário para análises dos parâmetros é maior que o TDH.

# Data frame de info inicial
dft <- df %>% 
  dplyr::mutate(dboin_pred = predict(rf_g, newdata = df),
                Datain = df_data$Datein) %>% 
  dplyr::select(Datain, dboin_pred, `BODin (ppm)`, 
                `CODin (ppm)`, `Colin(ppm)`,
                `FR(m3/dia)`, `Pulp(ton/dia)`, `Condin L1`,
                `Pap (ton/dia)`) %>% 
  dplyr::mutate(dboin_pred = dplyr::lag(dboin_pred, 2L),
                Datain  = dplyr::lag(Datain, 2L),
                Dataout = Datain + lubridate::days(2)) %>% 
  dplyr::rename("Cor_in"   = `Colin(ppm)`,
                "DQO_in"   = `CODin (ppm)`,
                "Vazao_in" = `FR(m3/dia)`,
                "Condutividade" = `Condin L1`,
                "Prod_polpa" = `Pulp(ton/dia)`) %>% 
  na.omit()

# Verificação das correlações de DBO de saída
saida2 <- saida %>% 
  dplyr::select(-Dateout)

# Checando NA's
colSums(is.na(saida2))

##   BODout(ppm)  CODout (ppm)    SSout(ppm)   Colout(ppm)      Tout(°C) 
##           194           192          1345           150           303 
##       Condout FRout(m3/dia) 
##           152           109

saida2 <- saida2 %>% 
  na.omit()

summary(saida2)

##   BODout(ppm)      CODout (ppm)     SSout(ppm)      Colout(ppm)   
##  Min.   : 18.00   Min.   :166.0   Min.   : 10.00   Min.   : 55.0  
##  1st Qu.: 70.75   1st Qu.:278.0   1st Qu.: 39.00   1st Qu.:370.8  
##  Median : 87.00   Median :315.0   Median : 52.00   Median :446.0  
##  Mean   : 87.38   Mean   :317.4   Mean   : 51.72   Mean   :422.2  
##  3rd Qu.:103.25   3rd Qu.:351.2   3rd Qu.: 62.00   3rd Qu.:511.5  
##  Max.   :159.00   Max.   :616.0   Max.   :107.00   Max.   :936.0  
##     Tout(°C)        Condout     FRout(m3/dia)   
##  Min.   :28.00   Min.   : 950   Min.   :  5913  
##  1st Qu.:33.88   1st Qu.:1358   1st Qu.: 62110  
##  Median :35.00   Median :1502   Median : 66257  
##  Mean   :34.92   Mean   :1548   Mean   : 67631  
##  3rd Qu.:37.00   3rd Qu.:1660   3rd Qu.: 76461  
##  Max.   :38.00   Max.   :2800   Max.   :106942

GGally::ggpairs(saida2) +
  labs(title = "Correlação, dispersão e distribuição para modelo de DBO out na ETE") + 
  theme(axis.text.x  = element_text(size = 5)) +
  theme_bw(base_size = 8)

A quantidade de NA do df de saída é bem menor que o df de entrada, e a variável DQO de saída é o parâmetro de qualidade da água que melhor se correlaciona com a DBO de saída e apresenta apenas 192 NA's.
Contudo, nenhum dos parâmetros observados (valores reais) de qualidade da água avaliados do efluente de saída da lagoa será utilizado para a predição da classe de DBO de saída.

# Montando data frame para predição de DBO de saída

saida2 <- dft %>% 
  dplyr::right_join(saida, by = c("Dataout" = "Dateout")) %>% 
  dplyr::select(Datain, Dataout, dboin_pred, `BODout(ppm)`, 
                `BODin (ppm)`, Cor_in, DQO_in, Vazao_in,
                `CODout (ppm)`, Prod_polpa, Condutividade,
                `Pap (ton/dia)`) %>%
  na.omit() %>% 
  dplyr::mutate(remocao_percent = (1-(`BODout(ppm)`/`BODin (ppm)`))*100,
                remocao = ifelse(remocao_percent > 60, "sim", "nao")) %>% 
  dplyr::rename("DBO_in"  = `BODin (ppm)`,
                "DQO_out" = `CODout (ppm)`,
                "DBO_out" = `BODout(ppm)`)

# remocao: % de remocao de DBO na ETE
# Sim == Atende a CONAMA 430 de padrão de lançamento de efluentes;
# Não == Não Atender a CONAMA 430.

DT::datatable(saida2, options = list(
  initComplete = JS(
    "function(settings, json) {",
    "$(this.api().table().header()).css({'font-size': '9px', 'background-color': '#c2d1f0', 'color': '#000'});",
    "}"))) %>%  
  formatStyle(columns = colnames(.$x$data), 
              `font-size` = '9px')

saidax <- saida2 %>% 
  dplyr::select(DBO_in, dboin_pred, DBO_out) %>% 
  dplyr::mutate(ord = seq(1, nrow(saida2), 1)) %>% 
  reshape2::melt(id.vars = c("ord"))

# levels(saidax$variable)
levels(saidax$variable) <- c("DBO in observada",
                             "DBO in Predita", 
                             "DBO out observada")

shiny::div(plotly::ggplotly(ggplot2::ggplot(saidax) +
  geom_line(aes(x   = ord,
                y   = value,
                col = variable),
            group   = 1) +
  facet_grid(~variable) +
  labs(title = " ",
       x = "Observação",
       y = "[ppm]",
       color = "Legenda") +
  theme_bw()), align = "center")

saidax <- saida2 %>%
  dplyr::select(-remocao, -Dataout, -Datain, -remocao_percent)

GGally::ggpairs(saidax) +
  labs(title = "Correlação, dispersão e distribuição para modelo de DBO out na ETE") + 
  theme(axis.text.x  = element_text(size = 5)) +
  theme_bw(base_size = 8)

Como a DBO de entrada predita está bem ajustada e reproduz bem o comportamento e padrão da DBO de entrada observada, e as variáveis aqui utilizadas são de mais fácil análise laboratorial ou mensuração, o modelo de classificação aqui proposto pressupõem uma possível automatização do sistema de avaliação de qualidade do parâmetro DBO, pois o tempo de detenção hidráulica é menor que o tempo necessário para avaliação laboratorial da DBO5. Logo, será predito a classificação da eficiência da ETE para o parâmetro DBO sem a necessidade de emprego do mesmo.
Para tal, ainda avaliando as correlações entre parâmetros de qualidade da água, a DQO de saída será predita para poder incrementar o modelo de DBO de saída.

# saidax

saidax <- saidax %>% 
  dplyr::select(-DBO_out, -DBO_in)

{
  set.seed(1)
  
  intrain  <- createDataPartition(y = saidax$DQO_out, p = 0.75, list = FALSE)
  training <- saidax[intrain,]
  testing  <- saidax[-intrain,]
  
  control  <- caret::trainControl(method  = "repeatedcv", 
                                  number  = 10, 
                                  repeats = 10,
                                  savePredictions = T)
  
  rf_dqo <- caret::train(DQO_out ~ ., 
                         data   = training, 
                         method = "rf",  
                         preProcess = c("center","scale"),
                         ntree = 1000, 
                         importance = TRUE, 
                         verbose = TRUE,
                         metric  = "RMSE",
                         trControl  = control)
  
  print(rf_dqo)
  rf_dqo$results
  }

## Random Forest 
## 
## 872 samples
##   7 predictor
## 
## Pre-processing: centered (7), scaled (7) 
## Resampling: Cross-Validated (10 fold, repeated 10 times) 
## Summary of sample sizes: 785, 786, 786, 784, 785, 783, ... 
## Resampling results across tuning parameters:
## 
##   mtry  RMSE      Rsquared   MAE     
##   2     50.89835  0.4771778  35.26798
##   4     51.30270  0.4695469  35.49111
##   7     52.20407  0.4552929  36.00623
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 2.

##   mtry     RMSE  Rsquared      MAE   RMSESD RsquaredSD    MAESD
## 1    2 50.89835 0.4771778 35.26798 8.892549 0.08284646 3.731508
## 2    4 51.30270 0.4695469 35.49111 9.143224 0.09248246 3.937532
## 3    7 52.20407 0.4552929 36.00623 9.519444 0.09973683 4.122038

xdqo <- varImp(rf_dqo)

xdqo <- as.data.frame(xdqo[["importance"]])

xdqo <- data.frame(Parametro = rownames(xdqo), Importancia = xdqo$Overall)

# levels(xdqo$Parametro)

levels(xdqo$Parametro) <- c("Condutividade Elétrica", " Cor in", "DQO in",
                            "Prod. polpa", "Vazão in", "Prod. Papel","DBO in predita")

# Plot de importância de variáveis
ggplot2::ggplot(xdqo) +
  geom_bar(aes(x = Parametro,
               y = Importancia,
               fill = Parametro),
           col  = "black",
           stat = "identity") +
  labs(title = "Importância das variáveis - Modelo geral",
       x = "Parâmetros Análisados", 
       y = "Scores de Importância", 
       color = "Legenda") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90))

final_dqo <- data.frame(Observado    = testing$DQO_out,
                        Predito_rf = predict(rf_dqo, newdata = testing))

final_dqo <- final_dqo %>%
  dplyr::mutate(ord = seq(1, nrow(final_dqo),1)) %>% 
  reshape2::melt(id.vars = c("ord"))
  
# Plot de ajuste ao longo do tempo
shiny::div(plotly::ggplotly(ggplot2::ggplot(final_dqo) + 
  geom_line(aes(x   = ord,
                y   = value, 
                col = variable),
            size  = 1,
            alpha = 0.4) +
  scale_color_discrete(labels = c("Observado (treinamento)", "Predito")) +
  labs(title = "DQO out [ppm]",
       x = "observação", 
       y = "DQO [ppm]", 
       color = "Legenda") +
  theme_bw()), align = "center")

Apesar do valor de R² não ser muito alto e o RMSE ser relativamente alto, o modelo preditor de DQO de saída apresenta um bom ajuste, principalmente quanto ao comportamento da variação do parâmetro ao longo do tempo. Deste modo, este modelo também servirá como input de informação independente para o modelo de predição da DBO de saída.
Este cenário de predição de valores de DBO de entrada e saída, apesar das incertezas que rondam o modelo aqui proposto, se apresenta como uma boa ferramenta para suporte a tomada de decisão dos gestores da planta industrial, pois existe a possibilidade de propor estratégias de controle e operação das lagoas aeradas a partir do momento que se tem uma estimativa do possível não enquadramento do efluente a ser lançado.

saidax <- saida2 %>% 
  dplyr::mutate(DQO_out_pred = predict(rf_dqo, newdata = saida2)) %>% 
  dplyr::select(-Datain, -Dataout, -DBO_in, -DQO_out, 
                -remocao_percent, -DBO_out) %>% 
  na.omit()

# saidax <- saida2 %>% 
#   dplyr::select(dboin_pred, DQO_in, Cor_in, Vazao_in, remocao)

paste0("Quantidade de objetos por classe de DBO de saída:",
       table(saidax$remocao))

## [1] "Quantidade de objetos por classe de DBO de saída:325"
## [2] "Quantidade de objetos por classe de DBO de saída:836"

table(saidax$remocao)

## 
## nao sim 
## 325 836

Observa-se que o modelo de classificação está desbalanceado, logo será empregado técnicas de resampling para refinar o melhor modelo de classificação pois a maioria dos algoritmos de classificação de aprendizado de máquina é sensível ao desequilíbrio nas classes preditoras.
Além do modelo considerando os dados originais (desbalanceados), os modelos de resampling aqui testados serão: Under-sampling, Oversampling e SMOTE.

# Modelo com dados originais [desbalanceados]

{
  set.seed(1)
  
  index <- createDataPartition(saidax$remocao, p = 0.75, list = FALSE) 
  train <- saidax[index, ]
  test  <- saidax[-index, ]

  ctrl <- caret::trainControl(method  = "repeatedcv", 
                              number  = 10, 
                              repeats = 10,
                              classProbs = T,
                              savePredictions = T)

  model_orig <- caret::train(remocao ~ .,
                             data   = train,
                             method = "rf",
                             preProcess = c("scale", "center"),
                             trControl  = ctrl,
                             ntree  = 1000)
  }

model_orig

## Random Forest 
## 
## 871 samples
##   8 predictor
##   2 classes: 'nao', 'sim' 
## 
## Pre-processing: scaled (8), centered (8) 
## Resampling: Cross-Validated (10 fold, repeated 10 times) 
## Summary of sample sizes: 783, 784, 783, 784, 784, 785, ... 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa    
##   2     0.7517793  0.2741141
##   5     0.7446681  0.2832799
##   8     0.7411005  0.2809963
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 2.

final_o <- data.frame(Real = test$remocao,
                      predict(model_orig, newdata = test, type = "prob"))

final_o$predict <- ifelse(final_o$nao > 0.5, "nao", "sim")

final_o$predict <- as.factor(final_o$predict)
test$remocao    <- as.factor(test$remocao)

# levels(final_o$predict)
# levels(test$remocao)

cm_orig <- confusionMatrix(final_o$predict, test$remocao)

# print(paste("Matriz de confusão do modelo desbalanceado"))
# cm_orig

No chunck abaixo está o modelo com subamostragem (under-sampling).

# Modelo com resampling under-sampling

{
  set.seed(1)
  
  index <- createDataPartition(saidax$remocao, p = 0.75, list = FALSE) 
  train <- saidax[index, ]
  test  <- saidax[-index, ]

  ctrl <- trainControl(method  = "repeatedcv", 
                       number  = 10, 
                       repeats = 10,
                       verboseIter = FALSE,
                       sampling   = "down",
                       classProbs = T, 
                       savePredictions = T)

  model_und <- caret::train(remocao ~ .,
                            data   = train,
                            method = "rf",
                            preProcess = c("scale", "center"),
                            trControl  = ctrl,
                            ntree  = 1000)
  }

model_und

## Random Forest 
## 
## 871 samples
##   8 predictor
##   2 classes: 'nao', 'sim' 
## 
## Pre-processing: scaled (8), centered (8) 
## Resampling: Cross-Validated (10 fold, repeated 10 times) 
## Summary of sample sizes: 783, 784, 783, 784, 784, 785, ... 
## Addtional sampling using down-sampling prior to pre-processing
## 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa    
##   2     0.6674173  0.2977857
##   5     0.6734856  0.3091695
##   8     0.6690066  0.3013622
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 5.

final_u <- data.frame(Real = test$remocao,
                      predict(model_und, newdata = test, type = "prob"))

final_u$predict <- ifelse(final_u$nao > 0.5, "nao", "sim")

final_u$predict <- as.factor(final_u$predict)
test$remocao    <- as.factor(test$remocao)

# levels(final_u$predict)
# levels(test$remocao)

cm_u <- confusionMatrix(final_u$predict, test$remocao)

# print(paste("Matriz de confusão do resampling uder-sampling"))
# cm_u

No chunck abaixo está o modelo com sobreamostragem (oversampling).

# Modelo com resampling oversampling

{
  set.seed(1)
  
  index <- createDataPartition(saidax$remocao, p = 0.75, list = FALSE) 
  train <- saidax[index, ]
  test  <- saidax[-index, ]

  ctrl <- trainControl(method  = "repeatedcv", 
                       number  = 10, 
                       repeats = 10,
                       verboseIter = FALSE,
                       sampling   = "up",
                       classProbs = T, 
                       savePredictions = T)

  model_ov <- caret::train(remocao ~ .,
                           data   = train,
                           method = "rf",
                           preProcess = c("scale", "center"),
                           trControl  = ctrl,
                           ntree  = 1000)
  }

model_ov

## Random Forest 
## 
## 871 samples
##   8 predictor
##   2 classes: 'nao', 'sim' 
## 
## Pre-processing: scaled (8), centered (8) 
## Resampling: Cross-Validated (10 fold, repeated 10 times) 
## Summary of sample sizes: 783, 784, 783, 784, 784, 785, ... 
## Addtional sampling using up-sampling prior to pre-processing
## 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa    
##   2     0.7397268  0.3056056
##   5     0.7325748  0.3111261
##   8     0.7291629  0.3066799
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 2.

final_ov <- data.frame(Real = test$remocao,
                       predict(model_ov, newdata = test, type = "prob"))

final_ov$predict <- ifelse(final_ov$nao > 0.5, "nao", "sim")

final_ov$predict <- as.factor(final_ov$predict)
test$remocao     <- as.factor(test$remocao)

# levels(final_ov$predict)
# levels(test$remocao)

cm_ov <- confusionMatrix(final_ov$predict, test$remocao)

# print(paste("Matriz de confusão do resampling oversampling"))
# cm_ov

No chunck abaixo está o modelo com SMOTE (Synthetic Minority Over-sampling Technique).

# Modelo com resampling SMOTE

{
  set.seed(1)
  
  index <- createDataPartition(saidax$remocao, p = 0.75, list = FALSE) 
  train <- saidax[index, ]
  test  <- saidax[-index, ]

  ctrl <- trainControl(method  = "repeatedcv", 
                       number  = 10, 
                       repeats = 10,
                       verboseIter = FALSE,
                       sampling   = "smote",
                       classProbs = T, 
                       savePredictions = T)

  model_s <- caret::train(remocao ~ .,
                          data   = train,
                          method = "rf",
                          preProcess = c("scale", "center"),
                          trControl  = ctrl,
                          ntree  = 1000)
  }

model_s

## Random Forest 
## 
## 871 samples
##   8 predictor
##   2 classes: 'nao', 'sim' 
## 
## Pre-processing: scaled (8), centered (8) 
## Resampling: Cross-Validated (10 fold, repeated 10 times) 
## Summary of sample sizes: 783, 784, 783, 784, 784, 785, ... 
## Addtional sampling using SMOTE prior to pre-processing
## 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa    
##   2     0.7006602  0.3075519
##   5     0.7018463  0.3078595
##   8     0.6952991  0.2986976
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 5.

final_s <- data.frame(Real = test$remocao,
                      predict(model_s, newdata = test, type = "prob"))

final_s$predict <- ifelse(final_s$nao > 0.5, "nao", "sim")

final_s$predict <- as.factor(final_s$predict)
test$remocao    <- as.factor(test$remocao)

# levels(final_s$predict)
# levels(test$remocao)

cm_s <- confusionMatrix(final_s$predict, test$remocao)

# print(paste("Matriz de confusão do resampling SMOTE"))
# cm_s

models <- list(original = model_orig,
               under = model_und,
               over  = model_ov,
               smote = model_s)

comparacao <- data.frame(Modelo = names(models),
                         Sensitividade  = rep(NA, length(models)),
                         Especificidade = rep(NA, length(models)),
                         Precisao = rep(NA, length(models)),
                         Recall   = rep(NA, length(models)),
                         F1       = rep(NA, length(models)))

comparacao <- comparacao %>% 
  dplyr::mutate(Sensitividade  = c(cm_orig$byClass[1],
                                   cm_u$byClass[1],
                                   cm_ov$byClass[1],
                                   cm_s$byClass[1]),
                Especificidade = c(cm_orig$byClass[2],                               
                                   cm_u$byClass[2],
                                   cm_ov$byClass[2],
                                   cm_s$byClass[2]),
                Precisao       = c(cm_orig$byClass[5],
                                   cm_u$byClass[5],
                                   cm_ov$byClass[5],
                                   cm_s$byClass[5]),
                Recall         = c(cm_orig$byClass[6],
                                   cm_u$byClass[6],
                                   cm_ov$byClass[6],
                                   cm_s$byClass[6]),
                F1             = c(cm_orig$byClass[7],
                                   cm_u$byClass[7],
                                   cm_ov$byClass[7],
                                   cm_s$byClass[7]),
                Acuracia       = c(max(model_orig$results[2]),
                                   max(model_und$results[2]),
                                   max(model_ov$results[2]),
                                   max(model_s$results[2]))) %>% 
  reshape2::melt(id.vars = c("Modelo"))

# levels(comparacao$variable)
levels(comparacao$variable) <- c("Sensitividade", "Especificidade", "Precisão",
                                 "Recall", "F1", "Acurácia")

shiny::div(plotly::ggplotly(ggplot2::ggplot(comparacao,
                                            (aes(x = variable,
                                                 y = value, 
                                                 color = Modelo))) +
                              geom_jitter(width = 0.2, alpha = 0.5, size = 3) +
                              scale_y_continuous(limits = c(0, 1)) +
                              labs(title = "Modelo de classificação para DBO de saída",
                                   x = "Índice",
                                   y = "Valor",
                                   color = "Legenda") +
                              theme_bw()), align = "center")

Como esperado, pode-se observar que os modelos propostos sofrem influência das técnicas de reamostragem. E o melhor modelo para classificação de DBO de saída é com resampling SMOTE (Synthetic Minority Over sampling Technique), comparado com as outras situações, apresenta bom desempenho e equilíbrio nos índices aqui analisados.
Os outros métodos de resampling, no geral, também apresentam resultados satisfatórios, o modelo com dados originais apresenta um desequilíbrio, tendo valores de sensitividade , recall e F1 muito baixo. A escolha do modelo smote, mesmo apresentando menor acurácia, é indicado pois em problemas com classes desproporcionais, ela causa uma falsa impressão de bom desempenho.
Este algoritmo escolhido aplica a abordagem KNN onde ele seleciona K vizinhos mais próximos, une-os e cria as amostras sintéticas no espaço.

Discussão

A estratégia de predição aqui realizada teve como objetivo principal predizer/estimar valores do parâmetro de qualidade da água mais importantes para avaliação de condicionantes ambientais. A DBO de entrada e de saída do sistemas de lagoas aeradas demoram até cinco dias para serem analisadas em laboratório para estimativas com baixa incerteza, entretanto o tempo de detenção hidráulica do sistema em questão é muito menor que o tempo necessário para a análise, ou seja, considerando que o efluente que entrou no sistema no dia 1 só será possível analisar se o sistema de lagoas aeradas cumpriu com seu objetivo (redução de DBO em 60%) no dia 7. Pois apenas no dia 5 será possível obter resultados da DBO5 do efluente de entrada e no dia 7 (2 dias de tempo de detenção hidráulica + 5 dias em análise) a DBO5 do efluente de saída da lagoa aerada. Logo existe um delay/janela de 6 dias de vulnerabilidade do sistema, onde a empresa pode estar lançando efluente de forma inadequada, fora do padrão exigido.
O modelo de predição de DBO5 de entrada e DBO5 de saída, suplementado pelo preditor de DQO de saída (bom parâmetro de correlação com DBO) possibilitam que a empresa planeje e execute possíveis ajustes no sistema de aeração da lagoa ou implemente o uso de algum agente químico para oxidação. Pois os sistemas de aeração têm, basicamente, a função de incrementar oxigênio à água em tratamento e essa aeração é responsável pelo fornecimento do oxigênio necessário ao desenvolvimento das reações para o tratamento de efluentes.
O conhecimento do regime de operação fabril, da vazão inflow, bem como da qualidade da água do efluente que chega nas lagoas é fundamental para a avaliação de possíveis mudanças abruptas no sistema, da eficiência do processo e do controle do oxigênio dissolvido, se a aeração for passível de controle e incremento.
Quanto aos agrupamentos, para os parâmetros de efluente de entrada, os mesmo demonstram potencial para melhorar os modelos preditores, entretanto a ausência de variáveis como Ocorrência e Volume diário de chuva, e Temperatura do efluente limitam o potencial de criação de clusters coesos numa abordagem não supervisionada (unsupervised learning). O modelo de saída de classificação em “sim” ou “não” quanto a adequabilidade de lançamento de efluente no corpo receptor já é uma abordagem supervisionada (supervised learning).
Em relação a possíveis trabalhos futuros, fica a sugestão de análise de capacidade produtiva da unidade fabril (produção de polpa de celulosa e produção de papel) tendo em vista o enquadramento do efluente para lançamento no corpo hídrico e modelo de controlo preditivo (MPC) baseado no modelo de classficicação aqui proposto.

---

Discentes:
Brenner Silva;
Herica Oliveira;
Lucas Mascarenhas.

Docentes:
Cristiano Fontes;
Karla Esquerre.