Diesel consumption anomalies:
What your tank reveals about your fleet

import pandas as pd
import numpy as np
from datetime import datetime, timedelta

np.random.seed(505)

LKW_COUNT = 150
WOCHEN = 52

# --- Fahrzeug-Basisprofile ---
fahrzeuge = []
for i in range(LKW_COUNT):
    baujahr = np.random.choice([2018,2019,2020,2021,2022,2023],
                                p=[0.08,0.12,0.20,0.25,0.20,0.15])
    motortyp = np.random.choice(['Euro6d', 'Euro6e'], p=[0.6, 0.4])
    basis_verbrauch = np.random.normal(30.5, 2.0)  # l/100km
    if baujahr <= 2019: basis_verbrauch += 1.5    # Ältere verbrauchen mehr

    # Anomalien: 12% der Fahrzeuge haben ein verstecktes Problem
    hat_anomalie = np.random.random() < 0.12
    anomalie_typ = None
    anomalie_start_woche = None
    if hat_anomalie:
        anomalie_typ = np.random.choice([
            'injektor_defekt',       # +3-5 l/100km
            'turbo_undicht',          # +2-4 l/100km
            'luftfilter_verstopft',   # +1.5-3 l/100km
            'reifendruck_niedrig',    # +1-2 l/100km
            'tempomat_defekt',        # +2-3 l/100km (ungleichmäßig)
        ])
        anomalie_start_woche = np.random.randint(10, 40)

    fahrzeuge.append({
        'lkw_id': f'LKW-{i+1:03d}',
        'baujahr': baujahr,
        'motortyp': motortyp,
        'basis_l100': basis_verbrauch,
        'hat_anomalie': hat_anomalie,
        'anomalie_typ': anomalie_typ,
        'anomalie_start': anomalie_start_woche,
    })

# --- Wochenweise Verbrauchsdaten ---
records = []
for fzg in fahrzeuge:
    for woche in range(1, WOCHEN + 1):
        km = np.random.normal(2200, 400)
        beladung = np.random.normal(18, 4)
        autobahn_anteil = np.random.uniform(0.3, 0.85)
        hoehenmeter = np.random.exponential(800)
        temperatur = 10 + 12 * np.sin((woche - 12) / 52 * 2 * np.pi)

        # Verbrauchsmodell: physikbasiert
        verbrauch = fzg['basis_l100']
        verbrauch += beladung * 0.12           # Beladungseffekt
        verbrauch -= autobahn_anteil * 2.5     # Autobahn effizienter
        verbrauch += hoehenmeter / 1000 * 1.8  # Höhenmeter-Effekt
        verbrauch += (20 - temperatur) * 0.04   # Kälte-Effekt
        verbrauch *= np.random.lognormal(0, 0.04)

        # Anomalie einfügen
        anomalie_aktiv = False
        if fzg['hat_anomalie'] and woche >= fzg['anomalie_start']:
            anomalie_aktiv = True
            wochen_seit_start = woche - fzg['anomalie_start']
            if fzg['anomalie_typ'] == 'injektor_defekt':
                verbrauch += 3.0 + wochen_seit_start * 0.15
            elif fzg['anomalie_typ'] == 'turbo_undicht':
                verbrauch += 2.5 + wochen_seit_start * 0.08
            elif fzg['anomalie_typ'] == 'luftfilter_verstopft':
                verbrauch += 1.5 + wochen_seit_start * 0.12
            elif fzg['anomalie_typ'] == 'reifendruck_niedrig':
                verbrauch += np.random.uniform(1.0, 2.0)
            elif fzg['anomalie_typ'] == 'tempomat_defekt':
                verbrauch += np.random.uniform(1.5, 4.0)

        liter_gesamt = verbrauch / 100 * km
        records.append({
            'kw': woche,
            'lkw_id': fzg['lkw_id'],
            'baujahr': fzg['baujahr'],
            'km': round(km),
            'beladung_t': round(beladung, 1),
            'autobahn_anteil': round(autobahn_anteil, 2),
            'hoehenmeter': round(hoehenmeter),
            'temperatur_c': round(temperatur, 1),
            'liter_gesamt': round(liter_gesamt, 1),
            'l_100km': round(verbrauch, 2),
            'anomalie_aktiv': anomalie_aktiv,
        })

df = pd.DataFrame(records)
print(f"Dataset: {len(df):,} Wochen-Records")
print(f"Fahrzeuge: {df['lkw_id'].nunique()}")
print(f"Ø Verbrauch: {df['l_100km'].mean():.1f} l/100km")
print(f"Anomalie-Wochen: {df['anomalie_aktiv'].sum():,}")

import xgboost as xgb
from sklearn.metrics import mean_absolute_error, r2_score

# Features die den "normalen" Verbrauch erklären
features = ['km', 'beladung_t', 'autobahn_anteil',
            'hoehenmeter', 'temperatur_c', 'baujahr', 'kw']

