Introdução

Este trabalho se dedica a analisar os dados de um sistema de lagoas aeradas utilizadas para tratar o resíduo da produção de papel e celulose de uma planta industrial. Diversos parâmetros de qualidade são medidos na entrada e saída do sistema, entretanto iremos remover as variáveis de entrada que possuem menos de 300 observações, uma vez que nosso intuito não é preencher dados faltantes. Assim, das 1430 observações originais permaneceram 971 completas.

##        Datein    FR(m3/dia)   BODin (ppm)   CODin (ppm)     SSin(ppm) 
##             0             0            89            89           862 
##          pHin    NAmin(ppm)    NNin (ppm)     Pin (ppm)    Colin(ppm) 
##            53           770          1151          1170            51 
##       Tin(°C)     Condin L1        RF(mm) Pulp(ton/dia) Pap (ton/dia) 
##           467            56           255           103            92

As variáveis de entrada restantes, bem como a variável de saída de interesse (\(\small BOD_{OUT}\)), são reunidas num único banco de dados. As variáveis a serem estudadas são:

\(\small FR_{IN}\): Vazão de efluente em \(\small m^{3}/dia\)
\(\small BOD_{IN}\), \(\small BOD_{OUT}\): A demanda bioquímica de oxigênio é definida como a quantidade de oxigênio consumida pela respiração de microorganismos presentes em uma amostra usualmente mantida a 20ºC por cinco dias. Comumente medida em mg de \(\small 0_{2}\) por litro. É um estimador da quantidade de matéria orgânica biodegradável presente na amostra
\(\small COD_{IN}\): Corresponde à quantidade de oxigênio necessário para oxidar a matéria orgânica numa amostra em meio ácido e aquecimento por meio de um agente oxidante forte. A maior parte das espécies orgânicas é oxidada, biodegradáveis ou não, todavia nem todas espécies são consumidas. A razão \(\small COD/BOD\) é uma medida de quão biodegradável é o resíduo líquido. Quão menor, mais biodegradável, sendo uma razão entre 1,5 e 2 um valor satisfatório
\(\small pH_{IN}\)
\(\small Col_{IN}\): Cor da amostra. Geralmente comparada relativamente a uma solução de cloroplatinado de potássio e cloreto de cobalto em água deionizada
\(\small RF_{IN}\): Precipitação em mm de coluna de água
\(\small Pulp\): Produção de celulose em ton/dia
\(\small pH_{IN}\): Produção de papel em ton/dia

Cientes dos objetos de nosso estudo, prosseguimos para a análise exploratória das variáveis.

Análise exploratória

Nesta etapa, para cada variável é apresentado um histograma, um boxplot, uma tabela com estatísticas amostrais relevantes, bem como o resultado do teste de normalidade de Shapiro-Wilk (cuja hipótese nula é que a distribuição amostral é normal).

Vazão de entrada

Estatísticas relevantes
Min.	1st Qu.	Median	Mean	3rd Qu.	Max.
5913	60495.5	64438	64906.64	71647	97850

Teste de normalidade de Shapiro-Wilk
W	p.valor	Resultado
0.9010198	0	Hipótese nula rejeitada

A priori nota-se um rígido regime da vazão de entrada no sistema, grande parte dos dados (70%) encontra-se entre 60.000 e 80.000 \(\small m^{3}/dia\). Todavia a série temporal parece sugerir que em verdade houve uma mudança de regime entre janeiro e julho de 1998 (destaque em azul).

Tal suspeita é ampliada ao separar o conjunto de dados entre antes de janeiro de 98 e após julho de 98. Os histogramas abaixo exibem certa diferença entre os regimes. Tal diferença é confirmada através de um teste de hipóteses (cujo p-valor é aproximadamente zero).

## 
##  Welch Two Sample t-test
## 
## data:  a$`FR(m3/dia)` and b$`FR(m3/dia)`
## t = 26.712, df = 294.77, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  15134.07 17541.47
## sample estimates:
## mean of x mean of y 
##  76995.69  60657.92

Se levanta então a questão: tal redução significativa na vazão de entrada no sistema de tratamento da lagoa aerada é decorrente de redução na produção?

Produção de celulose e de papel

Estatísticas relevantes
Min.	1st Qu.	Median	Mean	3rd Qu.	Max.
-48	865.2765	922.332	888.4425	967.191	1112.094

Teste de normalidade de Shapiro-Wilk
W	p.valor	Resultado
0.6772969	0	Hipótese nula rejeitada

Nota-se um regime um pouco mais flexível comparado ao regime de vazão de entrada de efluente, certa parte dos dados (60%) encontra-se entre 800 e 960 ton/dia. Já a série temporal não indica qualquer mudança de regime.

