ML Anomaly Detection for HTA/BT Transformers

Detect transformer faults 30-45 min early using Isolation Forest and Autoencoder.

Introduction

30 Second Summary

Every time you flip a light switch, electricity flows through a transformer sitting quietly in your neighborhood. When that transformer degrades without warning, entire streets lose power. Detecting problems early means fewer outages and safer equipment.

In this project, you will build a predictive monitoring system for a distribution transformer (20kV/400V, 630 kVA). You will train Machine Learning models to detect degradation 30 to 45 minutes before traditional protection systems react.

What You'll Build

You run your notebook, inject a fault scenario into a simulated transformer, and watch your ML models flag the anomaly while the legacy DGPT2 protection system stays silent.

By the end of this project, you'll have:

  • A synthetic data generator that simulates one month of sensor readings (currents, oil temperature, ambient temperature, dissolved H2) from a 630 kVA transformer, following the CEI 60076-7 thermal model.
  • An Isolation Forest and an Autoencoder trained only on normal operating data, capable of detecting four real-world fault scenarios before static thresholds trigger.
  • A composite Health Index that synthesizes temperature, load, and anomaly scores into a single value comparable across transformers.
  • Secret Mission: Simulate a fleet of 5 transformers and build a maintenance prioritization dashboard ranking them by Health Index, just like an Enedis exploitation engineer would.

Are there any prerequisites?

This project assumes you have Python installed and a working Jupyter Notebook environment. You should also be comfortable with basic ML concepts (training, scoring, unsupervised learning) and have foundational knowledge of electrical distribution networks.

Before We Start

Avant de plonger dans le code, prenons un moment pour cadrer ce que vous allez construire et pourquoi ce projet a un impact direct sur la maintenance du reseau de distribution electrique.

Set Up Your Environment

Le modele thermique CEI 60076-7 et les algorithmes de detection d'anomalies reposent sur des versions specifiques de librairies Python. Sans les bonnes versions, vous risquez des erreurs d'import silencieuses ou des resultats incorrects.

Dans cette etape, vous allez confirmer votre environnement Python, installer les dependances, creer votre notebook Jupyter, et valider que tout fonctionne ensemble avec un test rapide d'Isolation Forest.

In this step, get ready to:
  • Confirmer Python 3.10+ et installer toutes les dependances (numpy, pandas, scikit-learn, TensorFlow, matplotlib).
  • Creer le notebook du projet et executer tous les imports sans erreur.
  • Valider le setup avec un test rapide d'IsolationForest sur des donnees aleatoires.
Verify your Python environment and install dependencies

Python 3.10+ est deja installe sur votre machine. Confirmez la version avant d'installer les librairies ML.

  • Verifiez votre version de Python en executant cette commande dans votre terminal :
python --version

Vous devriez voir une version 3.10 ou superieure (par exemple Python 3.11.9).

  • Installez les dependances du projet en executant cette commande :
pip install numpy pandas scikit-learn tensorflow matplotlib

What does this command do?