# Training nur auf normalen Wochen
normal = df[~df['anomalie_aktiv']]
X = normal[features]
y = normal['l_100km']

# Train/Test split: Zeitlich
train_mask = normal['kw'] <= 40
X_train, y_train = X[train_mask], y[train_mask]
X_test, y_test = X[~train_mask], y[~train_mask]

model = xgb.XGBRegressor(
    n_estimators=200, max_depth=5, learning_rate=0.08,
    subsample=0.8, colsample_bytree=0.8,
    random_state=42
)
model.fit(X_train, y_train)

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

print(f"Erwartungswert-Modell:")
print(f"  MAE:  {mae:.2f} l/100km")
print(f"  R²:   {r2:.3f}")

# Residuen für ALLE Daten berechnen (inkl. Anomalien)
df['erwartet_l100'] = model.predict(df[features])
df['residuum'] = df['l_100km'] - df['erwartet_l100']
df['abweichung_%'] = df['residuum'] / df['erwartet_l100'] * 100

from sklearn.ensemble import IsolationForest

# Features für Anomalie-Erkennung
anom_features = [
    'residuum',           # Abweichung vom Erwartungswert
    'abweichung_%',       # Prozentuale Abweichung
    'l_100km',            # Absoluter Verbrauch
    'beladung_t',         # Kontext
    'autobahn_anteil',    # Kontext
    'hoehenmeter',        # Kontext
]

# Rolling-Features: Trend über 4 Wochen
df['residuum_rolling_4w'] = df.groupby('lkw_id')['residuum'].transform(
    lambda x: x.rolling(4, min_periods=2).mean()
)
df['residuum_trend'] = df.groupby('lkw_id')['residuum_rolling_4w'].transform(
    lambda x: x.diff()
)
anom_features += ['residuum_rolling_4w', 'residuum_trend']

X_anom = df[anom_features].dropna()

# Isolation Forest: contamination = erwarteter Anomalie-Anteil
iso = IsolationForest(
    n_estimators=300,
    contamination=0.08,
    max_samples=0.8,
    random_state=42
)
df.loc[X_anom.index, 'anomalie_score'] = iso.fit_predict(X_anom)
df['ist_anomalie'] = (df['anomalie_score'] == -1)

# Evaluation: Vergleich mit echten Anomalien
tp = ((df['ist_anomalie']) & (df['anomalie_aktiv'])).sum()
fp = ((df['ist_anomalie']) & (~df['anomalie_aktiv'])).sum()
fn = ((~df['ist_anomalie']) & (df['anomalie_aktiv'])).sum()
precision = tp / (tp + fp)
recall = tp / (tp + fn)
print(f"Precision: {precision:.3f}")
print(f"Recall:    {recall:.3f}")
print(f"Anomalien erkannt: {df['ist_anomalie'].sum()}")

# === EINSPARUNG DURCH ANOMALIE-ERKENNUNG ===

# 1. Diesel-Mehrverbrauch vermeiden (frühere Reparatur)
# Durchschnittlich 8 Wochen frühere Erkennung als Werkstatt
lkw_mit_anomalie = 18
avg_mehrverbrauch_l_woche = 55  # l/Woche Mehrverbrauch
wochen_frueher = 8
diesel_preis = 1.42
ersparnis_diesel = lkw_mit_anomalie * avg_mehrverbrauch_l_woche * wochen_frueher * diesel_preis

# 2. Folgeschäden vermeiden
folgeschaden_rate = 0.35   # 35% der Defekte führen zu Folgeschaden
avg_folgeschaden = 6500    # € pro Folgeschaden
vermiedene_schaeden = int(lkw_mit_anomalie * folgeschaden_rate)
ersparnis_folge = vermiedene_schaeden * avg_folgeschaden

# 3. Pannenvermeidung
pannen_rate = 0.20        # 20% der unerkannten Defekte → Panne
kosten_panne = 2800       # Pannendienst + Standzeit + Umplanung
vermiedene_pannen = int(lkw_mit_anomalie * pannen_rate)
ersparnis_pannen = vermiedene_pannen * kosten_panne

# 4. Eco-Driving-Effekte (Fahrverhalten-Anomalien)
fahrverhalten_lkw = 5     # Fahrzeuge mit Fahrverhalten-Anomalie
ersparnis_eco = fahrverhalten_lkw * 2400  # €/Jahr pro Fahrer nach Schulung

total = ersparnis_diesel + ersparnis_folge + ersparnis_pannen + ersparnis_eco

print(f"Diesel-Einsparung:      €{ersparnis_diesel:>10,.0f}")
print(f"Folgeschäden vermieden: €{ersparnis_folge:>10,.0f}")
print(f"Pannen vermieden:       €{ersparnis_pannen:>10,.0f}")
print(f"Eco-Driving:            €{ersparnis_eco:>10,.0f}")
print(f"{'─'*40}")
print(f"GESAMT:                 €{total:>10,.0f}")