Análise de Dados
de Processo

import pandas as pd
import numpy as np

# ── 1. Importar CSV do historian ──────────────────────────────
df = pd.read_csv(
    "historian_export.csv",
    parse_dates=["timestamp"],
    index_col="timestamp"
)
print(df.info())
print(f"Shape inicial: {df.shape}")

# ── 2. Remover registros com quality ruim ─────────────────────
df = df[df["quality"] == "GOOD"].drop(columns=["quality"])

# ── 3. Converter para numérico (coerce transforma erros em NaN)
for col in df.columns:
    df[col] = pd.to_numeric(df[col], errors="coerce")

# ── 4. Remover outliers grosseiros por faixa de engenharia ────
limites = {
    "PT101_bar" : (0,   100),   # transmissor 0–100 bar
    "TT201_C"   : (0,   400),   # termopar 0–400 °C
    "FT301_m3h" : (0,   50),    # vazão 0–50 m³/h
    "LT401_pct" : (0,   100),   # nível 0–100 %
}
for col, (lo, hi) in limites.items():
    df[col] = df[col].where(df[col].between(lo, hi))

# ── 5. Interpolação linear para NaN isolados (máx 3 min) ─────
df = df.interpolate(method="time", limit=3)

# ── 6. Reamostrar para 1 minuto (média) ──────────────────────
df = df.resample("1min").mean().round(3)

print(f"Shape final: {df.shape}")
print(df.head())

# Estatísticas descritivas completas
desc = df.describe(percentiles=[.05, .25, .5, .75, .95])
print(desc.round(3))

# Detectar assimetria e curtose (relevante para alarmes)
print("\nAssimetria (skewness):")
print(df.skew().round(3))

print("\nCurtose (kurtosis):")
print(df.kurt().round(3))

# Verificar operação fora da faixa normal (±2σ)
for col in df.columns:
    media  = df[col].mean()
    sigma  = df[col].std()
    fora   = ((df[col] < media - 2*sigma) |
              (df[col] > media + 2*sigma)).sum()
    pct    = fora / len(df) * 100
    print(f"{col}: {pct:.1f}% fora de ±2σ"
          f"  (μ={media:.2f}, σ={sigma:.3f})")

import numpy as np
from scipy import stats

col = "PT101_bar"
serie = df[col].dropna()

# ── Regressão linear (detecção de deriva) ─────────────────────
x = np.arange(len(serie))
slope, intercept, r, p, se = stats.linregress(x, serie.values)

print(f"Inclinação:  {slope*60:.4f} unid/hora")   # 60 amostras/hora
print(f"R²:          {r**2:.4f}")
print(f"p-valor:     {p:.6f} {'→ deriva significativa!' if p < 0.05 else '→ estável'}")

# ── Alerta se deriva diária superar 0.5% do span ──────────────
span       = 100    # bar (URV - LRV)
min_hora   = 60     # amostras por hora
deriva_dia = abs(slope) * min_hora * 24
limite     = span * 0.005

if deriva_dia > limite:
    print(f"⚠️  ALERTA: deriva de {deriva_dia:.3f} bar/dia"
          f" supera limite de {limite:.2f} bar/dia")
    print("   → Agendar calibração preventiva!")

# ── Média móvel (suavização de ruído) ─────────────────────────
df["PT101_MA10"] = df[col].rolling(window=10, center=True).mean()
df["PT101_MA60"] = df[col].rolling(window=60, center=True).mean()

# ── Banda de controle ±2σ (janela deslizante 1h) ──────────────
df["PT101_upper"] = df["PT101_MA60"] + 2 * df[col].rolling(60).std()
df["PT101_lower"] = df["PT101_MA60"] - 2 * df[col].rolling(60).std()

import seaborn as sns
import matplotlib.pyplot as plt

# ── Matriz de correlação de Pearson ───────────────────────────
corr_matrix = df.corr(method="pearson")
print(corr_matrix.round(3))

# ── Heatmap visual ────────────────────────────────────────────
plt.figure(figsize=(8, 6))
sns.heatmap(
    corr_matrix,
    annot=True, fmt=".2f",
    cmap="RdBu_r", vmin=-1, vmax=1,
    square=True, linewidths=0.5
)
plt.title("Correlação entre variáveis do Reator")
plt.tight_layout()
plt.savefig("correlacao_reator.png", dpi=150)