Cette commande installe (ou met a jour) cinq librairies en une seule fois.

  • numpy et pandas pour la manipulation de donnees numeriques et de series temporelles.
  • scikit-learn fournit l'IsolationForest et les outils de preprocessing.
  • tensorflow contient Keras pour construire l'Autoencoder.
  • matplotlib sert aux visualisations metier (courbes de charge, zones d'alerte).

L'installation peut prendre 1 a 3 minutes selon votre connexion. Vous verrez une serie de messages Successfully installed... ou Requirement already satisfied a la fin.

Seeing an error during pip install?

  • Si vous voyez une erreur de permissions, relancez avec pip install --user suivi des memes packages.
  • Si TensorFlow echoue avec un message de compatibilite, verifiez que votre Python est bien 3.10 ou 3.11 (TensorFlow ne supporte pas toutes les versions mineures).
  • help me debug my pip install error for tensorflow
Create the notebook and import all libraries

Toutes les librairies sont installees. Creez maintenant le notebook qui contiendra l'ensemble du pipeline de detection.

  • Lancez Jupyter en executant cette commande :
jupyter notebook

Votre navigateur s'ouvre sur l'interface Jupyter. Si le navigateur ne s'ouvre pas automatiquement, copiez l'URL affichee dans le terminal.

  • Cliquez sur New en haut a droite, puis selectionnez Python 3.
  • Renommez le notebook en cliquant sur Untitled en haut de la page, puis tapez anomaly_detection_hta_bt et confirmez.
  • Dans la premiere cellule du notebook, collez le code suivant :
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.ensemble import IsolationForest
from sklearn.preprocessing import StandardScaler
from tensorflow.keras.models import Model as KerasModel
from tensorflow.keras.layers import Input, Dense
from tensorflow.keras.optimizers import Adam
  • Executez la cellule avec Shift+Enter.

La cellule doit s'executer sans aucune erreur. Vous verrez un numero [1] apparaitre a gauche de la cellule.

What do these imports do?

  • IsolationForest est l'algorithme de detection d'anomalies non-supervise que vous entrainerez sur les donnees normales du transformateur.
  • StandardScaler normalise les features pour que chaque variable ait le meme poids dans le modele.
  • KerasModel, Input, Dense sont les briques pour construire l'Autoencoder (reseau de neurones en sablier).
  • Adam est l'optimiseur qui ajustera les poids du reseau pendant l'entrainement.

Import error on tensorflow.keras?

  • Si vous voyez ModuleNotFoundError: No module named 'tensorflow', relancez pip install tensorflow dans un terminal et redemarrez le kernel Jupyter (menu Kernel puis Restart).
  • Si l'import est lent (30+ secondes), c'est normal au premier chargement de TensorFlow.
  • help me fix my TensorFlow import error in Jupyter
Validate the setup with a quick IsolationForest test

Les imports passent, mais est-ce que les algorithmes fonctionnent vraiment ? Un test rapide sur des donnees aleatoires confirme que l'ensemble du pipeline est operationnel.

  • Dans une nouvelle cellule (cellule 2), collez ce code de validation :
# Quick validation: IsolationForest on random data
rng = np.random.default_rng(42)
X_test = rng.normal(size=(100, 3))

# Train on normal data
clf = IsolationForest(n_estimators=50, random_state=42)
clf.fit(X_test)
predictions = clf.predict(X_test)

# Confirm fit and predict work
print(f"Shape des predictions: {predictions.shape}")
print(f"Labels uniques: {np.unique(predictions)}")
print("Setup valide - IsolationForest OK!")
  • Executez la cellule avec Shift+Enter.

Vous devriez voir trois lignes de sortie confirmant que les predictions ont une shape (100,) et que les labels sont [-1 1] (anomalies et points normaux).

What does this test prove?

Ce test confirme trois choses en un seul run.

  • numpy genere correctement des donnees aleatoires (la matrice 100x3).
  • fit() entraine l'arbre d'isolation sans erreur sur ces donnees.
  • predict() retourne des labels (-1 pour anomalie, 1 pour normal) de la bonne dimension.

Seeing a NameError on np?

  • Verifiez que vous avez bien execute la cellule 1 (imports) avant la cellule 2.
  • Si le kernel a redemarre, re-executez toutes les cellules avec le menu Cell puis Run All.
  • help me fix NameError in my Jupyter notebook

A vous de predire : quand vous lancerez la commande de verification des versions dans le terminal, quels numeros de version vous attendez-vous a voir pour scikit-learn et TensorFlow ?

  • Ouvrez un nouveau terminal (sans fermer Jupyter) et executez cette commande :
python -c "import sklearn; import tensorflow; print(sklearn.__version__, tensorflow.__version__)"

Vous devriez voir deux numeros de version. Confirmez que scikit-learn affiche 1.5 ou superieur et que TensorFlow affiche 2.15 ou superieur.

✔️ Awesome, I've got everything!

Parfait. Verifiez que votre notebook anomaly_detection_hta_bt.ipynb est bien sauvegarde (Ctrl+S ou Cmd+S).

ⓧ I'd like to double check the full code

Voici le contenu complet de vos deux cellules :

Cellule 1 (imports) :

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.ensemble import IsolationForest
from sklearn.preprocessing import StandardScaler
from tensorflow.keras.models import Model as KerasModel
from tensorflow.keras.layers import Input, Dense
from tensorflow.keras.optimizers import Adam

Cellule 2 (test rapide) :

# Quick validation: IsolationForest on random data
rng = np.random.default_rng(42)
X_test = rng.normal(size=(100, 3))

# Train on normal data
clf = IsolationForest(n_estimators=50, random_state=42)
clf.fit(X_test)
predictions = clf.predict(X_test)

# Confirm fit and predict work
print(f"Shape des predictions: {predictions.shape}")
print(f"Labels uniques: {np.unique(predictions)}")
print("Setup valide - IsolationForest OK!")

Votre environnement est pret et valide. Dans la prochaine etape, vous allez generer un mois complet de donnees capteurs en utilisant le modele thermique CEI 60076-7 pour simuler un transformateur 630 kVA en fonctionnement reel.

Generate Sensor Data Using the IEC 60076-7 Thermal Model

Your environment is ready and your libraries are loaded. Now you need data to train your anomaly detection models.

Real-time monitoring data from Enedis transformers (collected by GridWatch or Sentinel FRTU sensors) is not publicly available. To build a meaningful ML pipeline, you need synthetic time series that respect the actual thermal physics of a transformer. Otherwise, your model would learn artificial patterns with no operational value.

In this step, you will implement the IEC 60076-7 thermal model (difference equations method from Clause 8.2.3) and generate one month of realistic sensor measurements for a 630 kVA ONAN distribution transformer.

In this step, get ready to:
  • Implement the simplified CEI 60076-7 thermal model to calculate oil and hot-spot temperatures.
  • Generate a realistic load profile, ambient temperature, three-phase currents, and dissolved H2 concentration.
  • Assemble a full DataFrame of 30 days at 10-minute intervals and visualize the time series.
Implement the IEC 60076-7 thermal model

The IEC 60076-7 standard describes how a transformer's temperature evolves over time as a function of load and ambient conditions. The difference equations method (Clause 8.2.3) discretizes the thermal differential equations into time steps, making it possible to compute the oil temperature and hot-spot temperature iteratively using Euler's method.

You will start by defining the ONAN thermal constants for a 630 kVA distribution transformer, then implement the simulation function.

  • In your notebook, create a new cell below your import cells and define the ONAN thermal constants by pasting this code:
# Constantes ONAN pour transformateur de distribution 630 kVA
# (valeurs IEC 60076-7 Table 4, refroidissement ONAN)
DELTA_THETA_OR = 52    # Echauffement huile nominal (K)
DELTA_THETA_HR = 26    # Gradient hot-spot nominal (K)
TAU_0 = 150            # Constante de temps huile (min)
TAU_W = 7              # Constante de temps enroulement (min)
R = 6                  # Ratio pertes charge / pertes a vide
X_OIL = 0.8            # Exposant huile (ONAN)
Y_WINDING = 1.3        # Exposant enroulement (ONAN)
K11 = 0.5              # Constante modele thermique
K21 = 2.0              # Constante modele thermique
K22 = 2.0              # Constante modele thermique

What do these constants represent?

These constants define the thermal behavior of an oil-immersed transformer with natural convection cooling (ONAN). DELTA_THETA_OR is the steady-state oil temperature rise at full load. TAU_0 is the oil thermal time constant (how fast the oil heats up).

The ratio R = 6 means that load losses at full load are 6 times the no-load losses. The exponents X_OIL and Y_WINDING model the non-linear relationship between load and temperature rise.

  • In the next cell, implement the thermal simulation function by pasting this code:
def simulate_thermal_model(load_profile, ambient_temp, dt=10, tau_0_override=None):
    """Simule le modele thermique CEI 60076-7 (equations aux differences)."""
    tau_0 = tau_0_override if tau_0_override is not None else TAU_0
    n = len(load_profile)
    # Conditions initiales (regime permanent au premier pas)
    K_0 = load_profile[0]
    theta_a_0 = ambient_temp[0]
    theta_0 = ((1 + K_0**2 * R) / (1 + R))**X_OIL * DELTA_THETA_OR + theta_a_0
    delta_h1 = K21 * K_0**Y_WINDING * DELTA_THETA_HR
    delta_h2 = (K21 - 1) * K_0**Y_WINDING * DELTA_THETA_HR
    theta_oil = np.zeros(n)
    theta_hs = np.zeros(n)
    for i in range(n):
        K_i = load_profile[i]
        theta_a_i = ambient_temp[i]
        # Mise a jour temperature huile (top-oil rise)
        d_theta_0 = (dt / (K11 * tau_0)) * (
            ((1 + K_i**2 * R) / (1 + R))**X_OIL * DELTA_THETA_OR
            - (theta_0 - theta_a_i))
        theta_0 = theta_0 + d_theta_0
        # Mise a jour gradient hot-spot (composante rapide)
        d_h1 = (dt / (K22 * TAU_W)) * (K21 * DELTA_THETA_HR * K_i**Y_WINDING - delta_h1)
        delta_h1 = delta_h1 + d_h1
        # Mise a jour gradient hot-spot (composante lente - flux huile)
        d_h2 = (dt / ((1 / K22) * tau_0)) * (
            (K21 - 1) * DELTA_THETA_HR * K_i**Y_WINDING - delta_h2)
        delta_h2 = delta_h2 + d_h2
        theta_oil[i] = theta_0
        theta_hs[i] = theta_0 + (delta_h1 - delta_h2)
    return theta_oil, theta_hs

What does this function do?

At each 10-minute time step, the function updates three quantities using Euler's method:

  • The top-oil temperature rise (d_theta_0), which depends on the current load factor K and how far the oil is from its steady-state target.
  • The fast hot-spot gradient (d_h1), driven by the winding heating directly.
  • The slow hot-spot gradient (d_h2), driven by oil flow changes in the cooling circuit.

The tau_0_override parameter lets you simulate degraded cooling (clogged radiators) later in the project.

  • Run both cells (constants and function) by pressing Shift+Enter on each.
  • In a new cell, run a quick sanity check by pasting and executing this code:
# Test rapide du modele thermique
test_load = np.full(144, 0.6)  # 1 jour a charge constante 60%
test_ambient = np.full(144, 10.0)  # 10 degC constant
theta_oil, theta_hs = simulate_thermal_model(test_load, test_ambient, dt=10)
print(f"Theta_oil final: {theta_oil[-1]:.1f} degC")
print(f"Theta_hs final: {theta_hs[-1]:.1f} degC")

You should see Theta_oil final around 30-35 degrees C and Theta_hs final a few degrees higher. These values confirm the model is computing realistic temperatures for a lightly loaded transformer in winter.

Seeing a NameError?

Make sure you ran the constants cell before the function cell. The function references DELTA_THETA_OR, TAU_0, and the other constants directly.

If you see np is not defined, run your import cell (Cell 1) again first.

Why is my thermal model returning a NameError?

Generate a realistic load profile and auxiliary measurements

A distribution transformer in a residential area shows a specific daily pattern. The evening peak (18h-21h) reaches 70-90% load factor as households cook and heat. The night trough (02h-06h) drops to 20-30%. Weekends have lower overall consumption.

  • In a new cell, create the load profile generator by pasting this code:
def generate_load_profile(n_days=30, dt_min=10, seed=42):
    """Genere un profil de charge realiste pour zone pavillonnaire."""
    rng = np.random.default_rng(seed)
    n_points = n_days * 24 * 60 // dt_min
    timestamps = pd.date_range("2025-01-01", periods=n_points, freq=f"{dt_min}min")
    hours = timestamps.hour + timestamps.minute / 60
    day_of_week = timestamps.dayofweek
    # Profil journalier type (creux nuit, pointe soir)
    base_load = np.where(
        (hours >= 2) & (hours < 6), 0.25,
        np.where((hours >= 6) & (hours < 9), 0.5,
            np.where((hours >= 9) & (hours < 17), 0.45,
                np.where((hours >= 17) & (hours < 21), 0.75,
                    np.where((hours >= 21) & (hours < 23), 0.55, 0.35)))))
    # Reduction weekend (consommation plus faible)
    weekend_factor = np.where(day_of_week >= 5, 0.85, 1.0)
    load = base_load * weekend_factor
    # Bruit aleatoire pour realisme
    load = load + rng.normal(0, 0.03, n_points)
    load = np.clip(load, 0.1, 1.0)
    return timestamps, load

What does this load profile look like?

The nested np.where calls create a piecewise daily curve. The peak value of 0.75 (75% of rated power) happens between 17h and 21h. The weekend factor of 0.85 reduces all values by 15% on Saturdays and Sundays.

The Gaussian noise (standard deviation 0.03) adds realistic measurement jitter. The seed=42 ensures reproducibility across runs.

  • Run the cell, then verify it works in a new cell:
# Verifier le profil de charge
timestamps, load = generate_load_profile(n_days=30, dt_min=10, seed=42)
print(f"Points generes: {len(load)}")
print(f"Charge min: {load.min():.3f}, max: {load.max():.3f}, moyenne: {load.mean():.3f}")

You should see 4320 points generated, with a minimum around 0.1, maximum around 0.85, and mean around 0.45.

  • Now add the three auxiliary measurement generators in a new cell:
def generate_ambient_temp(timestamps, seed=42):
    """Genere temperature ambiante avec cycle jour/nuit (hiver)."""
    rng = np.random.default_rng(seed)
    hours = timestamps.hour + timestamps.minute / 60
    day_of_year = timestamps.dayofyear
    # Cycle saisonnier (hiver ~ 5 degC moyen)
    seasonal = 5 + 3 * np.sin(2 * np.pi * (day_of_year - 80) / 365)
    # Cycle jour/nuit (+/- 4 degC)
    diurnal = 4 * np.sin(2 * np.pi * (hours - 6) / 24)
    noise = rng.normal(0, 1, len(timestamps))
    return seasonal + diurnal + noise


def generate_currents(load_profile, I_nominal=909, desequilibre_normal=0.08, seed=42):
    """Genere courants triphases avec desequilibre normal."""
    rng = np.random.default_rng(seed)
    n = len(load_profile)
    I_mean = load_profile * I_nominal
    I1 = I_mean * (1 + rng.uniform(-desequilibre_normal/2, desequilibre_normal/2, n))
    I2 = I_mean * (1 + rng.uniform(-desequilibre_normal/2, desequilibre_normal/2, n))
    I3 = I_mean * (1 + rng.uniform(-desequilibre_normal/2, desequilibre_normal/2, n))
    return I1, I2, I3


def generate_h2(n_points, base_ppm=7, noise_std=2, seed=42):
    """Genere concentration H2 dissoute normale."""
    rng = np.random.default_rng(seed)
    return np.clip(base_ppm + rng.normal(0, noise_std, n_points), 0, None)

Why these specific values?

The nominal current I_nominal=909 A corresponds to a 630 kVA transformer on a 400V secondary (630000 / (400 * sqrt(3)) = 909 A). The 8% normal imbalance reflects typical residential phase distribution.

The H2 baseline of 7 ppm with noise represents healthy oil. Dissolved hydrogen is produced by internal defects. A healthy transformer shows 5-10 ppm with random fluctuation.

  • Run this cell to define the three functions.

Seeing an IndentationError?

Make sure the blank lines between the three functions are truly empty (no trailing spaces). Python is sensitive to inconsistent indentation inside function bodies.

Help me fix my IndentationError in the generator functions.

Assemble the DataFrame and visualize normal operation

With all the individual generators ready, you can now assemble them into a single function that produces the complete sensor dataset for one month of normal operation.

  • In a new cell, create the assembly function by pasting this code:
def build_normal_dataset(n_days=30, dt_min=10, seed=42):
    """Construit le dataset complet en fonctionnement normal."""
    timestamps, load = generate_load_profile(n_days, dt_min, seed)
    ambient = generate_ambient_temp(timestamps, seed)
    theta_oil, theta_hs = simulate_thermal_model(load, ambient, dt_min)
    I1, I2, I3 = generate_currents(load, seed=seed)
    h2 = generate_h2(len(timestamps), seed=seed)
    # Calcul du desequilibre en pourcentage
    I_mean = (I1 + I2 + I3) / 3
    desequilibre = np.max(
        np.abs(np.column_stack([I1, I2, I3]) - I_mean.reshape(-1, 1)), axis=1
    ) / (I_mean + 1e-6) * 100
    df = pd.DataFrame({
        "timestamp": timestamps, "K": load,
        "I1": I1, "I2": I2, "I3": I3,
        "theta_oil": theta_oil, "theta_hs": theta_hs,
        "theta_amb": ambient, "H2_ppm": h2,
        "desequilibre_pct": desequilibre
    })
    return df

What does build_normal_dataset assemble?

This function chains all your generators together. It produces a DataFrame with 10 columns representing exactly what a GridWatch or Sentinel FRTU sensor would measure every 10 minutes: load factor, three-phase currents, oil temperature, hot-spot temperature, ambient temperature, dissolved H2, and phase imbalance percentage.

  • Run the cell, then generate your dataset and inspect it in a new cell:
# Generer les donnees normales (1 mois)
df_normal = build_normal_dataset(n_days=30, dt_min=10, seed=42)
print(f"Shape: {df_normal.shape}")
print(f"Theta_oil max: {df_normal['theta_oil'].max():.1f} degC")
print(f"Theta_oil min: {df_normal['theta_oil'].min():.1f} degC")
df_normal.head()

You should see a shape of (4320, 10). That is 30 days multiplied by 144 measurements per day (one every 10 minutes). The maximum oil temperature should be between 65 and 75 degrees C.

Shape is wrong or theta_oil looks unrealistic?

If your shape shows fewer rows, check that n_days=30 is passed correctly. If theta_oil exceeds 90 degrees C in normal mode, double-check that your ONAN constants cell was run with the correct values (DELTA_THETA_OR = 52, not 72).

My DataFrame shape or temperature values look wrong. Help me debug.

Now visualize the time series to confirm the data follows expected physical patterns. Before you run the plot, think about what you expect to see. What shape should the oil temperature curve take over 30 days if the transformer follows a residential load cycle?

  • In a new cell, create the visualization by pasting and running this code:
# Visualisation des series temporelles normales
fig, axes = plt.subplots(3, 1, figsize=(14, 8), sharex=True)
# Temperature huile
axes[0].plot(df_normal["timestamp"], df_normal["theta_oil"], color="orangered", linewidth=0.5)
axes[0].set_ylabel("Theta_oil (degC)")
axes[0].set_title("Temperature huile - 1 mois de fonctionnement normal")
axes[0].axhline(y=95, color="orange", linestyle="--", alpha=0.7, label="Seuil alarme DGPT2")
axes[0].legend()
# Taux de charge
axes[1].plot(df_normal["timestamp"], df_normal["K"], color="steelblue", linewidth=0.5)
axes[1].set_ylabel("K (taux de charge)")
axes[1].set_title("Profil de charge")
# H2 dissous
axes[2].plot(df_normal["timestamp"], df_normal["H2_ppm"], color="green", linewidth=0.5)
axes[2].set_ylabel("H2 (ppm)")
axes[2].set_title("Hydrogene dissous")
axes[2].set_xlabel("Date")
plt.tight_layout()
plt.show()

What should I see?

The oil temperature curve should show clear day/night oscillations, with peaks in the evening (when load is highest) and troughs at night. The maximum oil temperature stays well below the DGPT2 alarm threshold of 95 degrees C (the orange dashed line).

The load profile should show daily evening peaks between 17h and 21h, with visible weekend dips. The H2 plot should look like random noise centered around 7 ppm.

  • Confirm that your theta_oil curve follows the day/night cycles and evening load peaks (19h-21h), with a maximum between 65 and 75 degrees C.

Plot is empty or shows a flat line?

If the plot is empty, make sure plt.show() is included at the end. In some Jupyter configurations, you need %matplotlib inline at the top of your notebook.

If theta_oil shows a flat line, the thermal model loop may not be iterating. Verify that your generate_load_profile returns varying values (not all zeros).

My plot is empty or the temperature curve looks flat. Help me debug.

✔️ Awesome, I've got everything!

Great. Double check you have saved your notebook with all cells executed in order.

ⓧ I'd like to double check the full code

Scroll through your notebook and confirm you have these cells in order:

  • Cell 1: All imports (numpy, pandas, matplotlib, sklearn, tensorflow).
  • Cell 2: IsolationForest quick-test (from Step 1).
  • Cell 3: ONAN constants (12 constants from DELTA_THETA_OR through K22).
  • Cell 4: simulate_thermal_model function definition.
  • Cell 5: Quick thermal test (printing theta_oil and theta_hs final values).
  • Cell 6: generate_load_profile function definition.
  • Cell 7: Load profile verification (printing 4320 points, min, max, mean).
  • Cell 8: generate_ambient_temp, generate_currents, and generate_h2 function definitions.
  • Cell 9: build_normal_dataset function definition.
  • Cell 10: df_normal = build_normal_dataset(...) with shape and theta_oil printout.
  • Cell 11: 3-panel matplotlib visualization.

You now have a physically realistic synthetic dataset that mirrors what a real GridWatch sensor would capture on an Enedis distribution transformer. Next up, you will apply the DGPT2 static threshold detection to this data and discover its limitations on progressive faults.

Apply Static Threshold Detection (DGPT2 Approach)

Your synthetic dataset is generating realistic transformer behavior. Now it is time to understand what already protects these transformers in the field before adding any ML.

The DGPT2 (Detecteur de Gaz, Pression et Temperature a 2 seuils) is the standard protection device on oil-filled transformers in France. It triggers when a fixed threshold is crossed. In this step, you will implement that approach, inject realistic fault scenarios, and observe its blind spots on progressive degradations.

In this step, get ready to:
  • Define the DGPT2 static thresholds and code a detection function.
  • Create 4 anomaly injection functions covering real-world fault scenarios.
  • Observe how static thresholds fail on progressive degradations.
Define the DGPT2 thresholds and detection function

The DGPT2 uses fixed temperature and gas concentration thresholds. When oil temperature exceeds 95 degrees C, an alarm is raised. At 105 degrees C, the transformer is disconnected from the HTA network. For dissolved hydrogen, 100 ppm triggers an alarm.

  • In a new cell in your notebook, define the DGPT2 constants and detection function by pasting this code:
# =============================================================================
# DETECTION PAR SEUILS DGPT2
# =============================================================================

SEUIL_ALARME_TEMP = 95          # degC - seuil alarme (signalisation)
SEUIL_DECLENCHEMENT_TEMP = 105  # degC - seuil declenchement (coupure HTA)
SEUIL_H2_ALARME = 100           # ppm - seuil alarme gaz dissous

def detect_dgpt2(df):
    """Detecte les franchissements de seuils DGPT2."""
    alarme_temp = df["theta_oil"] > SEUIL_ALARME_TEMP
    declenchement_temp = df["theta_oil"] > SEUIL_DECLENCHEMENT_TEMP
    alarme_h2 = df["H2_ppm"] > SEUIL_H2_ALARME

    return pd.DataFrame({
        "alarme_temp": alarme_temp,
        "declenchement_temp": declenchement_temp,
        "alarme_h2": alarme_h2
    }, index=df.index)

What does this code do?

The three constants mirror the real DGPT2 protection logic: a first threshold raises an alarm (the operator is notified), and a second threshold forces disconnection (the transformer is de-energized to prevent catastrophic failure).

The detect_dgpt2 function simply compares each row against these fixed values and returns a boolean DataFrame. This is exactly how static protection works: no memory, no trend analysis, just "is the value above the line right now?"

  • Run the cell, then verify the function works on your normal data by running this in a new cell:
# Verifier sur donnees normales : aucune alarme attendue
dgpt2_normal = detect_dgpt2(df_normal)
print(f"Alarmes temperature sur donnees normales: {dgpt2_normal['alarme_temp'].sum()}")
print(f"Declenchements sur donnees normales: {dgpt2_normal['declenchement_temp'].sum()}")
print(f"Alarmes H2 sur donnees normales: {dgpt2_normal['alarme_h2'].sum()}")

You should see 0 alarms across all three columns. Your normal data stays well within safe operating limits.

Seeing non-zero alarm counts?

If you see alarms on normal data, your thermal simulation from the previous step may have diverged. Check that df_normal["theta_oil"].max() is below 80 degrees C. If it is higher, re-run the build_normal_dataset cell with seed=42 to ensure reproducibility.

Help me debug non-zero alarms on normal data

Inject the 4 fault scenarios

Real transformer failures are rarely instantaneous. They develop over minutes to hours. You will now create 4 injection functions that simulate realistic progressive faults on top of normal operating data.

  • In a new cell, define the surcharge (overload) injection function:
def inject_surcharge(df, start_idx, duration_steps=12):
    """Surcharge progressive : K passe de 0.8 a 1.3 sur duration_steps."""
    df_anom = df.copy()
    end_idx = start_idx + duration_steps
    ramp = np.linspace(0, 0.5, duration_steps)
    df_anom.loc[start_idx:end_idx-1, "K"] = df_anom.loc[start_idx:end_idx-1, "K"].values + ramp

    # Recalculer thermique
    theta_oil, theta_hs = simulate_thermal_model(
        df_anom["K"].values, df_anom["theta_amb"].values, dt=10
    )
    df_anom["theta_oil"] = theta_oil
    df_anom["theta_hs"] = theta_hs
    return df_anom

What does this code do?

This simulates a progressive overload. The load factor K ramps from its current value up by 0.5 pu over 12 time steps (2 hours at 10-minute intervals). The thermal model is fully recalculated so the temperature response is physically consistent with the increased load.

  • Run the cell, then in a new cell define the phase imbalance injection:
def inject_desequilibre(df, start_idx, duration_steps=18, deviation=0.30):
    """Desequilibre de phases : I1 diverge de deviation (30%)."""
    df_anom = df.copy()
    end_idx = start_idx + duration_steps
    ramp = np.linspace(0, deviation, duration_steps)
    I_nominal = 909
    df_anom.loc[start_idx:end_idx-1, "I1"] = (
        df_anom.loc[start_idx:end_idx-1, "I1"].values
        + ramp * I_nominal * df_anom.loc[start_idx:end_idx-1, "K"].values
    )
    # Recalculer desequilibre
    I_cols = df_anom[["I1", "I2", "I3"]].values[start_idx:end_idx]
    I_mean_local = I_cols.mean(axis=1)
    deseq = np.max(np.abs(I_cols - I_mean_local.reshape(-1, 1)), axis=1) / (I_mean_local + 1e-6) * 100
    df_anom.loc[start_idx:end_idx-1, "desequilibre_pct"] = deseq
    return df_anom

What does this code do?

This simulates an abusive single-phase connection (like a large industrial client overloading one phase). Phase I1 progressively diverges by 30% over 18 steps (3 hours). The percent imbalance metric is recalculated to reflect this drift.

  • Run the cell. Verify the surcharge injection works by running this quick test:
# Test rapide de l'injection de surcharge
df_test_surcharge = inject_surcharge(df_normal, start_idx=2000, duration_steps=12)
print(f"Theta_oil max apres surcharge: {df_test_surcharge['theta_oil'].max():.1f} degC")
print(f"K max pendant surcharge: {df_test_surcharge['K'].max():.2f}")

You should see a theta_oil max above 90 degrees C and a K max around 1.3. The overload is pushing the transformer toward its thermal limits.

  • In a new cell, define the air-cooler clogging injection:
def inject_colmatage(df, start_idx, duration_steps=24):
    """Colmatage aero-refrigerant : tau_0 augmente (refroidissement degrade)."""
    df_anom = df.copy()
    # Simuler avec tau_0 croissant
    tau_0_values = np.full(len(df), TAU_0, dtype=float)
    ramp = np.linspace(TAU_0, TAU_0 * 3, duration_steps)
    end_idx = min(start_idx + duration_steps, len(df))
    tau_0_values[start_idx:end_idx] = ramp[:end_idx - start_idx]
    tau_0_values[end_idx:] = TAU_0 * 3

    # Recalculer thermique avec tau_0 variable (approximation par segments)
    theta_oil, theta_hs = simulate_thermal_model(
        df_anom["K"].values, df_anom["theta_amb"].values, dt=10,
        tau_0_override=tau_0_values.mean()  # simplification pour demo
    )
    df_anom["theta_oil"] = theta_oil
    df_anom["theta_hs"] = theta_hs
    return df_anom

What does this code do?

This simulates a gradual clogging of the transformer's air-cooling system. As dust and debris accumulate on radiator fins, the thermal time constant tau_0 increases (the oil takes longer to cool down). Over 24 steps (4 hours), tau_0 triples from 150 to 450 minutes. This is the kind of slow degradation that is invisible to a simple temperature threshold until it is too late.

  • Run the cell, then in a new cell define the dissolved gas (H2) fault injection:
def inject_h2_defaut(df, start_idx, duration_steps=18, target_ppm=150):
    """Defaut naissant : H2 monte progressivement vers target_ppm."""
    df_anom = df.copy()
    end_idx = start_idx + duration_steps
    base = df_anom.loc[start_idx, "H2_ppm"]
    ramp = np.linspace(base, target_ppm, duration_steps)
    df_anom.loc[start_idx:end_idx-1, "H2_ppm"] = ramp
    return df_anom

What does this code do?

This simulates an incipient internal fault (partial discharge or hot spot in the windings). Dissolved hydrogen in the oil rises progressively from its baseline (~7 ppm) to 150 ppm over 3 hours. This pattern precedes catastrophic insulation failure by days or weeks in real transformers.

  • Run the cell. Verify the H2 injection works by running this test:
# Test rapide de l'injection H2
df_test_h2 = inject_h2_defaut(df_normal, start_idx=2000, duration_steps=18)
print(f"H2 avant injection (idx 1999): {df_test_h2.loc[1999, 'H2_ppm']:.1f} ppm")
print(f"H2 a la fin de l'injection (idx 2017): {df_test_h2.loc[2017, 'H2_ppm']:.1f} ppm")

You should see H2 at around 7 ppm before the injection point and 150 ppm at the end. The ramp simulates a 3-hour progressive gas buildup.

Seeing a KeyError on the loc index?

Make sure your df_normal uses default integer indexing (0 to 4319). If you reset the index with set_index("timestamp") earlier, the loc accessor expects datetime values instead of integers. Re-run build_normal_dataset without modifying the index.

Help me fix the index error in my injection function

Observe the limits of static thresholds

Your 4 injection functions are ready. Now comes the critical question: can the DGPT2's fixed thresholds catch all of these faults in time?

  • Generate all 4 anomaly scenarios at the same injection point by running this cell:
# Generer les 4 scenarios d'anomalie (injection a l'index 2000)
start = 2000
scenarios = {
    "Surcharge progressive": inject_surcharge(df_normal, start, duration_steps=12),
    "Desequilibre phases": inject_desequilibre(df_normal, start, duration_steps=18),
    "Colmatage aero-refrigerant": inject_colmatage(df_normal, start, duration_steps=24),
    "Defaut H2 naissant": inject_h2_defaut(df_normal, start, duration_steps=18),
}
print(f"4 scenarios generes, injection a l'index {start} (t = {df_normal.loc[start, 'timestamp']})")

Why inject at index 2000?

Index 2000 corresponds to roughly 14 days into the month of data. This ensures the thermal model has reached a stable periodic regime before the anomaly begins, so the fault signature is not confounded with initial transients.

Before you run the DGPT2 analysis, take a moment to predict: do you think the fixed thresholds (95 degrees C alarm, 105 degrees C trip) will detect all 4 scenarios equally well? Which ones might escape detection entirely?

  • Apply the DGPT2 detection to every scenario and measure detection delays by running this cell:
# Appliquer DGPT2 a chaque scenario et mesurer les delais
print("--- DETECTION PAR SEUILS DGPT2 ---")
print(f"{'Scenario':<30} {'Alarme 95C (min)':<20} {'Decl. 105C (min)':<20} {'Alarme H2'}")
print("-" * 85)

for name, df_anom in scenarios.items():
    dgpt2 = detect_dgpt2(df_anom)

    # Calculer delai depuis debut anomalie (start=2000, dt=10min)
    alarme_idx = dgpt2["alarme_temp"].idxmax() if dgpt2["alarme_temp"].any() else None
    decl_idx = dgpt2["declenchement_temp"].idxmax() if dgpt2["declenchement_temp"].any() else None
    h2_idx = dgpt2["alarme_h2"].idxmax() if dgpt2["alarme_h2"].any() else None

    alarme_min = (alarme_idx - start) * 10 if alarme_idx and alarme_idx >= start else "Jamais"
    decl_min = (decl_idx - start) * 10 if decl_idx and decl_idx >= start else "Jamais"
    h2_min = (h2_idx - start) * 10 if h2_idx and h2_idx >= start else "Jamais"

    print(f"{name:<30} {str(alarme_min):<20} {str(decl_min):<20} {str(h2_min)}")

Look at the output table. The surcharge scenario eventually triggers the temperature alarm after around 90 minutes. But look at the colmatage and desequilibre rows.

The DGPT2 blind spot is now visible

For progressive scenarios like air-cooler clogging and phase imbalance, the DGPT2 either triggers extremely late or prints "Jamais" (never). The temperature never reaches 95 degrees C because the degradation is subtle enough to stay below the static line. The H2 fault does eventually trigger its own alarm, but only once the concentration has already reached dangerous levels.

This is exactly the gap that predictive maintenance fills. A static threshold has no memory and no sense of trend. It cannot see that the temperature is rising faster than it should given the current load. In the next steps, you will train ML models that detect this drift 30-45 minutes earlier.

✔️ Awesome, I've got everything!

Great. Double check you have saved all cells in your notebook.

ⓧ I'd like to double check the full code

Scroll through your notebook and confirm these cells are present in order:

  • DGPT2 constants (SEUIL_ALARME_TEMP, SEUIL_DECLENCHEMENT_TEMP, SEUIL_H2_ALARME) and the detect_dgpt2 function.
  • Verification cell printing 0 alarms on df_normal.
  • inject_surcharge function (K ramp + thermal recalculation).
  • inject_desequilibre function (I1 diverges 30% over 18 steps).
  • Test cell printing theta_oil max for surcharge scenario.
  • inject_colmatage function (tau_0 triples over 24 steps).
  • inject_h2_defaut function (H2 ramp to 150 ppm).
  • Test cell printing H2 values before and after injection.
  • Scenario generation cell (start = 2000, 4 scenarios in a dictionary).
  • DGPT2 analysis loop printing delay table.

You have now quantified the exact gap in static protection. The DGPT2 is blind to progressive degradations that stay below its fixed line. In the next step, you will train an Isolation Forest that detects anomalies in the multivariate space and catches these drifts 30+ minutes earlier.

Detect Anomalies with Isolation Forest

Vous avez constate les limites du DGPT2 : sur les degradations progressives comme le colmatage ou le defaut H2, les seuils fixes ne declenchent pas ou declenchent trop tard.

L'Isolation Forest detecte les anomalies dans l'espace multivarie SANS avoir besoin d'exemples de pannes. Il apprend la "foret normale" et isole les points qui s'en ecartent. C'est exactement la realite operationnelle Enedis : beaucoup de donnees normales, tres peu d'historique de pannes labellise.

In this step, get ready to:
  • Entrainer un Isolation Forest uniquement sur les donnees de fonctionnement normal.
  • Scorer les 4 scenarios d'anomalies et comparer les delais de detection avec le DGPT2.
  • Calculer precision, rappel et F1 pour evaluer les forces et limites du modele.
Preparer les features et entrainer sur les donnees normales

L'Isolation Forest fonctionne en construisant des arbres de decision aleatoires. Les observations normales necessitent plus de coupures pour etre isolees (profondeur elevee).

Les anomalies, etant rares et differentes, sont isolees rapidement (faible profondeur). Le score d'anomalie est inversement proportionnel a la profondeur moyenne d'isolation.

  • Dans une nouvelle cellule de votre notebook, definissez la liste de features et la fonction d'entrainement en collant ce code :
# Features utilisees pour la detection d'anomalies
FEATURES = ["K", "theta_oil", "theta_hs", "H2_ppm", "desequilibre_pct", "I1", "I2"]

def train_isolation_forest(df_normal, features=FEATURES, contamination=0.01, n_estimators=200, seed=42):
    """Entraine un Isolation Forest sur les donnees normales."""
    scaler = StandardScaler()
    X_train = scaler.fit_transform(df_normal[features].values)

    iso_forest = IsolationForest(
        n_estimators=n_estimators,
        contamination=contamination,
        random_state=seed
    )
    iso_forest.fit(X_train)

    return iso_forest, scaler

# Entrainer sur les donnees normales
iso_forest, scaler_if = train_isolation_forest(df_normal)
print(f"Isolation Forest entraine sur {len(df_normal)} echantillons")
print(f"Features: {FEATURES}")

What does this code do?

FEATURES definit les 7 variables d'entree du modele : taux de charge, temperatures, H2, desequilibre et courants.

  • StandardScaler normalise chaque feature (moyenne 0, ecart-type 1) pour que toutes pesent egalement dans l'arbre.
  • contamination=0.01 fixe le percentile du score en dessous duquel un point est considere anomal (1% de faux positifs attendus).
  • n_estimators=200 utilise 200 arbres pour stabiliser le score d'isolation.
  • Executez la cellule.
  • Vous devriez voir Isolation Forest entraine sur 4320 echantillons confirmer que le modele a appris sur le mois complet de donnees normales.

Erreur sur StandardScaler ou IsolationForest?

  • Verifiez que vos imports du Step 1 incluent bien from sklearn.preprocessing import StandardScaler et from sklearn.ensemble import IsolationForest.
  • Confirmez que df_normal existe en executant df_normal.shape dans une cellule.

Pourquoi mon IsolationForest ne s'entraine pas?

  • Dans la cellule suivante, definissez la fonction de scoring et calculez le seuil d'alerte :
def score_isolation_forest(iso_forest, scaler, df, features=FEATURES):
    """Score les donnees avec l'Isolation Forest entraine."""
    X = scaler.transform(df[features].values)
    scores = iso_forest.score_samples(X)
    # Inverser : plus le score est negatif, plus c'est anormal
    anomaly_score = -scores
    return anomaly_score

# Calculer le seuil sur les donnees normales (percentile 99)
scores_normal = score_isolation_forest(iso_forest, scaler_if, df_normal)
if_threshold = np.percentile(scores_normal, 99)
print(f"Seuil IF (percentile 99): {if_threshold:.4f}")

Pourquoi inverser le score?

score_samples() de scikit-learn retourne des valeurs negatives pour les anomalies (plus negatif = plus anormal). On inverse avec -scores pour avoir un score d'anomalie positif croissant, plus intuitif pour la visualisation.

Le seuil (percentile 99) signifie que seulement 1% des points normaux depassent cette valeur. Tout point au-dessus est signale comme anomalie.

  • Executez la cellule.
  • Vous devriez voir le seuil IF affiche (une valeur autour de 0.10-0.15).
Analyser les scores et definir les seuils d'alerte

Le modele est entraine. Scorons maintenant les 4 scenarios d'anomalies et comparons la rapidite de detection de l'Isolation Forest face au DGPT2.

  • Construisez le tableau comparatif DGPT2 vs Isolation Forest en collant ce code dans une nouvelle cellule :
# Comparaison des delais de detection : DGPT2 vs Isolation Forest
start = 2000

print(f"{'Scenario':<30} {'DGPT2 (min)':<15} {'IF (min)':<12} {'Gain ML'}")
print("-" * 65)

for name, df_anom in scenarios.items():
    # Score Isolation Forest
    if_scores = score_isolation_forest(iso_forest, scaler_if, df_anom)
    if_anomalies = if_scores > if_threshold
    first_if = np.where(if_anomalies[start:])[0]
    first_if_min = first_if[0] * 10 if len(first_if) > 0 else None

    # DGPT2
    dgpt2_results = detect_dgpt2(df_anom)
    first_dgpt2 = dgpt2_results["alarme_temp"].idxmax() if dgpt2_results["alarme_temp"].any() else None
    dgpt2_min = (first_dgpt2 - start) * 10 if first_dgpt2 is not None and first_dgpt2 >= start else "Jamais"

    # Gain en minutes
    if isinstance(dgpt2_min, (int, float)) and first_if_min is not None:
        gain = f"{dgpt2_min - first_if_min} min"
    else:
        gain = "inf (DGPT2 ne detecte pas)"

    print(f"{name:<30} {str(dgpt2_min):<15} {str(first_if_min):<12} {gain}")

Comment lire ce tableau?

  • DGPT2 (min) : nombre de minutes entre le debut de l'anomalie et le franchissement du seuil de 95 degres C.
  • IF (min) : nombre de minutes entre le debut de l'anomalie et le premier point ou le score IF depasse le seuil.
  • Gain ML : avance de l'Isolation Forest sur le DGPT2 en minutes. "inf" signifie que le DGPT2 ne detecte jamais.
  • Executez la cellule et observez les delais de detection.
  • Pour la surcharge, l'IF devrait detecter 30+ minutes avant le DGPT2.
  • Pour le desequilibre et le defaut H2, le DGPT2 affiche "Jamais" car la temperature ne monte pas assez.
  • Visualisez la timeline de detection sur le scenario colmatage en collant ce code dans la cellule suivante :
# Visualisation : score IF vs temperature pour le colmatage
fig, axes = plt.subplots(2, 1, figsize=(12, 6), sharex=True)

df_plot = scenarios["Colmatage aero-refrigerant"]
if_scores_colm = score_isolation_forest(iso_forest, scaler_if, df_plot)

# Score IF dans le temps
axes[0].plot(if_scores_colm, linewidth=0.5)
axes[0].axhline(y=if_threshold, color='r', linestyle='--', label='Seuil IF (p99)')
axes[0].axvline(x=start, color='g', linestyle=':', label='Debut anomalie')
axes[0].set_ylabel("Score IF")
axes[0].set_title("Detection Isolation Forest - Colmatage aero-refrigerant")
axes[0].legend()

# Temperature huile avec seuils DGPT2
axes[1].plot(df_plot["theta_oil"].values, linewidth=0.5)
axes[1].axhline(y=95, color='orange', linestyle='--', label='Alarme DGPT2 (95C)')
axes[1].axhline(y=105, color='red', linestyle='--', label='Declenchement (105C)')
axes[1].axvline(x=start, color='g', linestyle=':', label='Debut anomalie')
axes[1].set_ylabel("Temperature huile (C)")
axes[1].set_xlabel("Pas de temps (10 min)")
axes[1].legend()

plt.tight_layout()
plt.show()
  • Executez la cellule.
  • Sur le graphe du haut, le score IF devrait depasser la ligne rouge pointillee (seuil) bien avant que la temperature (graphe du bas) n'atteigne la ligne orange du DGPT2.

Le graphe est vide ou le score reste plat?

  • Verifiez que votre dictionnaire scenarios contient bien la cle "Colmatage aero-refrigerant" (attention aux accents).
  • Confirmez que start = 2000 correspond bien au point d'injection dans vos donnees.

Aide-moi a debugger mon graphe Isolation Forest

Calculer precision, rappel et limites de l'Isolation Forest

Le tableau des delais montre l'avance temporelle de l'IF. Calculons maintenant la precision (combien d'alarmes sont de vraies anomalies) et le rappel (combien d'anomalies reelles sont detectees) pour quantifier la fiabilite du modele.

  • Calculez les metriques par scenario en collant ce code dans une nouvelle cellule :
# Precision / Rappel / F1 par scenario d'anomalie
from sklearn.metrics import precision_score, recall_score, f1_score

# Duree de chaque anomalie injectee (en pas de 10 min)
durations = {
    "Surcharge progressive": 12,
    "Desequilibre phases": 18,
    "Colmatage aero-refrigerant": 24,
    "Defaut H2 naissant": 18
}

print(f"{'Scenario':<30} {'Precision':<12} {'Rappel':<12} {'F1':<12}")
print("-" * 66)

for name, df_anom in scenarios.items():
    if_scores = score_isolation_forest(iso_forest, scaler_if, df_anom)
    y_pred = (if_scores > if_threshold).astype(int)

    # Verite terrain : 1 dans la zone d'injection, 0 ailleurs
    y_true = np.zeros(len(df_anom))
    y_true[start:start + durations[name]] = 1

    prec = precision_score(y_true, y_pred, zero_division=0)
    rec = recall_score(y_true, y_pred, zero_division=0)
    f1 = f1_score(y_true, y_pred, zero_division=0)

    print(f"{name:<30} {prec:<12.3f} {rec:<12.3f} {f1:<12.3f}")

What does this code do?

  • y_true contient 1 dans la zone d'injection de l'anomalie et 0 partout ailleurs (c'est la verite terrain).
  • y_pred contient 1 quand le score IF depasse le seuil (prediction du modele).
  • La precision mesure "quand l'IF alerte, a-t-il raison?". Le rappel mesure "quand il y a une anomalie reelle, l'IF la voit-il?".
  • Executez la cellule et observez les resultats.

Selon vous, sur quel scenario l'Isolation Forest a-t-il le meilleur F1? Et sur lequel le rappel est-il le plus faible?

L'IF devrait avoir un bon F1 sur la surcharge (anomalie franche dans l'espace des features). En revanche, le rappel sur le colmatage est souvent faible : la derive thermique est lente et subtile, et l'IF raisonne sur des coupures univariees qui manquent les correlations multi-variables.

Pourquoi l'IF manque le colmatage?

Le colmatage degrade le refroidissement progressivement. La temperature monte lentement, en restant dans la fourchette "plausible" feature par feature.

Ce que le DGPT2 et l'IF ratent tous les deux, c'est la correlation brisee : la temperature monte SANS que la charge augmente. C'est exactement ce qu'un Autoencoder peut capturer, car il apprend les correlations entre variables.

Revenons sur le scenario de colmatage pour la verification finale de ce step. Avant d'executer, a votre avis : l'IF detecte-t-il le colmatage au moins 30 minutes avant que la temperature n'atteigne 95 degres C?

  • Verifiez le delai precis en collant ce code :
# Verification : gain IF vs DGPT2 sur le colmatage
df_colm = scenarios["Colmatage aero-refrigerant"]
if_scores_v = score_isolation_forest(iso_forest, scaler_if, df_colm)
first_if_idx = np.where(if_scores_v[start:] > if_threshold)[0]

# Moment ou la temperature depasse 95 degC
first_95 = np.where(df_colm["theta_oil"].values[start:] > 95)[0]

if len(first_if_idx) > 0:
    print(f"IF detecte a t={first_if_idx[0] * 10} min apres debut anomalie")
if len(first_95) > 0:
    print(f"DGPT2 alarme a t={first_95[0] * 10} min apres debut anomalie")
    print(f"Gain IF: {(first_95[0] - first_if_idx[0]) * 10} minutes d'avance")
else:
    print("DGPT2 ne declenche pas sur ce scenario!")
    print(f"IF detecte quand meme a t={first_if_idx[0] * 10} min")

L'Isolation Forest devrait signaler le colmatage 30+ minutes avant le seuil de 95 degres C du DGPT2. C'est deja un gain significatif pour la maintenance predictive.

La limite que nous allons depasser

L'IF detecte bien les anomalies franches (surcharge, desequilibre), mais son rappel reste modeste sur les derives lentes multivarieees. Quand la charge est faible ET que le colmatage est progressif, l'IF peut passer a cote.

C'est la motivation directe pour l'Autoencoder du prochain step : il capture les correlations entre variables et explose en erreur de reconstruction quand ces correlations se brisent.

✔️ Awesome, I've got everything!

Parfait. Verifiez que vous avez bien sauvegarde votre notebook.

ⓧ I'd like to double check the full code

Remontez dans votre notebook et confirmez que ces cellules sont presentes, dans l'ordre :

  • Cellule 1 : definition de FEATURES + fonction train_isolation_forest + appel d'entrainement.
  • Cellule 2 : fonction score_isolation_forest + calcul de if_threshold.
  • Cellule 3 : boucle de comparaison DGPT2 vs IF pour les 4 scenarios.
  • Cellule 4 : visualisation timeline (2 subplots).
  • Cellule 5 : precision / rappel / F1 par scenario.
  • Cellule 6 : verification du gain IF vs DGPT2 sur le colmatage.

L'Isolation Forest vous donne deja une avance precieuse sur les seuils fixes. Dans le prochain step, vous allez construire un Autoencoder qui apprend les correlations entre variables et surpasse l'IF de 10-15 minutes sur les derives lentes.

Detect Anomalies with an Autoencoder

Dans l'etape precedente, l'Isolation Forest a detecte les anomalies franches (surcharge) mais a montre ses limites sur les derives lentes multivarieees comme le colmatage aero-refrigerant.

L'Autoencoder resout ce probleme en apprenant a compresser puis reconstruire les donnees normales. Quand une anomalie jamais vue arrive, la reconstruction echoue et l'erreur explose. Contrairement a l'Isolation Forest qui raisonne sur des coupures univariees, l'Autoencoder capture les correlations multi-variables.

Un exemple concret : si la temperature monte SANS que la charge augmente, cela signifie un refroidissement defaillant. L'Autoencoder detecte cette correlation brisee la ou l'IF ne voit que des valeurs individuelles encore dans la norme.

In this step, get ready to:
  • Construire et entrainer un Autoencoder sur les donnees normales.
  • Comparer les performances de detection DGPT2 vs Isolation Forest vs Autoencoder.
  • Construire un Health Index composite pour synthetiser l'etat du transformateur.
Construire et entrainer l'Autoencoder

Un Autoencoder est un reseau de neurones en forme de sablier. L'encodeur compresse les 7 features en un espace latent de dimension 3. Le decodeur tente de reconstruire les 7 features originales a partir de cette representation comprimee.

Entraine uniquement sur des donnees normales, le reseau apprend la structure interne des correlations normales. Une anomalie brise ces correlations et genere une erreur de reconstruction elevee.

  • Dans votre notebook, creez une nouvelle cellule et definissez la fonction build_autoencoder en copiant ce code :
def build_autoencoder(input_dim=7, encoding_dim=3):
    """Construit un Autoencoder pour detection d'anomalies."""
    input_layer = Input(shape=(input_dim,))
    encoded = Dense(16, activation="relu")(input_layer)
    encoded = Dense(8, activation="relu")(encoded)
    encoded = Dense(encoding_dim, activation="relu")(encoded)
    decoded = Dense(8, activation="relu")(encoded)
    decoded = Dense(16, activation="relu")(decoded)
    decoded = Dense(input_dim, activation="linear")(decoded)

    autoencoder = KerasModel(input_layer, decoded)
    autoencoder.compile(optimizer=Adam(learning_rate=0.001), loss="mse")

    return autoencoder

What does this architecture do?

  • L'architecture forme un sablier : Input(7) -> Dense(16) -> Dense(8) -> Dense(3) -> Dense(8) -> Dense(16) -> Dense(7). Le goulot d'etranglement a 3 neurones force le reseau a capturer l'essentiel.
  • La derniere couche utilise une activation linear (pas de relu) pour pouvoir reconstruire des valeurs negatives (features normalisees).
  • L'optimiseur Adam avec lr=0.001 et la loss MSE mesurent l'ecart entre l'entree et la sortie reconstruite.
  • Dans la cellule suivante, definissez la fonction d'entrainement train_autoencoder qui entraine sur les donnees normales et calcule le seuil d'anomalie :
def train_autoencoder(df_normal, features=FEATURES, epochs=50, batch_size=32, seed=42):
    """Entraine l'Autoencoder sur les donnees normales."""
    np.random.seed(seed)

    scaler = StandardScaler()
    X_train = scaler.fit_transform(df_normal[features].values)

    autoencoder = build_autoencoder(input_dim=len(features))
    autoencoder.fit(
        X_train, X_train,
        epochs=epochs,
        batch_size=batch_size,
        validation_split=0.1,
        verbose=0
    )

    # Calcul du seuil (percentile 99 sur donnees normales)
    X_pred = autoencoder.predict(X_train, verbose=0)
    mse = np.mean((X_train - X_pred)**2, axis=1)
    threshold = np.percentile(mse, 99)

    return autoencoder, scaler, threshold

Why train on X_train twice?

L'Autoencoder apprend a reconstruire son entree. L'input ET le target sont les memes donnees (X_train, X_train). La loss MSE penalise toute difference entre la sortie et l'entree originale.

Le seuil au percentile 99 signifie que 99% des points normaux ont une erreur inferieure a ce seuil. Tout point au-dessus est signale comme anomalie potentielle.

  • Lancez l'entrainement en executant cette cellule :
# Entrainer l'Autoencoder sur les donnees normales
print("Entrainement Autoencoder...")
autoencoder, scaler_ae, ae_threshold = train_autoencoder(df_normal)
print(f"Seuil AE (p99): {ae_threshold:.6f}")
  • Executez la cellule. L'entrainement prend environ 10-20 secondes. Vous devriez voir le seuil s'afficher (une valeur decimale autour de 0.01-0.05).

Erreur d'import TensorFlow?

Si vous voyez une erreur mentionnant KerasModel ou Input non defini, verifiez que vos imports du debut du notebook incluent bien les lignes TensorFlow/Keras.

help me debug my TensorFlow Keras import error in Jupyter

Calculer l'erreur de reconstruction et scorer les scenarios

L'erreur de reconstruction est la MSE (Mean Squared Error) entre l'input et l'output, calculee point par point. Tout point dont l'erreur depasse le seuil calcule a l'etape precedente est signale comme anomalie.

  • Definissez la fonction score_autoencoder dans une nouvelle cellule :
def score_autoencoder(autoencoder, scaler, df, features=FEATURES):
    """Calcule l'erreur de reconstruction (score d'anomalie)."""
    X = scaler.transform(df[features].values)
    X_pred = autoencoder.predict(X, verbose=0)
    mse = np.mean((X - X_pred)**2, axis=1)
    return mse

What does this function return?

Pour chaque instant t, la fonction retourne un scalaire : l'erreur de reconstruction MSE. Plus cette valeur est elevee, plus le point s'ecarte du comportement normal appris par le reseau.

  • Tracez l'erreur de reconstruction pour les 4 scenarios d'anomalie en executant cette cellule :
# Tracer l'erreur de reconstruction pour chaque scenario
fig, axes = plt.subplots(2, 2, figsize=(14, 8))
start = 2000

for ax, (name, df_anom) in zip(axes.flat, scenarios.items()):
    ae_scores = score_autoencoder(autoencoder, scaler_ae, df_anom)
    ax.plot(ae_scores, label="Erreur reconstruction", alpha=0.7)
    ax.axhline(y=ae_threshold, color="r", linestyle="--", label=f"Seuil p99")
    ax.axvline(x=start, color="g", linestyle=":", label="Debut anomalie")
    ax.set_title(name)
    ax.legend(fontsize=8)
    ax.set_xlabel("Temps (pas de 10 min)")
    ax.set_ylabel("MSE")

plt.tight_layout()
plt.show()
  • Executez la cellule. Vous devriez voir 4 graphiques ou l'erreur de reconstruction explose apres le debut de l'anomalie (ligne verte verticale), depassant le seuil (ligne rouge) bien avant le point d'injection + 90 minutes.

L'erreur ne depasse pas le seuil?

Verifiez que scaler_ae (et non scaler_if) est passe a score_autoencoder. Les deux scalers sont distincts car l'IF et l'AE normalisent independamment.

help me debug my autoencoder reconstruction error not detecting anomalies

Comparer les methodes et construire le Health Index

Il est temps de mettre les trois approches en face a face. Pour chaque scenario, vous allez mesurer le delai de detection et quantifier le gain du ML sur les seuils statiques.

  • Construisez le tableau comparatif final en executant cette cellule :
# Comparaison DGPT2 vs IF vs AE
print("--- COMPARAISON DES METHODES DE DETECTION ---")
print(f"{'Scenario':<30} {'DGPT2 (min)':<15} {'IF (min)':<12} {'AE (min)':<12} {'Gain ML'}")
print("-" * 80)

start = 2000
if_threshold = np.percentile(
    score_isolation_forest(iso_forest, scaler_if, df_normal), 99
)

for name, df_anom in scenarios.items():
    # DGPT2
    dgpt2_results = detect_dgpt2(df_anom)
    first_dgpt2 = dgpt2_results["alarme_temp"].idxmax() if dgpt2_results["alarme_temp"].any() else None
    dgpt2_min = (first_dgpt2 - start) * 10 if first_dgpt2 is not None and first_dgpt2 >= start else "Jamais"
    # Isolation Forest
    if_scores = score_isolation_forest(iso_forest, scaler_if, df_anom)
    if_anomalies = if_scores > if_threshold
    first_if = np.where(if_anomalies[start:])[0]
    first_if_min = first_if[0] * 10 if len(first_if) > 0 else None

    # Autoencoder
    ae_scores = score_autoencoder(autoencoder, scaler_ae, df_anom)
    ae_anomalies = ae_scores > ae_threshold
    first_ae = np.where(ae_anomalies[start:])[0]
    first_ae_min = first_ae[0] * 10 if len(first_ae) > 0 else None

    # Gain ML
    ml_best = min(x for x in [first_if_min, first_ae_min] if x is not None) if any(
        x is not None for x in [first_if_min, first_ae_min]) else "N/A"
    gain = f"{dgpt2_min - ml_best} min" if isinstance(dgpt2_min, (int, float)) and isinstance(
        ml_best, (int, float)) else "inf"

    print(f"{name:<30} {str(dgpt2_min):<15} {str(first_if_min):<12} {str(first_ae_min):<12} {gain}")

What does this comparison show?

Pour chaque scenario, le tableau affiche en minutes le delai entre le debut de l'anomalie et la premiere detection. Le "Gain ML" montre combien de minutes d'anticipation le meilleur modele ML offre par rapport au DGPT2.

"Jamais" signifie que le DGPT2 ne declenche pas du tout sur ce scenario. C'est exactement le cas ou le ML apporte le plus de valeur.

Avant d'executer cette cellule, a votre avis : sur quel scenario l'Autoencoder va-t-il le plus surpasser l'Isolation Forest?

  • Executez les deux cellules. Observez le tableau : sur le scenario de colmatage aero-refrigerant, l'erreur de reconstruction de l'Autoencoder depasse le seuil d'alerte 45+ minutes avant le DGPT2, surpassant l'Isolation Forest de 10-15 minutes.

L'Autoencoder excelle sur le colmatage car il detecte que la correlation temperature/charge est brisee. L'IF ne voyait que des valeurs individuelles encore dans la plage normale.

  • Maintenant, construisez le Health Index composite qui synthetise l'etat global du transformateur en copiant cette fonction :
def compute_health_index(if_score, ae_score, load_profile, ae_threshold,
                         w_if=0.3, w_ae=0.5, w_load=0.2):
    """
    Calcule le Health Index composite.
    HI = w_if * IF_norm + w_ae * AE_norm + w_load * K_norm
    Normalise entre 0 (sain) et 1 (critique).
    """
    # Normaliser IF score entre 0 et 1
    if_norm = (if_score - if_score.min()) / (if_score.max() - if_score.min() + 1e-6)
    # Normaliser AE score par rapport au seuil
    ae_norm = np.clip(ae_score / (ae_threshold * 3), 0, 1)
    # Normaliser charge (>1 = surcharge)
    load_norm = np.clip(load_profile / 1.2, 0, 1)

    hi = w_if * if_norm + w_ae * ae_norm + w_load * load_norm
    return np.clip(hi, 0, 1)

Why these weights?

Le poids w_ae=0.5 donne la priorite a l'Autoencoder car il capture les correlations multivarieees. Le poids w_if=0.3 renforce avec l'IF qui excelle sur les outliers francs. Le poids w_load=0.2 penalise la surcharge. Cette ponderation s'aligne sur la methodologie Enedis CIRED (scoring multi-critere pondere par la charge).

  • Tracez l'evolution du Health Index sur le scenario de colmatage en executant cette cellule :
# Health Index sur le scenario colmatage
df_test = scenarios["Colmatage aero-refrigerant"]
if_scores_test = score_isolation_forest(iso_forest, scaler_if, df_test)
ae_scores_test = score_autoencoder(autoencoder, scaler_ae, df_test)
hi = compute_health_index(if_scores_test, ae_scores_test, df_test["K"].values, ae_threshold)

print(f"HI moyen: {hi.mean():.3f}")
print(f"HI max: {hi.max():.3f}")
print(f"HI au point d'injection: {hi[start]:.3f}")
print(f"HI 1h apres injection: {hi[start + 6]:.3f}")
# Visualiser l'evolution du HI
fig, ax = plt.subplots(figsize=(12, 4))
ax.plot(hi, label="Health Index", color="purple", alpha=0.8)
ax.axhline(y=0.3, color="orange", linestyle="--", label="Seuil surveillance")
ax.axhline(y=0.6, color="red", linestyle="--", label="Seuil intervention")
ax.axvline(x=start, color="green", linestyle=":", label="Debut colmatage")
ax.set_xlabel("Temps (pas de 10 min)")
ax.set_ylabel("Health Index")
ax.set_title("Evolution du Health Index - Colmatage aero-refrigerant")
ax.legend()
ax.set_ylim(0, 1)
plt.tight_layout()
plt.show()
  • Executez les deux cellules. Vous devriez voir le HI rester bas (~0.1-0.2) en fonctionnement normal, puis monter en fleche apres le debut du colmatage, franchissant le seuil orange (0.3 = surveillance renforcee) puis le seuil rouge (0.6 = intervention prioritaire).

Le HI reste plat apres l'injection?

Verifiez que le dictionnaire scenarios contient bien la cle exacte "Colmatage aero-refrigerant" (avec accent). Verifiez aussi que start = 2000 correspond bien au point d'injection utilise lors de la creation des scenarios.

help me debug my Health Index computation that stays flat

✔️ Awesome, I've got everything!

Parfait. Verifiez que toutes vos cellules sont sauvegardees dans le notebook.

ⓧ I'd like to double check the full code

Remontez dans votre notebook et confirmez que ces blocs sont presents dans l'ordre :

  • La fonction build_autoencoder(input_dim=7, encoding_dim=3) avec l'architecture sablier 7->16->8->3->8->16->7.
  • La fonction train_autoencoder qui normalise, entraine (50 epochs, batch 32), et calcule le seuil p99.
  • La cellule d'entrainement produisant autoencoder, scaler_ae, et ae_threshold.
  • La fonction score_autoencoder retournant le MSE par point.
  • Les 4 graphiques d'erreur de reconstruction par scenario.
  • Le tableau comparatif DGPT2 vs IF vs AE.
  • La fonction compute_health_index avec poids w_if=0.3, w_ae=0.5, w_load=0.2.
  • Le graphique du HI pour le scenario colmatage avec les seuils 0.3 et 0.6.

Vous avez demontre qu'un pipeline ML non-supervise (Isolation Forest + Autoencoder) detecte une degradation progressive du transformateur 30 a 45 minutes avant les seuils statiques du DGPT2. C'est la preuve de concept de la maintenance predictive. Prochaine etape : generaliser cette approche a un parc de transformateurs dans la Secret Mission.

Secret mission

Tableau de Priorisation Maintenance d'un Parc de 5 Transformateurs

En exploitation Enedis, l'enjeu n'est pas de surveiller UN transformateur mais de prioriser les interventions sur un PARC entier. Construisez un tableau de ranking par Health Index pour 5 postes HTA/BT, exactement comme le ferait un ingenieur exploitation.

Clean Up Your Resources

Clean Up Your Resources

Decide whether to keep your resources running, pause them to come back later, or delete them entirely. This project runs entirely locally with no cloud resources, so there are no ongoing costs.

Resources you used:

  • Jupyter notebook anomaly_detection_hta_bt.ipynb.
  • Python file anomaly_detection_hta_bt_fleet.py (Secret Mission).
  • Generated image fleet_ranking.png (Secret Mission).
  • Python packages installed via pip (numpy, pandas, scikit-learn, tensorflow, matplotlib).

Keep everything running

No action needed. Choose this if you plan to extend your models or connect them to real GridWatch/Sentinel data.

  • Your notebook is self-contained and can be reused on real monitoring data.
  • The pip packages (numpy, pandas, scikit-learn, tensorflow, matplotlib) remain installed and available for other projects.
  • No processes remain running after you close Jupyter.

Pause - I'll come back to this later

Shut down Jupyter to free up memory. Your files stay on disk for when you return.

  • In the terminal where Jupyter is running, press Ctrl+C twice to stop the server.
  • Your notebook and Python files remain saved on your Desktop.
  • When you return, navigate to the project folder and run jupyter notebook to pick up where you left off.

Delete - I don't want to use this again

Remove all project files and optionally uninstall the Python packages.

  • Stop Jupyter by pressing Ctrl+C twice in the terminal.
  • Delete the project folder by running:

macOS

cd ~/Desktop
rm -rf anomaly_detection_hta_bt.ipynb anomaly_detection_hta_bt.py anomaly_detection_hta_bt_fleet.py fleet_ranking.png

Windows

cd ~\Desktop
Remove-Item anomaly_detection_hta_bt.ipynb, anomaly_detection_hta_bt.py, anomaly_detection_hta_bt_fleet.py, fleet_ranking.png
  • Optionally, uninstall the Python packages by running:
pip uninstall numpy pandas scikit-learn tensorflow matplotlib -y

Keep common packages

Packages like numpy, pandas, and matplotlib are used across many Python projects. Only uninstall them if you are certain no other project depends on them.

Nice Work!

Nice Work!

Vous venez de construire un systeme complet de surveillance predictive pour un transformateur de distribution HTA/BT 630 kVA. Votre pipeline detecte des anomalies 30 a 45 minutes avant les seuils statiques du DGPT2, validant la plus-value du Machine Learning pour la maintenance predictive du parc Enedis.

Vous avez appris a :

  • Implementer le modele thermique CEI 60076-7 en Python pour simuler un transformateur de distribution, et constater les limites de la detection par seuils fixes (DGPT2) sur les degradations progressives.
  • Entrainer un Isolation Forest et un Autoencoder en mode non-supervise (uniquement sur des donnees normales) et demontrer un gain de detection de 30 a 45 minutes sur les seuils statiques.
  • Construire un Health Index composite qui synthetise temperature, charge et scores d'anomalie en un indicateur unique de priorisation maintenance.
  • Secret Mission : Produire un tableau de priorisation maintenance pour un parc de 5 transformateurs avec ranking par Health Index et recommandations d'intervention, exactement comme le ferait l'outil de scoring deploye par Enedis dans le cadre de son programme CIRED de monitoring en ligne.

Pret a tester vos connaissances ?