Para confirmar que não há relação entre a mudança de regime de vazão de entrada de efluente e produção de celulose fazemos os histogramas dos períodos antes e depois da alteração, bem como efetuamos o teste de hipóteses. O teste confirma que não há diferença significativa, pois o p-valor é superior a 0,05.

## 
##  Welch Two Sample t-test
## 
## data:  a$`Pulp(ton/dia)` and b$`Pulp(ton/dia)`
## t = -0.91449, df = 252.14, p-value = 0.3613
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -37.78368  13.82114
## sample estimates:
## mean of x mean of y 
##  880.9722  892.9535

Prosseguimos, então, para o outro suspeito de provocar a mudança observada: a massa de produção de papel.

Estatísticas relevantes
Min.	1st Qu.	Median	Mean	3rd Qu.	Max.
382.4	1008.35	1060.4	1045.932	1107.65	1304.8

Teste de normalidade de Shapiro-Wilk
W	p.valor	Resultado
0.8963636	0	Hipótese nula rejeitada

Não existem grandes novidades nesta variável. A série temporal não indica qualquer mudança de regime. A quantidade de outliers é menor do que a encontrada para a massa de celulose produzida.

Para confirmar que não há relação entre a mudança de regime de vazão de entrada de efluente e produção de papel fazemos os histogramas dos períodos antes e depois da alteração, bem como efetuamos o teste de hipóteses. O teste aponta diferença significativa, pois o p-valor é superior a 0,05. Entretanto em termos práticos a diferença parece pequena se comparada a magnitude da produção (limite superior da diferença é 44,5 frente a uma média de 1060,4 ton/dia).

## 
##  Welch Two Sample t-test
## 
## data:  a$`Pap (ton/dia)` and b$`Pap (ton/dia)`
## t = -3.6488, df = 273.46, p-value = 0.0003155
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -44.59679 -13.33832
## sample estimates:
## mean of x mean of y 
##  1022.366  1051.333

DBO de entrada e de saída, DQO de entrada

Estatísticas relevantes
Min.	1st Qu.	Median	Mean	3rd Qu.	Max.
41	223	247	247.3164	272	428

Teste de normalidade de Shapiro-Wilk
W	p.valor	Resultado
0.9793191	0	Hipótese nula rejeitada

Estatísticas relevantes
Min.	1st Qu.	Median	Mean	3rd Qu.	Max.
16	72	86	86.39765	100	187

Teste de normalidade de Shapiro-Wilk
W	p.valor	Resultado
0.9877907	3e-07	Hipótese nula rejeitada

Estatísticas relevantes
Min.	1st Qu.	Median	Mean	3rd Qu.	Max.
136	528	585	585.2472	639	924

Teste de normalidade de Shapiro-Wilk
W	p.valor	Resultado
0.9810531	0	Hipótese nula rejeitada

É interessante notar que para as três variáveis a média e mediana não diferem muito (evidência de distribuição normal).

Os perfis temporais das variáveis apresentam semelhanças visuais, mas a correlação linear precisa de confirmação.

Correlação entre variáveis com lag
lag	BODin_x_BODout	BODin_x_CODin	BODout_x_CODin
0	0.4863904	0.6051251	0.4134975
1	0.4819040	0.4219760	0.3763221
2	0.4589478	0.3388424	0.3124284
3	0.3346562	0.2910248	0.2381219
4	0.3232135	0.2722827	0.2066645
5	0.2746676	0.2259565	0.1903808
6	0.2568760	0.2245563	0.1625922
7	0.2693774	0.2306913	0.1375835
8	0.2772512	0.2305730	0.1304655
9	0.2571072	0.2190677	0.1304306
10	0.2702836	0.1965503	0.1283663

A tabela indica que a máxima correlação entre as variáveis ocorre na ausência de lag, isto é, os valores pareados estão vinculados à mesma data. Todavia as correlações não são tão intensas. Outro objeto de interesse na análise inicial são as razões \(\small COD_{IN}/BOD_{IN}\) (uma medida da não-biodegradabilidade do efluente que entre no sistema) e \(\small (BOD_{IN}-BOD_{OUT})/BOD_{IN}\) (uma medida da eficiência do sistema de remoção de matéria orgânica).

Ambas as razões parecem ter um movimento sinuoso a partir de 1999, sendo que a não-biodegradabilidade parece aumentar aos poucos se afastando da zona para que um efluente líquido seja considerado biodegradável (destaque em vermelho). Além disso, mais uma vez, a mudança de vazão de efluente não parece se relacionar com variável alguma.

pH e coloração

Estatísticas relevantes
Min.	1st Qu.	Median	Mean	3rd Qu.	Max.
4.22	6.85	7.12	7.45402	7.46	12.53