# ── Correlação cruzada com lag (causa-efeito) ─────────────────
# Ex.: mudança de FT-301 afeta PT-101 após quantos minutos?
lags = range(-30, 31)   # ±30 minutos
xcorr = [
    df["FT301_m3h"].corr(df["PT101_bar"].shift(lag))
    for lag in lags
]
lag_max = lags[xcorr.index(max(xcorr))]
print(f"Máxima correlação FT→PT em lag = {lag_max} min")
print(f"Isso indica que FT-301 afeta PT-101 com {lag_max} min de atraso")

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

def grafico_IMR(serie, titulo=""):
    """Gráfico de controle Individual e Amplitude Móvel."""
    n  = len(serie)
    MR = serie.diff().abs()   # amplitude móvel

    # Constantes d2 e D4 para n=2 (tabela Shewhart)
    d2, D4 = 1.128, 3.267

    # Limites do gráfico I
    CL_I  = serie.mean()
    sigma = MR.mean() / d2
    UCL_I = CL_I + 3 * sigma
    LCL_I = CL_I - 3 * sigma

    # Limites do gráfico MR
    CL_MR  = MR.mean()
    UCL_MR = D4 * CL_MR

    # Plot duplo I-MR
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(14, 6), sharex=True)

    # Gráfico I
    ax1.plot(serie.values, color="#22c55e", lw=1.2)
    ax1.axhline(UCL_I, color="#ef4444", ls="--", label=f"UCL={UCL_I:.2f}")
    ax1.axhline(CL_I,  color="#f59e0b", ls="-",  label=f"CL={CL_I:.2f}")
    ax1.axhline(LCL_I, color="#ef4444", ls="--", label=f"LCL={LCL_I:.2f}")

    # Destacar pontos fora de controle
    oos = (serie > UCL_I) | (serie < LCL_I)
    ax1.scatter(serie.index[oos], serie[oos], color="#ef4444", s=40, zorder=5)
    ax1.set_title(f"Gráfico I — {titulo} ({oos.sum()} pontos OOS)")
    ax1.legend(fontsize=8)

    # Gráfico MR
    ax2.bar(range(n), MR.values, color="#3b82f6", alpha=0.6)
    ax2.axhline(UCL_MR, color="#ef4444", ls="--", label=f"UCL={UCL_MR:.2f}")
    ax2.axhline(CL_MR,  color="#f59e0b", ls="-",  label=f"CL={CL_MR:.2f}")
    ax2.set_title("Gráfico MR — Amplitude Móvel")
    ax2.legend(fontsize=8)

    plt.tight_layout()
    return fig, sigma, oos.sum()

fig, sigma, n_oos = grafico_IMR(df["PT101_bar"], "PT-101 Pressão")
print(f"σ estimado: {sigma:.4f} bar | Pontos OOS: {n_oos}")
plt.savefig("spc_PT101.png", dpi=150)

import numpy as np
from scipy.fft import fft, fftfreq
import matplotlib.pyplot as plt

# Sinal de vibração de uma bomba (acelerômetro, 1000 Hz)
fs     = 1000                        # frequência de amostragem (Hz)
signal = df_vib["accel_g"].values   # aceleração em g
N      = len(signal)

