Este trabalho apresenta um modelo de otimização de manutenção para uma frota de transformadores de alta potência. Utilizando dados históricos de falhas, modelamos o processo de degradação através de um Processo de Poisson Não Homogêneo (NHPP) na forma de Lei de Potência (PLP). Os parâmetros do modelo, estimados por Máxima Verossimilhança, indicam um claro processo de desgaste (\(\beta > 1\)). Com base neste modelo e em uma análise de custo, determinamos uma política de manutenção preventiva ótima que minimiza os custos operacionais a longo prazo, estabelecendo um intervalo de substituição de aproximadamente 6286 horas.
2 Introdução
O objetivo deste trabalho é determinar uma política ótima de manutenção para um conjunto de transformadores de alta potência de uma empresa de energia elétrica. Após uma falha, os custos de reparo são elevados e a companhia pode ser multada por interrupções não programadas no fornecimento de energia [@fonteEmpresaAno]. Uma estratégia de manutenção preventiva bem definida é crucial para otimizar recursos e aumentar a confiabilidade do sistema elétrico [@normaSetorAno].
Nesta análise, consideramos dois tipos de intervenção:
Reparo Mínimo (MR - Minimal Repair): Realizado após uma falha. Restaura a função do equipamento, mas o retorna à condição de confiabilidade exata de antes da falha, ou seja, “tão ruim quanto antes” (bad as old) [@barlow1975; @brown1981].
Manutenção Preventiva (PM - Preventive Maintenance): Ação planejada e completa que retorna o sistema a um estado de “como novo” (good as new) [@nakagawa2007]. Para encontrar o equilíbrio ideal entre essas intervenções, é fundamental modelar matematicamente como as falhas ocorrem ao longo do tempo.
3 Modelagem do Processo de Falhas
Para descrever as falhas que ocorrem entre as manutenções preventivas, utilizamos um modelo estatístico para sistemas reparáveis.
3.1 Modelo Adotado: Processo de Poisson Não Homogêneo (NHPP)
O NHPP é ideal para modelar eventos (falhas) que ocorrem ao longo do tempo com uma taxa que não é constante [@rigdon2000]. A taxa de falhas, ou função de intensidade, é denotada por \(\rho(t)\).
A forma específica utilizada é o Power Law Process (PLP), que é a implementação mais popular do NHPP para análise de confiabilidade de sistemas reparáveis [@crow1991; @asher1984]. A função de intensidade do PLP é dada por:
\[
\rho(t) = \frac{\beta}{\theta}\left(\frac{t}{\theta}\right)^{\beta-1}
\] - Os parâmetros \(\beta\) (forma) e \(\theta\) (escala) são análogos aos da distribuição de Weibull. - Significado de \(\beta\): Se \(\beta > 1\), a taxa de falhas aumenta com o tempo, indicando um processo de desgaste e justificando a necessidade de manutenção preventiva. Se \(\beta < 1\), a confiabilidade melhora com o tempo (mortalidade infantil). Se \(\beta = 1\), as falhas são constantes (aleatórias), caracterizando um Processo de Poisson Homogêneo.
4 Análise dos Dados
4.1 Estimação dos Parâmetros do Modelo
O primeiro passo é estimar os parâmetros \(\beta\) e \(\theta\) a partir dos dados históricos de falha dos transformadores. Utilizamos o Método da Máxima Verossimilhança (MLE), que encontra os valores dos parâmetros que maximizam a probabilidade de observar os dados coletados.
Parâmetros Estimados para o Modelo PLP via Máxima Verossimilhança.
Parâmetro
Valor.Estimado
Beta (Forma)
1.9955
Theta (Escala)
24363.2130
O valor de \(\beta\) estimado em aproximadamente 1.996 é crucial para esta análise. Por ser significativamente maior que 1, ele confirma estatisticamente que os transformadores estão em um processo de envelhecimento e desgaste, onde a probabilidade de falha aumenta à medida que o equipamento opera por mais tempo. É este comportamento que justifica economicamente a implementação de uma política de manutenção preventiva, em vez de simplesmente reparar os equipamentos conforme eles falham.
4.2 Avaliação da Adequação do Modelo
Para validar nosso modelo, comparamos a Mean Cumulative Function (MCF), que representa a média de falhas acumuladas ao longo do tempo. Comparamos a curva não-paramétrica (calculada diretamente dos dados via estimador de Nelson-Aalen) com a curva paramétrica (prevista pelo nosso modelo PLP com os parâmetros estimados). Se as curvas forem próximas, o modelo é considerado adequado.
A proximidade entre a curva paramétrica (vermelha), prevista por nosso modelo, e a curva não-paramétrica (azul), calculada diretamente dos dados, é notável. O modelo PLP consegue capturar com precisão a dinâmica das falhas: tanto o período inicial de baixa ocorrência quanto a aceleração do processo de desgaste após as 10.000 horas de operação. Essa validação gráfica nos dá um alto grau de confiança para usar este modelo na determinação da política ótima de manutenção.
5 Determinação da Política Ótima de Manutenção
O objetivo final é encontrar o intervalo de manutenção preventiva, \(\tau\), que minimiza o custo esperado por unidade de tempo, \(H(\tau)\). A função de custo é definida como [@gilardoni2007]:
\[H(\tau) = \frac{C_{PM} + C_{MR} \cdot M(\tau)}{\tau}
\]Onde \(M(\tau) = (\tau / \theta)^\beta\) é a MCF do nosso modelo e a razão de custos entre reparo mínimo e manutenção preventiva é definida como \(C_{MR}/C_{PM} = 15\). Para o modelo PLP, a solução ótima \(\tau^*\) que minimiza \(H(\tau)\) pode ser encontrada analiticamente pela fórmula:
Aplicando a fórmula com os parâmetros estimados, encontramos o intervalo ótimo e visualizamos a função de custo para confirmar o ponto de mínimo.
Código
C_PM <-1; C_MR <-15tau_otimo <- theta_estimado_correto * (C_PM / ((beta_estimado_correto -1) * C_MR))^(1/ beta_estimado_correto)funcao_custo_H <-function(tau, beta, theta, cpm, cmr) { (cpm + cmr * (tau / theta)^beta) / tau }custo_otimo <-funcao_custo_H(tau_otimo, beta_estimado_correto, theta_estimado_correto, C_PM, C_MR)intervalo_tau <-seq(tau_otimo -800, tau_otimo +800, length.out =500)dados_custo <-data.frame(tau = intervalo_tau, custo =funcao_custo_H(intervalo_tau, beta_estimado_correto, theta_estimado_correto, C_PM, C_MR))ggplot(dados_custo, aes(x = tau, y = custo)) +geom_line(color ="blue", size =1) +geom_vline(xintercept = tau_otimo, color ="red", linetype ="dashed", size =1) +annotate("point", x = tau_otimo, y = custo_otimo, color ="red", size =4) +annotate("text", x = tau_otimo, y = custo_otimo + (max(dados_custo$custo, na.rm =TRUE) -min(dados_custo$custo, na.rm =TRUE))*0.05, label =paste0("Ponto Ótimo\nτ* ≈ ", round(tau_otimo), " horas\nH(τ*) ≈ ", format(custo_otimo, scientific =FALSE, digits =5)),hjust =0.5, vjust =-0.5, size=3.5 ) +labs(title ="Função de Custo da Manutenção Preventiva (Focado no Ótimo)",x ="Intervalo de Manutenção Preventiva τ (horas)",y ="Custo por Unidade de Tempo H(τ)") +theme_bw() +theme(panel.grid.major =element_line(linetype ='dotted'),panel.grid.minor =element_line(linetype ='dotted'))
Política Ótima de Manutenção:
Política de Manutenção Preventiva Ótima.
Métrica
Valor Recomendado
Intervalo Ótimo (τ*) em Horas
6286
Intervalo Ótimo (τ*) em Dias
262
O intervalo ótimo de 6286 horas representa o ponto de equilíbrio econômico ideal. Realizar a manutenção antes deste marco seria prematuro e implicaria em custos desnecessários com a substituição de equipamentos ainda confiáveis. Por outro lado, adiar a manutenção para além deste ponto aumenta exponencialmente o risco de falhas em serviço, cujos custos de reparo e multas superam a economia obtida ao estender a vida útil do equipamento. Portanto, este valor não é apenas um resultado estatístico, mas uma recomendação gerencial direta.
6 Conclusão
O processo de falha dos transformadores é bem descrito por um modelo de desgaste (PLP com \(\beta > 1\)), validando a necessidade de manutenção preventiva.
A análise fornece uma política de manutenção orientada por dados: realizar uma manutenção preventiva completa a cada 6286 horas (ou aproximadamente 262 dias).
Adotar este intervalo ótimo permite à empresa minimizar os custos operacionais totais, equilibrando os gastos com manutenções planejadas e os custos mais elevados de reparos corretivos não planejados.
Código fonte
---title: "Análise de Sobrevivência e Otimização de Manutenção para Transformadores"format: html: toc: true toc-title: "Conteúdo" number-sections: true self-contained: true code-fold: true code-tools: true revealjs: slide-level: 2 theme: sky self-contained: true pdf: toc: true number-sections: truelang: pt---## ResumoEste trabalho apresenta um modelo de otimização de manutenção para uma frota de transformadores de alta potência. Utilizando dados históricos de falhas, modelamos o processo de degradação através de um Processo de Poisson Não Homogêneo (NHPP) na forma de Lei de Potência (PLP). Os parâmetros do modelo, estimados por Máxima Verossimilhança, indicam um claro processo de desgaste ($\beta > 1$). Com base neste modelo e em uma análise de custo, determinamos uma política de manutenção preventiva ótima que minimiza os custos operacionais a longo prazo, estabelecendo um intervalo de substituição de aproximadamente 6286 horas.```{r setup, include=FALSE}library(dplyr)library(ggplot2)library(survival)knitr::opts_chunk$set(echo =TRUE, message =FALSE, warning =FALSE, fig.align ='center')```## IntroduçãoO objetivo deste trabalho é determinar uma política ótima de manutenção para um conjunto de transformadores de alta potência de uma empresa de energia elétrica. Após uma falha, os custos de reparo são elevados e a companhia pode ser multada por interrupções não programadas no fornecimento de energia [@fonteEmpresaAno]. Uma estratégia de manutenção preventiva bem definida é crucial para otimizar recursos e aumentar a confiabilidade do sistema elétrico [@normaSetorAno].Nesta análise, consideramos dois tipos de intervenção: - **Reparo Mínimo (MR - *Minimal Repair*):** Realizado após uma falha. Restaura a função do equipamento, mas o retorna à condição de confiabilidade exata de antes da falha, ou seja, "tão ruim quanto antes" (*bad as old*) [@barlow1975; @brown1981]. - **Manutenção Preventiva (PM - *Preventive Maintenance*):** Ação planejada e completa que retorna o sistema a um estado de "como novo" (*good as new*) [@nakagawa2007]. Para encontrar o equilíbrio ideal entre essas intervenções, é fundamental modelar matematicamente como as falhas ocorrem ao longo do tempo.## Modelagem do Processo de FalhasPara descrever as falhas que ocorrem entre as manutenções preventivas, utilizamos um modelo estatístico para sistemas reparáveis.### Modelo Adotado: Processo de Poisson Não Homogêneo (NHPP)O NHPP é ideal para modelar eventos (falhas) que ocorrem ao longo do tempo com uma taxa que não é constante [@rigdon2000]. A taxa de falhas, ou **função de intensidade**, é denotada por $\rho(t)$.A forma específica utilizada é o **Power Law Process (PLP)**, que é a implementação mais popular do NHPP para análise de confiabilidade de sistemas reparáveis [@crow1991; @asher1984]. A função de intensidade do PLP é dada por:$$\rho(t) = \frac{\beta}{\theta}\left(\frac{t}{\theta}\right)^{\beta-1}$$ - Os parâmetros $\beta$ (forma) e $\theta$ (escala) são análogos aos da distribuição de Weibull.- **Significado de $\beta$**: Se $\beta > 1$, a taxa de falhas aumenta com o tempo, indicando um processo de desgaste e justificando a necessidade de manutenção preventiva. Se $\beta < 1$, a confiabilidade melhora com o tempo (mortalidade infantil). Se $\beta = 1$, as falhas são constantes (aleatórias), caracterizando um Processo de Poisson Homogêneo.## Análise dos Dados### Estimação dos Parâmetros do ModeloO primeiro passo é estimar os parâmetros $\beta$ e $\theta$ a partir dos dados históricos de falha dos transformadores. Utilizamos o **Método da Máxima Verossimilhança (MLE)**, que encontra os valores dos parâmetros que maximizam a probabilidade de observar os dados coletados.```{r estimacao, eval=TRUE}dados_trafo <-data.frame(Equi =c(1, 1, 1, 2, 2, 2, 3, 3, 4, 5, 6, 7, 7, 8, 8, 9, 9, 10, 11, 12, 12, 12, 13, 14, 14, 14, 15, 16, 16, 17, 17, 18, 19, 20, 20, 20, 21, 22, 23, 23, 24, 24, 25, 26, 27, 27, 28, 28, 29, 30, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40),Failure_Time =c(8839, 17057, 21887, 9280, 16442, 21887, 10445, 13533, 7902, 8414, 13331, 17156, 21887, 16305, 21887, 16802, 21887, 4881, 16625, 7396, 7541, 19590, 2211, 15821, 19746, 19877, 1927, 15813, 21886, 15524, 21886, 21440, 369, 11664, 17031, 21857, 7544, 6039, 2168, 6698, 18840, 21879, 2288, 2499, 10668, 16838, 15550, 21887, 1616, 14041, 20004, 21888, 21888, 21888, 21888, 21888, 21888, 21888, 21888, 21888, 21888),Censura =c(1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))dados_lista <- dados_trafo %>%group_by(Equi) %>%summarise(falhas =list(Failure_Time[Censura ==1]),tempo_final =max(Failure_Time) )neg_log_likelihood_plp <-function(params, dados) { beta <- params[1] theta <- params[2]# Parâmetros devem ser positivosif (beta <=0|| theta <=0) return(Inf)# sapply itera sobre cada sistema (cada linha de 'dados_lista') total_log_likelihood <-sum(sapply(1:nrow(dados), function(i) { tempos_falha <-unlist(dados$falhas[[i]]) tempo_max <- dados$tempo_final[i]# Termo da verossimilhança para as falhas observadas# log(intensidade) = log(dweibull/pweibull) log_intensidade <-0if (length(tempos_falha) >0) { log_intensidade <-sum(dweibull(tempos_falha, shape = beta, scale = theta, log =TRUE) -pweibull(tempos_falha, shape = beta, scale = theta, lower.tail =FALSE, log.p =TRUE)) }# Termo da verossimilhança para o período total de observação# Intensidade acumulada M(T) = (T / theta)^beta intensidade_acumulada <- (tempo_max / theta)^beta# Log-verossimilhança para este sistemareturn(log_intensidade - intensidade_acumulada) }))# A função optim minimiza, então retornamos o valor negativoreturn(-total_log_likelihood)}valores_iniciais <-c(2.0, 20000)otimizacao_corrigida <-optim(par = valores_iniciais, fn = neg_log_likelihood_plp, dados = dados_lista,method ="BFGS",control =list(parscale = valores_iniciais), # Conforme final.rhessian =TRUE# Mantido para informações adicionais)beta_estimado_correto <- otimizacao_corrigida$par[1]theta_estimado_correto <- otimizacao_corrigida$par[2]```**Resultados da Estimação:**```{r print_estimacao, echo=FALSE}resultado_params <-data.frame( Parâmetro =c("Beta (Forma)", "Theta (Escala)"),`Valor Estimado`=c(beta_estimado_correto, theta_estimado_correto))knitr::kable( resultado_params,caption ="Parâmetros Estimados para o Modelo PLP via Máxima Verossimilhança.",digits =4)```O valor de $\beta$ estimado em aproximadamente `r round(beta_estimado_correto, 3)` é crucial para esta análise. Por ser significativamente maior que 1, ele confirma estatisticamente que os transformadores estão em um processo de envelhecimento e desgaste, onde a probabilidade de falha aumenta à medida que o equipamento opera por mais tempo. É este comportamento que justifica economicamente a implementação de uma política de manutenção preventiva, em vez de simplesmente reparar os equipamentos conforme eles falham.### Avaliação da Adequação do ModeloPara validar nosso modelo, comparamos a **Mean Cumulative Function (MCF)**, que representa a média de falhas acumuladas ao longo do tempo. Comparamos a curva não-paramétrica (calculada diretamente dos dados via estimador de Nelson-Aalen) com a curva paramétrica (prevista pelo nosso modelo PLP com os parâmetros estimados). Se as curvas forem próximas, o modelo é considerado adequado.```{r mcf_plot}dados_intervalo <- dados_trafo %>%group_by(Equi) %>%mutate(tempo_inicio =lag(Failure_Time, default=0), evento = Censura, tempo_fim = Failure_Time)surv_obj <-Surv(time=dados_intervalo$tempo_inicio, time2=dados_intervalo$tempo_fim, event=dados_intervalo$evento)fit_na <-survfit(surv_obj ~1)mcf_na_data <-data.frame(tempo = fit_na$time, mcf = fit_na$cumhaz)mcf_parametrica <-function(t, beta, theta) { (t / theta)^beta }ggplot(mcf_na_data, aes(x = tempo, y = mcf)) +geom_step(aes(color ="Dados (Nelson-Aalen)"), size =1.2) +stat_function(fun = mcf_parametrica, args =list(beta = beta_estimado_correto, theta = theta_estimado_correto),aes(color ="Modelo (Weibull)"), size =1) +labs(title ="Avaliação do Modelo",x ="Tempo (horas)", y ="Média de Falhas Acumuladas (MCF)", color ="Estimador") +scale_color_manual(values =c("Dados (Nelson-Aalen)"="dodgerblue", "Modelo (Weibull)"="red")) +theme_bw() +theme(legend.position ="bottom") +coord_cartesian(xlim =c(0, 22000))```A proximidade entre a curva paramétrica (vermelha), prevista por nosso modelo, e a curva não-paramétrica (azul), calculada diretamente dos dados, é notável. O modelo PLP consegue capturar com precisão a dinâmica das falhas: tanto o período inicial de baixa ocorrência quanto a aceleração do processo de desgaste após as 10.000 horas de operação. Essa validação gráfica nos dá um alto grau de confiança para usar este modelo na determinação da política ótima de manutenção.## Determinação da Política Ótima de ManutençãoO objetivo final é encontrar o intervalo de manutenção preventiva, $\tau$, que minimiza o custo esperado por unidade de tempo, $H(\tau)$. A função de custo é definida como [@gilardoni2007]:$$H(\tau) = \frac{C_{PM} + C_{MR} \cdot M(\tau)}{\tau}$$Onde $M(\tau) = (\tau / \theta)^\beta$ é a MCF do nosso modelo e a razão de custos entre reparo mínimo e manutenção preventiva é definida como $C_{MR}/C_{PM} = 15$. Para o modelo PLP, a solução ótima $\tau^*$ que minimiza $H(\tau)$ pode ser encontrada analiticamente pela fórmula:$$\tau^* = \theta \left[ \frac{C_{PM}}{C_{MR}(\beta-1)} \right]^{1/\beta}$$### Resultados da OtimizaçãoAplicando a fórmula com os parâmetros estimados, encontramos o intervalo ótimo e visualizamos a função de custo para confirmar o ponto de mínimo.```{r custo_plot}C_PM <-1; C_MR <-15tau_otimo <- theta_estimado_correto * (C_PM / ((beta_estimado_correto -1) * C_MR))^(1/ beta_estimado_correto)funcao_custo_H <-function(tau, beta, theta, cpm, cmr) { (cpm + cmr * (tau / theta)^beta) / tau }custo_otimo <-funcao_custo_H(tau_otimo, beta_estimado_correto, theta_estimado_correto, C_PM, C_MR)intervalo_tau <-seq(tau_otimo -800, tau_otimo +800, length.out =500)dados_custo <-data.frame(tau = intervalo_tau, custo =funcao_custo_H(intervalo_tau, beta_estimado_correto, theta_estimado_correto, C_PM, C_MR))ggplot(dados_custo, aes(x = tau, y = custo)) +geom_line(color ="blue", size =1) +geom_vline(xintercept = tau_otimo, color ="red", linetype ="dashed", size =1) +annotate("point", x = tau_otimo, y = custo_otimo, color ="red", size =4) +annotate("text", x = tau_otimo, y = custo_otimo + (max(dados_custo$custo, na.rm =TRUE) -min(dados_custo$custo, na.rm =TRUE))*0.05, label =paste0("Ponto Ótimo\nτ* ≈ ", round(tau_otimo), " horas\nH(τ*) ≈ ", format(custo_otimo, scientific =FALSE, digits =5)),hjust =0.5, vjust =-0.5, size=3.5 ) +labs(title ="Função de Custo da Manutenção Preventiva (Focado no Ótimo)",x ="Intervalo de Manutenção Preventiva τ (horas)",y ="Custo por Unidade de Tempo H(τ)") +theme_bw() +theme(panel.grid.major =element_line(linetype ='dotted'),panel.grid.minor =element_line(linetype ='dotted'))```**Política Ótima de Manutenção:**```{r print_otimo, echo=FALSE}resultado_politica <-data.frame( Métrica =c("Intervalo Ótimo (τ*) em Horas", "Intervalo Ótimo (τ*) em Dias"),Valor =c(round(tau_otimo), round(tau_otimo /24)))knitr::kable( resultado_politica,caption ="Política de Manutenção Preventiva Ótima.",col.names =c("Métrica", "Valor Recomendado"))```O intervalo ótimo de **`r round(tau_otimo)` horas** representa o ponto de equilíbrio econômico ideal. Realizar a manutenção antes deste marco seria prematuro e implicaria em custos desnecessários com a substituição de equipamentos ainda confiáveis. Por outro lado, adiar a manutenção para além deste ponto aumenta exponencialmente o risco de falhas em serviço, cujos custos de reparo e multas superam a economia obtida ao estender a vida útil do equipamento. Portanto, este valor não é apenas um resultado estatístico, mas uma recomendação gerencial direta.## Conclusão- O processo de falha dos transformadores é bem descrito por um modelo de **desgaste** (PLP com $\beta > 1$), validando a necessidade de manutenção preventiva.- A análise fornece uma **política de manutenção orientada por dados**: realizar uma manutenção preventiva completa a cada **`r round(tau_otimo)` horas (ou aproximadamente `r round(tau_otimo / 24)` dias)**.- Adotar este intervalo ótimo permite à empresa minimizar os custos operacionais totais, equilibrando os gastos com manutenções planejadas e os custos mais elevados de reparos corretivos não planejados.