Teste de normalidade de Shapiro-Wilk
W	p.valor	Resultado
0.7480285	0	Hipótese nula rejeitada

Estatísticas relevantes
Min.	1st Qu.	Median	Mean	3rd Qu.	Max.
41	395	480	462.4552	533	911

Teste de normalidade de Shapiro-Wilk
W	p.valor	Resultado
0.9765562	0	Hipótese nula rejeitada

o pH tem uma forte constância na região neutra, mas rotineiramente acontencem picos de alcalinidade (observados como outliers no boxplot). Já a coloração apresenta uma mudança de regime a partir de fevereiro/2000.

Precipitação

Estatísticas relevantes
Min.	1st Qu.	Median	Mean	3rd Qu.	Max.
0	0	0.1	4.398146	2.8	175.4

Teste de normalidade de Shapiro-Wilk
W	p.valor	Resultado
0.4331416	0	Hipótese nula rejeitada

Como esperado, a precipitação é sazonal, tendo picos entre janeiro e março. Entretanto, as chuvas não são muito volumosas.

Random Forest

Serão testados modelos de árvores de decisões, mais especificamente, random forest. O trabalho se dá em três etapas principais:

1 - Agrupamento dos dados utilizando da random forest não supervisionada (variando o número de agrupamentos);

2 - Validação e aprendizado do agrupamento acima, utilizando como label o agrupamento sugerido anteriormente (variando o mtry)

3 - Cruzamento com a variável Dbo de saída, considerando como tempo de detenção 2 dias. Desejamos aqui que o número de grupos e os parâmetros da random forest sejam escolhidos de forma a criar grupos, com maior nível de pureza possível, e que distinguam de forma satisfatória os níveis de DBO de saida.

Desta forma, busca-se compreender, também, como funciona o reforçamento do aprendizado de máquinas. Uma vez que a variável que define a adequalidade do modelo não está presente na etapa de treino e teste, sendo, então, adicionada após o aprendizado como se fosse uma condição de contorno externa ao modelo.

comp <- matrix(ncol = 4, nrow = 4)

for(n in 2:5){ 
  
## Random Forest não supervisionada 

  
clust <- rf.unsupervised(x[,-c(9:10)], n=n, proximity = TRUE)

for(m in 2:5){ 
  
new_x <- cbind(x, clust$k)  

## Random Forest supervisionada 

{set.seed(3)

mtry = m
  
control <- trainControl(method="repeatedcv", number=10, repeats=3)

tunegrid <- expand.grid(.mtry=mtry)

g_new <- train(x=new_x[,-c(9:11)], 
                      y=as.factor(new_x[,c(11)]), 
               method="rf", 
               trControl=control,
               tuneGrid=tunegrid) 

g_new


prob <- predict(g_new, new_x[,-c(9:11)], type="raw")

aux <- table(prob)

}
}
}

Cruzando com DBO out

ind <- new_x$`BODout(ppm)`

dbo_out <- cbind(prob,ind) %>% 
  as.tibble() %>% 
  dplyr::arrange() %>% 
  mutate(dbo_class = 0 )

## Warning: `as.tibble()` is deprecated, use `as_tibble()` (but mind the new semantics).
## This warning is displayed once per session.

# Criando n grupos sequenciados do menor para o maior e com números de elementos igual a proporção de classificação

if(n == 2){
  
  dbo_out <- dbo_out %>% 
    mutate(dbo_class = c(rep(1,aux[1]), rep(2, aux[2])))
 
  }else if (n == 3) {
  
  dbo_out <- dbo_out %>% 
    mutate(dbo_class = c(rep(1,aux[1]), rep(2, aux[2]),rep(3, aux[3])))

  }else if (n == 4) {
    
    
  dbo_out <- dbo_out %>% 
      mutate(dbo_class = c(rep(1,aux[1]), rep(2, aux[2]),rep(3, aux[3]), rep(4, aux[4]))) 
  
  }else{ 
    
    dbo_out <- dbo_out %>% 
      mutate(dbo_class = c(rep(1,aux[1]), rep(2, aux[2]),rep(3, aux[3]), rep(4, aux[4]),rep(5, aux[5])))
    
    }

colnames(dbo_out) <- c("cluster", 'value', "class")

comparation <- funtimes::purity(as.factor(dbo_out$cluster), as.factor(dbo_out$class))

comp[n-1,m-1] <- comparation$pur

Escolha

comp

##      [,1] [,2] [,3]     [,4]
## [1,]   NA   NA   NA       NA
## [2,]   NA   NA   NA       NA
## [3,]   NA   NA   NA       NA
## [4,]   NA   NA   NA 0.569516