# ── FFT ───────────────────────────────────────────────────────
yf   = fft(signal)
xf   = fftfreq(N, 1/fs)[:N//2]     # frequências positivas
amp  = (2.0/N) * np.abs(yf[:N//2])  # amplitude real

# ── Frequências características da bomba ─────────────────────
RPM        = 1450           # rotação nominal
f_rotor    = RPM / 60       # = 24,2 Hz → 1× rotação
f_palhetas = f_rotor * 6    # = 145 Hz → 6 palhetas (BPF)
f_rolamento= 87.3           # frequência de falha BPFO (calculada)

# ── Extrair amplitude nas frequências de interesse ────────────
def amp_em(freq, tol=2.0):
    mask = np.abs(xf - freq) < tol
    return amp[mask].max() if mask.any() else 0

print(f"1× rotação ({f_rotor:.1f} Hz):    {amp_em(f_rotor):.4f} g")
print(f"BPF        ({f_palhetas:.1f} Hz): {amp_em(f_palhetas):.4f} g")
print(f"BPFO       ({f_rolamento:.1f} Hz): {amp_em(f_rolamento):.4f} g")

# Alarme: amplitude 1× > 1g indica desbalanceamento
if amp_em(f_rotor) > 1.0:
    print("⚠️  ALERTA: Desbalanceamento! → Verificar balanceamento dinâmico")

from sklearn.ensemble import GradientBoostingRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_absolute_error, r2_score
from sklearn.preprocessing import StandardScaler
import joblib

# ── Sensor virtual: prever PT-101 a partir de FT, TT e LT ────
# Útil quando PT-101 está em manutenção ou indica erro

features = ["TT201_C", "FT301_m3h", "LT401_pct"]
target   = "PT101_bar"

X = df[features].dropna()
y = df.loc[X.index, target]

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, shuffle=False  # não embaralhar séries temporais!
)

# Normalização
scaler  = StandardScaler()
X_train = scaler.fit_transform(X_train)
X_test  = scaler.transform(X_test)

# ── Modelo GBR ───────────────────────────────────────────────
modelo = GradientBoostingRegressor(
    n_estimators=200, max_depth=4,
    learning_rate=0.05, random_state=42
)
modelo.fit(X_train, y_train)

y_pred = modelo.predict(X_test)
mae    = mean_absolute_error(y_test, y_pred)
r2     = r2_score(y_test, y_pred)

print(f"Sensor virtual PT-101:")
print(f"  MAE = {mae:.4f} bar")
print(f"  R²  = {r2:.4f} ({'excelente' if r2>0.95 else 'bom' if r2>0.85 else 'revisar'})")

# Detectar discrepância: instrumento real vs sensor virtual
diferenca = abs(df[target] - modelo.predict(scaler.transform(df[features])))
if diferenca.iloc[-1] > 2.0:  # 2 bar de diferença
    print("⚠️  Divergência entre PT-101 real e virtual → Verificar instrumento!")

# Salvar modelo treinado
joblib.dump((modelo, scaler), "sensor_virtual_PT101.pkl")

import pandas as pd
import smtplib
from email.mime.multipart import MIMEMultipart
from email.mime.text import MIMEText
from email.mime.base import MIMEBase
from email import encoders
from datetime import date
import openpyxl

def gerar_relatorio_excel(df, anomalias, arquivo):
    """Gera relatório diário em Excel com múltiplas abas."""
    with pd.ExcelWriter(arquivo, engine="openpyxl") as writer:
        # Aba 1 — Resumo estatístico
        df.describe().round(3).to_excel(writer, sheet_name="Resumo")

        # Aba 2 — Dados brutos
        df.to_excel(writer, sheet_name="Dados")

        # Aba 3 — Anomalias detectadas
        anomalias.to_excel(writer, sheet_name="Anomalias")

    print(f"✅ Relatório salvo: {arquivo}")

def enviar_alerta(destinatario, assunto, corpo, anexo=None):
    """Envia alerta por e-mail com relatório em anexo."""
    msg = MIMEMultipart()
    msg["From"]    = "planta@empresa.com.br"
    msg["To"]      = destinatario
    msg["Subject"] = assunto
    msg.attach(MIMEText(corpo, "html"))

    if anexo:
        with open(anexo, "rb") as f:
            part = MIMEBase("application", "octet-stream")
            part.set_payload(f.read())
            encoders.encode_base64(part)
            part.add_header("Content-Disposition",
                            f"attachment; filename={anexo}")
            msg.attach(part)

    with smtplib.SMTP("mail.empresa.com.br", 587) as smtp:
        smtp.starttls()
        smtp.login("planta@empresa.com.br", "senha")
        smtp.send_message(msg)
    print(f"📧 Alerta enviado para {destinatario}")

# ── Execução diária (agendar com cron ou Task Scheduler) ──────
hoje     = date.today().isoformat()
arquivo  = f"relatorio_{hoje}.xlsx"

gerar_relatorio_excel(df, anomalias, arquivo)

if len(anomalias) > 0:
    corpo = f"""
    <h2>⚠️ Relatório de Processo — {hoje}</h2>
    <p><b>{len(anomalias)} anomalias</b> detectadas nas últimas 24h.</p>
    <p>Ver relatório em anexo para detalhes.</p>
    """
    enviar_alerta(
        "instrumentacao@empresa.com.br",
        f"⚠️ Anomalias Processo — {hoje}",
        corpo, arquivo
    )