Com 2 grupos e mtry = 2, obtemos o melhor agrupamento, devido a relação de grau de pureza, sendo, então, a combinação mais simples que produz uma melhor separação da DBO de saída. Temos, assim, dois grupos, que apresentam uma boa separação em dois estados de DBO de saída. Podendo, prever, então, dois estados de operação com riscos associados (sugestão de trabalhos futuros).

Refazendo o modelo

clust <- rf.unsupervised(x[,-c(9:10)], n=2, proximity = TRUE, mtry = 2)

x <- cbind(x, clust$k)

mtry = 2
  
control <- trainControl(method="repeatedcv", number=10, repeats=3)

tunegrid <- expand.grid(.mtry=mtry)

g_new <- train(x=x[,-c(9:11)], 
                      y=as.factor(x[,c(11)]), 
               method="rf", 
               trControl=control,
               tuneGrid=tunegrid) 


g_new$finalModel$importance

##               MeanDecreaseGini
## FR(m3/dia)           35.673661
## BODin (ppm)          25.327493
## CODin (ppm)          39.965427
## pHin                  5.648882
## Colin(ppm)           49.095141
## RF(mm)                4.549629
## Pulp(ton/dia)         7.308811
## Pap (ton/dia)         4.914173

A árvore de decisão média, composta pelas árvores de decisão distintas, apresentaram como variáveis mais pertinentes: Colin(ppm); FR(m³/dia); CODin(ppm) e BODin(ppm). É interessante, também, realizar a ligação destas variáveis de importância com a classificação final como é proposto no artigo: “Why Should I Trust You?” dos atores: Marco Túlio Ribeiro, Sammer Singh e Carlos Guestrin.

Histograma da DBO de saida

DT::datatable(x, options = list(
  initComplete = JS(
    "function(settings, json) {",
    "$(this.api().table().header()).css({'background-color': '#000', 'color': '#fff'});",
    "}")
))

x %>% 
  ggplot(aes(x =`BODout(ppm)`)) +
  geom_density(aes(fill = as.factor(`clust$k`)), alpha = 0.4) +
  theme_bw()

Relação entre Saida e Entrada

x %>% 
  ggplot(aes(x =`BODout(ppm)`, y = (x$`BODout(ppm)`/x$`BODin (ppm)`))) +
  geom_point(aes(col = as.factor(`clust$k`)), alpha = 0.4) +
  coord_flip()+
  theme_bw()

# Observamos que no grupo 2 estão as menores eficiencias para os mesmos valores de DBO

Estas observações, sugerem uma boa separação entre os grupos, sendo a príncipio um bom ponto de partida para a divisão da operação em dois estados e escolha daquele que se opera com maior confiabilidade. É interessante, também, observar a separação após métodos de redução de dimensionalidade, reduzindo a dimensão das variáveis de entrada. Para isto, têm-se classicamente a Análise de Componentes Principais (PCA), mas é interessante que se invista, também, em algoritmos de redução de dimensionalidade mais novos como o t- SNE (http://www.jmlr.org/papers/v9/vandermaaten08a.html).

Algumas outras observações

x %>% 
  ggplot(aes(x =x$`BODout(ppm)`/x$`FR(m3/dia)`)) +
  geom_density(aes(fill = as.factor(`clust$k`)), alpha = 0.4) +
  theme_bw()

Temos uma separação interessante da relação entre a DBO de saída e a vazão de entrada, gerando dois agrupamentos MUITO CLARAMENTE!

Sugestões para estudos futuros

x %>% 
  ggplot(aes(x = x$`Pap (ton/dia)`/x$`FR(m3/dia)`)) +
  geom_density(aes(fill = as.factor(`clust$k`)), alpha = 0.4) +
  theme_bw()

Precisamos verificar melhor, mas os dados tem sugerido que temos um grupo com menos elementos que possuem as maiores produções conseguemente maior DBO, porém baixas eficiencias. Pode-se pensar em estudar este módulo de operação com maior risco.

Referências

WOODARD, Frank. Industrial waste treatment handbook. Elsevier, 2001.
VALENTE, José Pedro Serra; PADILHA, Pedro Magalhães; SILVA, Assunta Maria Marques. Oxigênio dissolvido (OD), demanda bioquímica de oxigênio (DBO) e demanda química de oxigênio (DQO) como parâmetros de poluição no ribeirão Lavapés/Botucatu-SP. Eclética Química, p. 49-66, 1997.
JOUANNEAU, S. et al. Methods for assessing biochemical oxygen demand (BOD): A review. Water research, v. 49, p. 62-82, 2014.

Trabalho DBO - DQO

Icaro Bernardes, Isabela Almeida, Tarssio Barreto & Vitor Quintella