Automotive Ml

Anomaly Detection

Anomaly Detection for Automotive Systems

Detect unusual vehicle behavior across battery systems, sensors, and drivetrain components using unsupervised ML techniques.

Use Cases

1. Battery Anomalies: Cell voltage drift, thermal runaway precursors, SOC inconsistencies
2. Sensor Failures: LiDAR/radar malfunction, camera degradation, IMU drift
3. Drivetrain Issues: Motor vibration anomalies, inverter faults, cooling system failures
4. Charging Anomalies: Abnormal charging curves, connector issues, grid irregularities

Algorithm Selection

Isolation Forest

Best for: High-dimensional sensor data with mixed feature types

Pros:

• Handles non-Gaussian distributions
• Efficient for large datasets
• No assumptions about normal behavior shape
• Low memory footprint

Cons:

• Sensitive to feature scaling
• May struggle with local anomalies

Use cases: Real-time battery monitoring, sensor fault detection

Autoencoder (Deep Learning)

Best for: Complex time-series patterns, image-based anomalies

Pros:

• Learns compressed representation
• Excellent for time-series sequences
• Handles multi-modal data
• Can detect subtle pattern deviations

Cons:

• Requires significant training data
• Computationally expensive
• Black-box interpretation

Use cases: Camera degradation, LiDAR point cloud anomalies, battery degradation patterns

Local Outlier Factor (LOF)

Best for: Local density-based anomalies

Pros:

• Detects local outliers in varying density regions
• No global threshold needed
• Good for spatial data

Cons:

• Computationally intensive for large datasets
• Requires careful k-neighbor selection

Use cases: Geographic anomalies (GPS data), fleet-wide comparison

One-Class SVM

Best for: Small, well-defined normal behavior regions

Pros:

• Kernel trick for non-linear boundaries
• Robust to outliers in training set
• Theoretical foundation

Cons:

• Difficult hyperparameter tuning
• Slow on large datasets
• Memory intensive

Use cases: Safety-critical systems with narrow normal operating ranges

Feature Engineering

Battery Systems

python
features = [
    # Static features
    'cell_voltage_mean', 'cell_voltage_std', 'cell_voltage_range',
    'temperature_mean', 'temperature_std', 'temperature_max',
    'soc_value', 'soh_estimate', 'current_value',

    # Derived features
    'voltage_imbalance_ratio',  # max_voltage / min_voltage
    'thermal_gradient',          # max_temp - min_temp
    'charge_acceptance',         # dSOC / dTime at constant current

    # Time-series features (sliding window)
    'voltage_trend_5min',        # Linear regression slope
    'temperature_volatility',    # Rolling std deviation
    'current_spike_count',       # Threshold crossings

    # Cross-domain features
    'power_efficiency',          # Output power / Input power
    'thermal_power_ratio'        # Temperature rise / Power delivered
]

Sensor Systems

python
features = [
    # LiDAR
    'point_cloud_density', 'max_range', 'noise_level',
    'ground_plane_fit_error', 'object_count',

    # Camera
    'brightness_mean', 'contrast_std', 'edge_density',
    'motion_blur_score', 'lens_distortion_coefficient',

    # IMU
    'accel_magnitude_mean', 'gyro_drift_rate',
    'vibration_frequency_peak', 'calibration_offset',

    # GPS
    'hdop', 'satellite_count', 'position_jump_distance',
    'velocity_consistency_score'
]

Implementation Example: Isolation Forest for Battery

python
import numpy as np
import pandas as pd
from sklearn.ensemble import IsolationForest
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
import joblib
from typing import Dict, List, Tuple
import logging

logger = logging.getLogger(__name__)

class BatteryAnomalyDetector:
    """
    Isolation Forest-based anomaly detection for battery systems.

    Detects:
    - Cell voltage imbalances
    - Thermal anomalies
    - Charge/discharge irregularities
    - SOH degradation patterns
    """

    def __init__(
        self,
        contamination: float = 0.01,
        n_estimators: int = 100,
        max_samples: int = 256,
        random_state: int = 42
    ):
        """
        Args:
            contamination: Expected proportion of outliers (0.01 = 1%)
            n_estimators: Number of isolation trees
            max_samples: Samples per tree (higher = more global detection)
            random_state: Reproducibility seed
        """
        self.pipeline = Pipeline([
            ('scaler', StandardScaler()),
            ('detector', IsolationForest(
                contamination=contamination,
                n_estimators=n_estimators,
                max_samples=max_samples,
                random_state=random_state,
                n_jobs=-1
            ))
        ])
        self.feature_names = None
        self.trained = False

    def engineer_features(self, df: pd.DataFrame) -> pd.DataFrame:
        """
        Create derived features from raw battery telemetry.

        Expected columns:
        - cell_voltages: list of voltages [V]
        - cell_temperatures: list of temps [C]
        - soc: State of charge [0-100]
        - soh: State of health [0-100]
        - current: Pack current [A]
        - timestamp: datetime
        """
        features = pd.DataFrame()

        # Voltage statistics
        voltages = np.array(df['cell_voltages'].tolist())
        features['voltage_mean'] = voltages.mean(axis=1)
        features['voltage_std'] = voltages.std(axis=1)
        features['voltage_min'] = voltages.min(axis=1)
        features['voltage_max'] = voltages.max(axis=1)
        features['voltage_range'] = features['voltage_max'] - features['voltage_min']
        features['voltage_imbalance'] = features['voltage_max'] / (features['voltage_min'] + 1e-6)

        # Temperature statistics
        temps = np.array(df['cell_temperatures'].tolist())
        features['temp_mean'] = temps.mean(axis=1)
        features['temp_std'] = temps.std(axis=1)
        features['temp_max'] = temps.max(axis=1)
        features['thermal_gradient'] = temps.max(axis=1) - temps.min(axis=1)

        # State features
        features['soc'] = df['soc']
        features['soh'] = df['soh']
        features['current'] = df['current']
        features['power'] = df['current'] * features['voltage_mean']

        # Time-series features (requires sorted by timestamp)
        df_sorted = df.sort_values('timestamp')
        features['voltage_trend_5min'] = (
            features['voltage_mean']
            .rolling(window=30, min_periods=5)  # 30 samples @ 10s = 5min
            .apply(lambda x: np.polyfit(range(len(x)), x, 1)[0])
        )

        features['temp_volatility'] = (
            features['temp_mean']
            .rolling(window=12, min_periods=3)
            .std()
        )

        # Efficiency metrics
        features['thermal_power_ratio'] = (
            features['thermal_gradient'] / (abs(features['power']) + 1.0)
        )

        # Handle NaN from rolling windows
        features = features.fillna(method='bfill').fillna(0)

        self.feature_names = features.columns.tolist()
        return features

    def fit(self, df: pd.DataFrame) -> 'BatteryAnomalyDetector':
        """
        Train detector on normal operating data.

        Args:
            df: DataFrame with battery telemetry (normal conditions only)

        Returns:
            self for chaining
        """
        features = self.engineer_features(df)

        logger.info(f"Training on {len(features)} samples with {len(self.feature_names)} features")
        self.pipeline.fit(features)
        self.trained = True

        # Compute baseline statistics
        predictions = self.pipeline.predict(features)
        scores = self.pipeline.decision_function(features)

        logger.info(f"Training complete. Anomaly rate: {(predictions == -1).mean():.2%}")
        logger.info(f"Score range: [{scores.min():.3f}, {scores.max():.3f}]")

        return self

    def predict(
        self,
        df: pd.DataFrame,
        return_scores: bool = False
    ) -> np.ndarray | Tuple[np.ndarray, np.ndarray]:
        """
        Detect anomalies in new data.

        Args:
            df: DataFrame with battery telemetry
            return_scores: If True, return (labels, scores) tuple

        Returns:
            labels: 1 for normal, -1 for anomaly
            scores: (optional) Anomaly scores (lower = more anomalous)
        """
        if not self.trained:
            raise ValueError("Detector must be trained before prediction")

        features = self.engineer_features(df)
        labels = self.pipeline.predict(features)

        if return_scores:
            scores = self.pipeline.decision_function(features)
            return labels, scores

        return labels

    def explain_anomaly(self, df: pd.DataFrame, index: int) -> Dict[str, float]:
        """
        Explain why a specific sample was flagged as anomalous.

        Returns feature contributions sorted by absolute deviation from mean.
        """
        if not self.trained:
            raise ValueError("Detector must be trained before explanation")

        features = self.engineer_features(df)
        sample = features.iloc[index]

        # Get feature means from training (after scaling)
        scaler = self.pipeline.named_steps['scaler']
        feature_means = scaler.mean_
        feature_stds = scaler.scale_

        # Compute z-scores
        z_scores = {}
        for i, feature_name in enumerate(self.feature_names):
            z_score = abs((sample[feature_name] - feature_means[i]) / feature_stds[i])
            z_scores[feature_name] = z_score

        # Sort by deviation magnitude
        sorted_features = sorted(z_scores.items(), key=lambda x: x[1], reverse=True)

        return dict(sorted_features)

    def save(self, path: str):
        """Save trained model to disk."""
        if not self.trained:
            raise ValueError("Cannot save untrained model")

        joblib.dump({
            'pipeline': self.pipeline,
            'feature_names': self.feature_names
        }, path)
        logger.info(f"Model saved to {path}")

    @classmethod
    def load(cls, path: str) -> 'BatteryAnomalyDetector':
        """Load trained model from disk."""
        data = joblib.load(path)
        detector = cls()
        detector.pipeline = data['pipeline']
        detector.feature_names = data['feature_names']
        detector.trained = True
        logger.info(f"Model loaded from {path}")
        return detector


# Example usage
if __name__ == "__main__":
    # Load training data (normal operation only)
    train_df = pd.read_parquet('battery_telemetry_normal.parquet')

    # Train detector
    detector = BatteryAnomalyDetector(contamination=0.01)
    detector.fit(train_df)

    # Save model
    detector.save('battery_anomaly_detector.pkl')

    # Detect anomalies in new data
    test_df = pd.read_parquet('battery_telemetry_recent.parquet')
    labels, scores = detector.predict(test_df, return_scores=True)

    # Analyze anomalies
    anomaly_indices = np.where(labels == -1)[0]
    print(f"Detected {len(anomaly_indices)} anomalies ({len(anomaly_indices)/len(test_df):.2%})")

    for idx in anomaly_indices[:5]:  # Show first 5
        print(f"\nAnomaly at index {idx} (score: {scores[idx]:.3f})")
        explanation = detector.explain_anomaly(test_df, idx)
        print("Top contributing features:")
        for feature, z_score in list(explanation.items())[:5]:
            print(f"  {feature}: z-score = {z_score:.2f}")

Autoencoder for Time-Series Anomalies

python
import torch
import torch.nn as nn
import torch.optim as optim
from torch.utils.data import Dataset, DataLoader
import numpy as np
import pandas as pd
from typing import Tuple

class BatterySequenceDataset(Dataset):
    """Dataset for time-series battery data."""

    def __init__(self, sequences: np.ndarray):
        """
        Args:
            sequences: (N, seq_len, n_features) array
        """
        self.sequences = torch.FloatTensor(sequences)

    def __len__(self):
        return len(self.sequences)

    def __getitem__(self, idx):
        return self.sequences[idx]


class LSTMAutoencoder(nn.Module):
    """
    LSTM-based autoencoder for battery time-series anomaly detection.

    Architecture:
    - Encoder: LSTM -> reduces sequence to latent vector
    - Decoder: LSTM -> reconstructs sequence from latent
    - Loss: MSE between input and reconstruction
    """

    def __init__(
        self,
        n_features: int,
        seq_len: int,
        latent_dim: int = 16,
        hidden_dim: int = 64,
        n_layers: int = 2,
        dropout: float = 0.2
    ):
        super().__init__()

        self.n_features = n_features
        self.seq_len = seq_len
        self.latent_dim = latent_dim

        # Encoder: (batch, seq_len, n_features) -> (batch, latent_dim)
        self.encoder_lstm = nn.LSTM(
            input_size=n_features,
            hidden_size=hidden_dim,
            num_layers=n_layers,
            batch_first=True,
            dropout=dropout if n_layers > 1 else 0
        )
        self.encoder_fc = nn.Linear(hidden_dim, latent_dim)

        # Decoder: (batch, latent_dim) -> (batch, seq_len, n_features)
        self.decoder_fc = nn.Linear(latent_dim, hidden_dim)
        self.decoder_lstm = nn.LSTM(
            input_size=hidden_dim,
            hidden_size=hidden_dim,
            num_layers=n_layers,
            batch_first=True,
            dropout=dropout if n_layers > 1 else 0
        )
        self.output_fc = nn.Linear(hidden_dim, n_features)

    def encode(self, x: torch.Tensor) -> torch.Tensor:
        """Encode sequence to latent vector."""
        # x: (batch, seq_len, n_features)
        _, (hidden, _) = self.encoder_lstm(x)
        # hidden: (n_layers, batch, hidden_dim) -> take last layer
        latent = self.encoder_fc(hidden[-1])  # (batch, latent_dim)
        return latent

    def decode(self, latent: torch.Tensor) -> torch.Tensor:
        """Decode latent vector to sequence."""
        # latent: (batch, latent_dim)
        hidden = self.decoder_fc(latent)  # (batch, hidden_dim)

        # Repeat hidden for each time step
        hidden_seq = hidden.unsqueeze(1).repeat(1, self.seq_len, 1)
        # hidden_seq: (batch, seq_len, hidden_dim)

        decoder_out, _ = self.decoder_lstm(hidden_seq)
        reconstruction = self.output_fc(decoder_out)
        # reconstruction: (batch, seq_len, n_features)

        return reconstruction

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        """Forward pass: encode then decode."""
        latent = self.encode(x)
        reconstruction = self.decode(latent)
        return reconstruction


class AutoencoderAnomalyDetector:
    """Autoencoder-based anomaly detection for time-series data."""

    def __init__(
        self,
        n_features: int,
        seq_len: int,
        latent_dim: int = 16,
        device: str = 'cpu'
    ):
        self.model = LSTMAutoencoder(
            n_features=n_features,
            seq_len=seq_len,
            latent_dim=latent_dim
        ).to(device)

        self.device = device
        self.threshold = None
        self.scaler_mean = None
        self.scaler_std = None

    def train(
        self,
        sequences: np.ndarray,
        epochs: int = 50,
        batch_size: int = 64,
        lr: float = 0.001
    ):
        """
        Train autoencoder on normal sequences.

        Args:
            sequences: (N, seq_len, n_features) array of normal data
            epochs: Training epochs
            batch_size: Batch size
            lr: Learning rate
        """
        # Normalize
        self.scaler_mean = sequences.mean(axis=(0, 1))
        self.scaler_std = sequences.std(axis=(0, 1)) + 1e-6
        sequences_norm = (sequences - self.scaler_mean) / self.scaler_std

        # Create data loader
        dataset = BatterySequenceDataset(sequences_norm)
        loader = DataLoader(dataset, batch_size=batch_size, shuffle=True)

        # Training
        optimizer = optim.Adam(self.model.parameters(), lr=lr)
        criterion = nn.MSELoss()

        self.model.train()
        for epoch in range(epochs):
            total_loss = 0
            for batch in loader:
                batch = batch.to(self.device)

                optimizer.zero_grad()
                reconstruction = self.model(batch)
                loss = criterion(reconstruction, batch)
                loss.backward()
                optimizer.step()

                total_loss += loss.item()

            avg_loss = total_loss / len(loader)
            if (epoch + 1) % 10 == 0:
                print(f"Epoch {epoch+1}/{epochs}, Loss: {avg_loss:.6f}")

        # Compute threshold (95th percentile of training reconstruction errors)
        self.model.eval()
        with torch.no_grad():
            all_errors = []
            for batch in DataLoader(dataset, batch_size=batch_size):
                batch = batch.to(self.device)
                reconstruction = self.model(batch)
                errors = torch.mean((batch - reconstruction) ** 2, dim=(1, 2))
                all_errors.extend(errors.cpu().numpy())

            self.threshold = np.percentile(all_errors, 95)
            print(f"Anomaly threshold set to: {self.threshold:.6f}")

    def predict(self, sequences: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
        """
        Detect anomalies in sequences.

        Returns:
            labels: 0 for normal, 1 for anomaly
            errors: Reconstruction error per sequence
        """
        # Normalize
        sequences_norm = (sequences - self.scaler_mean) / self.scaler_std

        dataset = BatterySequenceDataset(sequences_norm)
        loader = DataLoader(dataset, batch_size=64)

        self.model.eval()
        all_errors = []
        with torch.no_grad():
            for batch in loader:
                batch = batch.to(self.device)
                reconstruction = self.model(batch)
                errors = torch.mean((batch - reconstruction) ** 2, dim=(1, 2))
                all_errors.extend(errors.cpu().numpy())

        errors = np.array(all_errors)
        labels = (errors > self.threshold).astype(int)

        return labels, errors

Deployment Strategy

Edge Deployment (Vehicle ECU)

yaml
# Model optimization for embedded systems
optimization:
  quantization: INT8  # Reduce model size by 4x
  pruning: 0.3        # Remove 30% least important weights
  inference_engine: ONNX Runtime
  target_latency: <50ms

hardware:
  platform: NVIDIA Jetson Xavier NX / Intel Myriad X
  memory: <100MB per model
  power_budget: <5W

monitoring:
  - Track inference latency
  - Log anomaly rate (should be <1%)
  - Alert on detector failures

Cloud Deployment (Fleet-wide)

yaml
architecture:
  ingestion: Apache Kafka
  preprocessing: Apache Spark Structured Streaming
  inference: TensorFlow Serving / Ray Serve
  storage: TimescaleDB (anomalies) + S3 (raw data)

scalability:
  - Horizontal scaling based on vehicle count
  - Batch processing for non-real-time analysis
  - A/B testing for model updates

mlops:
  - Model versioning with MLflow
  - Automated retraining on new anomaly patterns
  - Canary deployments (5% -> 25% -> 100%)
  - Rollback on performance degradation

Evaluation Metrics

For Labeled Test Sets

python
from sklearn.metrics import (
    precision_score, recall_score, f1_score,
    roc_auc_score, confusion_matrix, classification_report
)

def evaluate_detector(y_true: np.ndarray, y_pred: np.ndarray, scores: np.ndarray):
    """
    Comprehensive evaluation of anomaly detector.

    Args:
        y_true: Ground truth labels (1=anomaly, 0=normal)
        y_pred: Predicted labels (1=anomaly, 0=normal)
        scores: Anomaly scores (continuous)
    """
    # Convert Isolation Forest labels (-1, 1) to (1, 0)
    y_pred_binary = (y_pred == -1).astype(int)

    # Metrics
    precision = precision_score(y_true, y_pred_binary)
    recall = recall_score(y_true, y_pred_binary)
    f1 = f1_score(y_true, y_pred_binary)

    # AUC-ROC (requires scores)
    # Negate scores since Isolation Forest: lower = more anomalous
    auc = roc_auc_score(y_true, -scores)

    print(f"Precision: {precision:.3f}")
    print(f"Recall: {recall:.3f}")
    print(f"F1-Score: {f1:.3f}")
    print(f"AUC-ROC: {auc:.3f}")
    print("\nConfusion Matrix:")
    print(confusion_matrix(y_true, y_pred_binary))
    print("\nClassification Report:")
    print(classification_report(y_true, y_pred_binary,
                                target_names=['Normal', 'Anomaly']))

For Unlabeled Production Data

python
def monitor_detector_health(predictions: np.ndarray, scores: np.ndarray):
    """
    Monitor detector behavior in production (no ground truth).

    Alerts:
    - Anomaly rate too high (>5%) -> possible drift
    - Anomaly rate too low (<0.1%) -> detector not sensitive
    - Score distribution shift -> model degradation
    """
    anomaly_rate = (predictions == -1).mean()
    score_mean = scores.mean()
    score_std = scores.std()

    print(f"Anomaly Rate: {anomaly_rate:.2%}")
    print(f"Score Mean: {score_mean:.3f} +/- {score_std:.3f}")

    # Alerts
    if anomaly_rate > 0.05:
        print("WARNING: High anomaly rate detected. Possible data drift.")
    elif anomaly_rate < 0.001:
        print("WARNING: Low anomaly rate. Detector may not be sensitive enough.")

Production Checklist

• ] Model trained on representative normal data (>10k samples)
• ] Threshold calibrated on validation set (target: 1-2% false positive rate)
• ] Feature engineering validated for edge cases (missing sensors, extreme temps)
• ] Model compressed for edge deployment (ONNX + quantization)
• ] Inference latency tested (<100ms p99)
• ] Monitoring dashboards configured (Grafana + Prometheus)
• ] Alert thresholds set (anomaly rate, detector health)
• ] Retraining pipeline automated (weekly batch on new normal data)
• ] A/B testing framework ready
• ] Rollback procedure documented

Driver Behavior Analysis

Driver Behavior Analysis and Safety Scoring

Analyze driving patterns, compute safety scores, and provide personalized feedback for driver improvement. Focus on risk assessment, insurance optimization, and training recommendations.

Use Cases

1. Safety Scoring: Quantitative assessment of driver risk (0-100 scale)
2. Insurance Optimization: Usage-based insurance (UBI) premiums
3. Driver Training: Personalized improvement recommendations
4. Fleet Management: Identify high-risk drivers for intervention
5. Accident Prevention: Predict and prevent risky behaviors

Feature Engineering

python
import pandas as pd
import numpy as np
from typing import Dict, Tuple
from scipy.stats import zscore

class DriverBehaviorFeatureEngineer:
    """
    Extract driver behavior features from vehicle telemetry.

    Feature categories:
    - Acceleration patterns (smooth vs aggressive)
    - Braking behavior (harsh braking frequency)
    - Speed management (speeding, speed variance)
    - Cornering (lateral acceleration)
    - Anticipation (time-to-collision events)
    - Distraction indicators (steering wheel anomalies)
    """

    @staticmethod
    def extract_trip_features(trip_data: pd.DataFrame) -> Dict[str, float]:
        """
        Extract behavior features from a single trip.

        Expected columns:
        - timestamp, speed_kmh, acceleration_ms2, lateral_accel_ms2,
        - brake_pressure, steering_angle, distance_to_vehicle_m,
        - speed_limit_kmh (if available)
        """
        features = {}

        # Trip metadata
        features['duration_min'] = (
            (trip_data['timestamp'].max() - trip_data['timestamp'].min())
            .total_seconds() / 60
        )
        features['distance_km'] = trip_data['speed_kmh'].sum() / 3600  # Simplified

        # Acceleration analysis
        accel = trip_data['acceleration_ms2']
        features['accel_mean'] = accel.mean()
        features['accel_std'] = accel.std()
        features['accel_max'] = accel.max()

        # Harsh acceleration events (> 2.5 m/s²)
        harsh_accel = accel > 2.5
        features['harsh_accel_count'] = harsh_accel.sum()
        features['harsh_accel_rate_per_100km'] = (
            features['harsh_accel_count'] / (features['distance_km'] + 1e-3) * 100
        )

        # Braking analysis
        # Assume deceleration is negative acceleration
        decel = accel[accel < 0]
        features['decel_mean'] = decel.mean() if len(decel) > 0 else 0
        features['decel_std'] = decel.std() if len(decel) > 0 else 0

        # Harsh braking events (< -3.0 m/s²)
        harsh_brake = accel < -3.0
        features['harsh_brake_count'] = harsh_brake.sum()
        features['harsh_brake_rate_per_100km'] = (
            features['harsh_brake_count'] / (features['distance_km'] + 1e-3) * 100
        )

        # Brake pressure statistics
        if 'brake_pressure' in trip_data.columns:
            brake = trip_data['brake_pressure']
            features['brake_pressure_mean'] = brake.mean()
            features['brake_pressure_max'] = brake.max()

            # Emergency braking (brake pressure > 0.8)
            emergency_brake = brake > 0.8
            features['emergency_brake_count'] = emergency_brake.sum()

        # Speed management
        speed = trip_data['speed_kmh']
        features['speed_mean'] = speed.mean()
        features['speed_std'] = speed.std()
        features['speed_max'] = speed.max()

        # Speed variance (indication of smooth vs erratic driving)
        features['speed_variance'] = speed.var()

        # Speeding analysis
        if 'speed_limit_kmh' in trip_data.columns:
            speeding = speed > trip_data['speed_limit_kmh']
            features['speeding_pct'] = (speeding.sum() / len(trip_data)) * 100

            # Excessive speeding (>20 km/h over limit)
            excessive_speeding = (speed - trip_data['speed_limit_kmh']) > 20
            features['excessive_speeding_pct'] = (
                excessive_speeding.sum() / len(trip_data)
            ) * 100
        else:
            features['speeding_pct'] = 0
            features['excessive_speeding_pct'] = 0

        # Cornering analysis
        if 'lateral_accel_ms2' in trip_data.columns:
            lateral_accel = trip_data['lateral_accel_ms2']
            features['lateral_accel_mean'] = abs(lateral_accel).mean()
            features['lateral_accel_max'] = abs(lateral_accel).max()

            # Harsh cornering (lateral accel > 4 m/s²)
            harsh_corner = abs(lateral_accel) > 4.0
            features['harsh_corner_count'] = harsh_corner.sum()
            features['harsh_corner_rate_per_100km'] = (
                features['harsh_corner_count'] / (features['distance_km'] + 1e-3) * 100
            )

        # Anticipation and safety margins
        if 'distance_to_vehicle_m' in trip_data.columns:
            # Time-to-collision (TTC) at current speed
            # TTC = distance / speed
            ttc = trip_data['distance_to_vehicle_m'] / (
                trip_data['speed_kmh'] / 3.6 + 1e-3
            )  # Convert to m/s

            # Critical TTC events (<2 seconds)
            critical_ttc = ttc < 2.0
            features['critical_ttc_count'] = critical_ttc.sum()
            features['critical_ttc_rate_per_100km'] = (
                features['critical_ttc_count'] / (features['distance_km'] + 1e-3) * 100
            )

            features['avg_following_distance_m'] = trip_data['distance_to_vehicle_m'].mean()

        # Steering behavior (distraction indicator)
        if 'steering_angle' in trip_data.columns:
            steering = trip_data['steering_angle']
            features['steering_variance'] = steering.diff().abs().var()

            # Erratic steering (high frequency oscillations)
            steering_changes = abs(steering.diff())
            features['erratic_steering_events'] = (steering_changes > 10).sum()

        # Idle time (inefficient driving)
        idle = (speed < 1) & (accel.abs() < 0.1)
        features['idle_time_min'] = idle.sum() / 60  # Assuming 1 Hz data

        return features

    @staticmethod
    def aggregate_driver_features(
        trip_features_list: list
    ) -> pd.DataFrame:
        """
        Aggregate features across multiple trips for driver profiling.

        Args:
            trip_features_list: List of feature dictionaries (one per trip)

        Returns:
            DataFrame with aggregated driver statistics
        """
        df = pd.DataFrame(trip_features_list)

        # Aggregate statistics
        agg_features = {}

        # Mean metrics
        mean_cols = [
            'harsh_accel_rate_per_100km',
            'harsh_brake_rate_per_100km',
            'harsh_corner_rate_per_100km',
            'speeding_pct',
            'excessive_speeding_pct',
            'speed_variance',
            'critical_ttc_rate_per_100km'
        ]

        for col in mean_cols:
            if col in df.columns:
                agg_features[f'{col}_mean'] = df[col].mean()
                agg_features[f'{col}_std'] = df[col].std()

        # Count metrics
        count_cols = [
            'harsh_accel_count',
            'harsh_brake_count',
            'emergency_brake_count',
            'critical_ttc_count'
        ]

        for col in count_cols:
            if col in df.columns:
                agg_features[f'{col}_total'] = df[col].sum()

        # Totals
        agg_features['total_trips'] = len(df)
        agg_features['total_distance_km'] = df['distance_km'].sum()
        agg_features['total_duration_min'] = df['duration_min'].sum()

        return pd.DataFrame([agg_features])


# Example usage
if __name__ == "__main__":
    # Load trip data
    trip_df = pd.read_parquet('trip_001.parquet')

    # Extract features
    engineer = DriverBehaviorFeatureEngineer()
    features = engineer.extract_trip_features(trip_df)

    print("Trip Features:")
    for key, value in features.items():
        print(f"  {key}: {value:.3f}")

Safety Scoring Model

python
import numpy as np
import pandas as pd
from sklearn.preprocessing import MinMaxScaler
from typing import Dict, Tuple
import joblib

class DriverSafetyScoringModel:
    """
    Compute comprehensive driver safety score (0-100).

    Methodology:
    - Weighted composite of sub-scores
    - Penalize risky behaviors
    - Reward safe driving patterns
    - Calibrated against industry benchmarks
    """

    def __init__(self):
        """
        Initialize scoring model with component weights.
        """
        # Component weights (must sum to 1.0)
        self.weights = {
            'acceleration': 0.20,
            'braking': 0.25,
            'speed': 0.25,
            'cornering': 0.15,
            'anticipation': 0.15
        }

        # Benchmark values (50th percentile of safe drivers)
        self.benchmarks = {
            'harsh_accel_rate_per_100km': 5.0,
            'harsh_brake_rate_per_100km': 3.0,
            'harsh_corner_rate_per_100km': 2.0,
            'speeding_pct': 10.0,
            'excessive_speeding_pct': 1.0,
            'speed_variance': 150.0,
            'critical_ttc_rate_per_100km': 2.0
        }

        # Penalty factors (multiplier for scores below threshold)
        self.penalties = {
            'emergency_brake_count': 5.0,  # -5 points per event
            'critical_ttc_count': 3.0       # -3 points per event
        }

    def compute_acceleration_score(self, features: Dict) -> float:
        """
        Score acceleration behavior (0-100).

        Lower harsh acceleration rate = higher score.
        """
        rate = features.get('harsh_accel_rate_per_100km_mean', 0)
        benchmark = self.benchmarks['harsh_accel_rate_per_100km']

        # Exponential decay: score = 100 * exp(-rate / benchmark)
        score = 100 * np.exp(-rate / benchmark)

        return np.clip(score, 0, 100)

    def compute_braking_score(self, features: Dict) -> float:
        """
        Score braking behavior (0-100).

        Penalize harsh braking and emergency stops.
        """
        rate = features.get('harsh_brake_rate_per_100km_mean', 0)
        benchmark = self.benchmarks['harsh_brake_rate_per_100km']

        base_score = 100 * np.exp(-rate / benchmark)

        # Apply penalty for emergency braking
        emergency_count = features.get('emergency_brake_count_total', 0)
        penalty = emergency_count * self.penalties['emergency_brake_count']

        score = base_score - penalty

        return np.clip(score, 0, 100)

    def compute_speed_score(self, features: Dict) -> float:
        """
        Score speed management (0-100).

        Penalize speeding and erratic speed changes.
        """
        speeding_pct = features.get('speeding_pct_mean', 0)
        excessive_pct = features.get('excessive_speeding_pct_mean', 0)
        speed_var = features.get('speed_variance_mean', 0)

        # Speed limit compliance
        speeding_benchmark = self.benchmarks['speeding_pct']
        speeding_score = 100 * np.exp(-speeding_pct / speeding_benchmark)

        # Excessive speeding (heavy penalty)
        excessive_benchmark = self.benchmarks['excessive_speeding_pct']
        excessive_penalty = 50 * (excessive_pct / (excessive_benchmark + 1))

        # Speed smoothness
        variance_benchmark = self.benchmarks['speed_variance']
        smoothness_score = 100 * np.exp(-speed_var / variance_benchmark)

        # Weighted average
        score = (
            0.5 * speeding_score +
            0.3 * smoothness_score -
            0.2 * excessive_penalty
        )

        return np.clip(score, 0, 100)

    def compute_cornering_score(self, features: Dict) -> float:
        """
        Score cornering behavior (0-100).

        Reward smooth cornering, penalize harsh turns.
        """
        rate = features.get('harsh_corner_rate_per_100km_mean', 0)
        benchmark = self.benchmarks['harsh_corner_rate_per_100km']

        score = 100 * np.exp(-rate / benchmark)

        return np.clip(score, 0, 100)

    def compute_anticipation_score(self, features: Dict) -> float:
        """
        Score anticipation and safety margins (0-100).

        Penalize tailgating and critical TTC events.
        """
        ttc_rate = features.get('critical_ttc_rate_per_100km_mean', 0)
        benchmark = self.benchmarks['critical_ttc_rate_per_100km']

        base_score = 100 * np.exp(-ttc_rate / benchmark)

        # Apply penalty for critical TTC events
        ttc_count = features.get('critical_ttc_count_total', 0)
        penalty = ttc_count * self.penalties['critical_ttc_count']

        score = base_score - penalty

        return np.clip(score, 0, 100)

    def compute_overall_score(self, features: Dict) -> Tuple[float, Dict[str, float]]:
        """
        Compute weighted overall safety score.

        Returns:
            overall_score: Composite score (0-100)
            component_scores: Dictionary of sub-scores
        """
        # Compute component scores
        components = {
            'acceleration': self.compute_acceleration_score(features),
            'braking': self.compute_braking_score(features),
            'speed': self.compute_speed_score(features),
            'cornering': self.compute_cornering_score(features),
            'anticipation': self.compute_anticipation_score(features)
        }

        # Weighted average
        overall = sum(
            components[key] * self.weights[key]
            for key in components
        )

        return overall, components

    def classify_risk_level(self, score: float) -> str:
        """
        Classify driver into risk category.

        Args:
            score: Safety score (0-100)

        Returns:
            Risk category string
        """
        if score >= 85:
            return 'Low Risk'
        elif score >= 70:
            return 'Medium Risk'
        elif score >= 50:
            return 'High Risk'
        else:
            return 'Critical Risk'

    def generate_feedback(
        self,
        components: Dict[str, float]
    ) -> list:
        """
        Generate personalized feedback for driver improvement.

        Args:
            components: Component scores dictionary

        Returns:
            List of improvement recommendations
        """
        feedback = []

        # Identify weakest areas
        sorted_components = sorted(components.items(), key=lambda x: x[1])

        for component, score in sorted_components[:2]:  # Focus on 2 worst
            if score < 70:
                if component == 'acceleration':
                    feedback.append(
                        "⚠️ Reduce harsh acceleration. Gradually increase speed to "
                        "improve fuel efficiency and safety."
                    )
                elif component == 'braking':
                    feedback.append(
                        "⚠️ Anticipate stops earlier to avoid harsh braking. "
                        "Maintain safe following distance."
                    )
                elif component == 'speed':
                    feedback.append(
                        "⚠️ Adhere to speed limits and maintain consistent speed. "
                        "Reduce speeding incidents."
                    )
                elif component == 'cornering':
                    feedback.append(
                        "⚠️ Slow down before corners. Smooth steering reduces tire wear "
                        "and improves stability."
                    )
                elif component == 'anticipation':
                    feedback.append(
                        "⚠️ Increase following distance. Avoid tailgating to reduce "
                        "collision risk."
                    )

        # Positive reinforcement for strong areas
        best_component = max(components.items(), key=lambda x: x[1])
        if best_component[1] >= 90:
            feedback.append(
                f"✅ Excellent {best_component[0]} behavior! Keep it up."
            )

        return feedback


# Example usage
if __name__ == "__main__":
    # Sample driver features (aggregated from multiple trips)
    driver_features = {
        'harsh_accel_rate_per_100km_mean': 7.2,
        'harsh_brake_rate_per_100km_mean': 4.5,
        'harsh_corner_rate_per_100km_mean': 1.8,
        'speeding_pct_mean': 15.0,
        'excessive_speeding_pct_mean': 2.5,
        'speed_variance_mean': 180.0,
        'critical_ttc_rate_per_100km_mean': 3.0,
        'emergency_brake_count_total': 2,
        'critical_ttc_count_total': 5
    }

    # Compute safety score
    scorer = DriverSafetyScoringModel()
    overall_score, components = scorer.compute_overall_score(driver_features)
    risk_level = scorer.classify_risk_level(overall_score)

    print(f"\nOverall Safety Score: {overall_score:.1f} / 100")
    print(f"Risk Level: {risk_level}\n")

    print("Component Scores:")
    for component, score in components.items():
        print(f"  {component.capitalize()}: {score:.1f}")

    print("\nPersonalized Feedback:")
    feedback = scorer.generate_feedback(components)
    for item in feedback:
        print(f"  {item}")

Driver Clustering

python
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

class DriverClustering:
    """
    Cluster drivers into behavioral profiles.

    Profiles:
    - Conservative: Low speeds, smooth acceleration/braking
    - Aggressive: High speeds, harsh maneuvers
    - Efficient: Optimized for fuel/energy efficiency
    - Distracted: Erratic patterns, low anticipation
    """

    def __init__(self, n_clusters: int = 4):
        self.n_clusters = n_clusters
        self.scaler = StandardScaler()
        self.kmeans = KMeans(n_clusters=n_clusters, random_state=42)

    def prepare_features(self, driver_data: pd.DataFrame) -> pd.DataFrame:
        """
        Select and normalize features for clustering.

        Expected columns: driver_id, harsh_accel_rate, harsh_brake_rate,
        speeding_pct, speed_variance, critical_ttc_rate
        """
        feature_cols = [
            'harsh_accel_rate_per_100km_mean',
            'harsh_brake_rate_per_100km_mean',
            'speeding_pct_mean',
            'speed_variance_mean',
            'critical_ttc_rate_per_100km_mean'
        ]

        X = driver_data[feature_cols]
        X_scaled = self.scaler.fit_transform(X)

        return X_scaled, feature_cols

    def fit_predict(self, driver_data: pd.DataFrame) -> np.ndarray:
        """
        Cluster drivers and return labels.

        Returns:
            Cluster labels (0 to n_clusters-1)
        """
        X_scaled, _ = self.prepare_features(driver_data)
        labels = self.kmeans.fit_predict(X_scaled)

        return labels

    def visualize_clusters(
        self,
        driver_data: pd.DataFrame,
        labels: np.ndarray
    ):
        """
        Visualize driver clusters using pair plot.
        """
        # Add cluster labels
        driver_data_labeled = driver_data.copy()
        driver_data_labeled['cluster'] = labels

        # Select key features for visualization
        viz_features = [
            'harsh_accel_rate_per_100km_mean',
            'harsh_brake_rate_per_100km_mean',
            'speeding_pct_mean',
            'cluster'
        ]

        sns.pairplot(
            driver_data_labeled[viz_features],
            hue='cluster',
            palette='viridis',
            diag_kind='kde'
        )
        plt.suptitle('Driver Behavior Clusters', y=1.02)
        plt.show()

    def describe_clusters(
        self,
        driver_data: pd.DataFrame,
        labels: np.ndarray
    ) -> pd.DataFrame:
        """
        Describe cluster characteristics.

        Returns:
            DataFrame with cluster profiles
        """
        driver_data_labeled = driver_data.copy()
        driver_data_labeled['cluster'] = labels

        profiles = driver_data_labeled.groupby('cluster').agg({
            'harsh_accel_rate_per_100km_mean': 'mean',
            'harsh_brake_rate_per_100km_mean': 'mean',
            'speeding_pct_mean': 'mean',
            'speed_variance_mean': 'mean',
            'critical_ttc_rate_per_100km_mean': 'mean'
        }).round(2)

        return profiles

    def assign_profile_names(
        self,
        profiles: pd.DataFrame
    ) -> Dict[int, str]:
        """
        Assign interpretable names to clusters.
        """
        names = {}

        for cluster_id in profiles.index:
            accel = profiles.loc[cluster_id, 'harsh_accel_rate_per_100km_mean']
            brake = profiles.loc[cluster_id, 'harsh_brake_rate_per_100km_mean']
            speeding = profiles.loc[cluster_id, 'speeding_pct_mean']

            # Rule-based naming
            if accel < 3 and brake < 2 and speeding < 5:
                names[cluster_id] = "Conservative"
            elif accel > 8 and brake > 5:
                names[cluster_id] = "Aggressive"
            elif speeding > 20:
                names[cluster_id] = "Speeder"
            else:
                names[cluster_id] = "Average"

        return names

Deployment Architecture

yaml
# Driver behavior analytics pipeline
pipeline:
  data_collection:
    source: Vehicle telemetry (CAN bus) + GPS
    frequency: 10 Hz (acceleration, speed), 1 Hz (location)
    storage: TimescaleDB (raw) + Parquet (trips)

  trip_segmentation:
    algorithm: Speed-based (ignition on/off or idle detection)
    min_duration: 2 minutes
    min_distance: 1 km

  feature_extraction:
    engine: Apache Spark (batch) or Python (streaming)
    frequency: Per-trip (real-time) + weekly aggregation
    storage: PostgreSQL (features table)

  scoring:
    model: Rule-based composite score
    update_frequency: Daily (driver profile)
    storage: PostgreSQL (scores table)

  dashboard:
    framework: Streamlit / React
    update_frequency: Real-time (trip-level), daily (aggregates)
    features: Score trends, feedback, leaderboard

monitoring:
  metrics:
    - Score distribution (fleet-wide)
    - Score volatility (per driver)
    - Feature extraction success rate
    - Dashboard latency

alerts:
  - Critical risk driver detected (score < 50)
  - Repeated risky behaviors (3+ harsh events per trip)
  - Score degradation (>10 point drop in 7 days)

Production Checklist

• ] Trip segmentation validated (accuracy > 95%)
• ] Feature extraction handles edge cases (short trips, missing sensors)
• ] Score calibration validated against insurance claims data
• ] Feedback messages reviewed for clarity and positivity
• ] Dashboard responsive (<2s load time)
• ] Privacy compliance (GDPR, CCPA) - anonymization
• ] Driver consent workflow implemented
• ] Alert thresholds configured (avoid alert fatigue)
• ] A/B testing for scoring algorithm changes
• ] Driver appeal process documented

Energy Optimization

Energy Optimization and Route Planning

Optimize energy consumption through intelligent route planning, charging strategies, and predictive algorithms. Focus on electric vehicles with battery-electric and plug-in hybrid powertrains.

Use Cases

1. Route Optimization: Minimize energy consumption for given destination
2. Charging Strategy: Optimize charging times and locations for cost/time
3. Range Prediction: Accurate remaining range estimation
4. Eco-Routing: Balance time vs energy trade-offs
5. Fleet Electrification: Optimal EV deployment planning

Route Optimization with Reinforcement Learning

python
import numpy as np
import pandas as pd
import gym
from gym import spaces
from typing import Tuple, Dict, List
import torch
import torch.nn as nn
import torch.optim as optim
from collections import deque
import random

class EVRoutingEnvironment(gym.Env):
    """
    OpenAI Gym environment for EV route optimization.

    State: Current location, SOC, traffic conditions, weather
    Action: Next waypoint selection
    Reward: -(energy_consumed + time_penalty + charging_cost)
    """

    def __init__(
        self,
        road_network: pd.DataFrame,
        charging_stations: pd.DataFrame,
        battery_capacity_kwh: float = 75.0
    ):
        """
        Args:
            road_network: Graph with edges (src, dst, distance, elevation, speed_limit)
            charging_stations: Locations with (lat, lon, power_kw, cost_per_kwh)
            battery_capacity_kwh: Vehicle battery capacity
        """
        super().__init__()

        self.road_network = road_network
        self.charging_stations = charging_stations
        self.battery_capacity = battery_capacity_kwh

        # State space: [current_node_id, soc, time_of_day, traffic_level]
        self.observation_space = spaces.Box(
            low=np.array([0, 0, 0, 0]),
            high=np.array([len(road_network), 100, 24, 10]),
            dtype=np.float32
        )

        # Action space: [next_node_id, charge_decision (0=no, 1=yes)]
        self.action_space = spaces.MultiDiscrete([len(road_network), 2])

        self.reset()

    def reset(self) -> np.ndarray:
        """
        Reset environment to initial state.

        Returns:
            Initial observation
        """
        self.current_node = 0  # Start node
        self.destination_node = len(self.road_network) - 1
        self.soc = 80.0  # Start with 80% SOC
        self.time = 8.0  # Start at 8 AM
        self.total_energy = 0
        self.total_cost = 0
        self.trajectory = [self.current_node]

        return self._get_observation()

    def _get_observation(self) -> np.ndarray:
        """Get current state observation."""
        traffic_level = self._get_traffic_level(self.time)
        return np.array([
            self.current_node,
            self.soc,
            self.time,
            traffic_level
        ], dtype=np.float32)

    def _get_traffic_level(self, time: float) -> float:
        """
        Model traffic congestion by time of day.

        Peak hours: 7-9 AM, 5-7 PM -> High traffic
        """
        if (7 <= time <= 9) or (17 <= time <= 19):
            return 8.0  # High traffic
        elif (6 <= time <= 10) or (16 <= time <= 20):
            return 5.0  # Medium traffic
        else:
            return 2.0  # Low traffic

    def _compute_energy_consumption(
        self,
        distance_km: float,
        elevation_gain_m: float,
        speed_kmh: float,
        traffic_level: float
    ) -> float:
        """
        Estimate energy consumption for segment.

        Model: Base consumption + elevation + traffic penalty
        """
        # Base consumption (kWh per 100 km)
        base_consumption = 18.0

        # Elevation impact (100m climb ≈ 1 kWh)
        elevation_energy = max(0, elevation_gain_m / 100)

        # Speed efficiency (optimal around 60 km/h)
        speed_factor = 1 + 0.01 * abs(speed_kmh - 60)

        # Traffic penalty (stop-and-go increases consumption)
        traffic_factor = 1 + 0.05 * traffic_level

        total_energy = (
            base_consumption * distance_km / 100 * speed_factor * traffic_factor +
            elevation_energy
        )

        return total_energy

    def step(self, action: Tuple[int, int]) -> Tuple[np.ndarray, float, bool, Dict]:
        """
        Execute action and return next state.

        Args:
            action: (next_node_id, charge_decision)

        Returns:
            observation, reward, done, info
        """
        next_node, charge_decision = action

        # Get edge info
        edge = self.road_network[
            (self.road_network['src'] == self.current_node) &
            (self.road_network['dst'] == next_node)
        ]

        if edge.empty:
            # Invalid action (no edge)
            return self._get_observation(), -1000, True, {'error': 'invalid_edge'}

        edge = edge.iloc[0]
        distance = edge['distance_km']
        elevation = edge['elevation_gain_m']
        speed_limit = edge['speed_limit_kmh']

        # Traffic-adjusted speed
        traffic_level = self._get_traffic_level(self.time)
        avg_speed = speed_limit * (1 - 0.05 * traffic_level)  # Slower in traffic

        # Energy consumption
        energy_consumed = self._compute_energy_consumption(
            distance, elevation, avg_speed, traffic_level
        )

        # Update SOC
        self.soc -= (energy_consumed / self.battery_capacity) * 100

        # Time elapsed
        time_elapsed = distance / avg_speed  # hours
        self.time += time_elapsed

        # Charging decision
        charging_cost = 0
        charging_time = 0

        if charge_decision == 1:
            # Find nearest charging station
            station = self.charging_stations.iloc[0]  # Simplified
            charge_amount_kwh = (80 - self.soc) / 100 * self.battery_capacity
            charging_time = charge_amount_kwh / station['power_kw']
            charging_cost = charge_amount_kwh * station['cost_per_kwh']

            self.soc = 80.0  # Charge to 80%
            self.time += charging_time
            self.total_cost += charging_cost

        # Update trajectory
        self.current_node = next_node
        self.trajectory.append(next_node)
        self.total_energy += energy_consumed

        # Check if done
        done = False
        if self.current_node == self.destination_node:
            done = True
        elif self.soc < 10:
            # Ran out of battery
            done = True
            reward = -1000
            return self._get_observation(), reward, done, {'error': 'battery_depleted'}

        # Compute reward
        # Minimize: energy + time + charging cost
        reward = -(
            energy_consumed +
            10 * time_elapsed +  # Time penalty
            charging_cost
        )

        # Bonus for reaching destination
        if done and self.current_node == self.destination_node:
            reward += 100

        info = {
            'energy_consumed': energy_consumed,
            'time_elapsed': time_elapsed,
            'charging_cost': charging_cost,
            'soc': self.soc
        }

        return self._get_observation(), reward, done, info


class DQN(nn.Module):
    """
    Deep Q-Network for route optimization.

    Architecture: Fully connected layers with ReLU activations.
    """

    def __init__(self, state_dim: int, action_dim: int):
        super().__init__()

        self.fc = nn.Sequential(
            nn.Linear(state_dim, 128),
            nn.ReLU(),
            nn.Linear(128, 128),
            nn.ReLU(),
            nn.Linear(128, action_dim)
        )

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        return self.fc(x)


class DQNAgent:
    """
    DQN agent for EV route optimization.

    Uses experience replay and target network for stable learning.
    """

    def __init__(
        self,
        state_dim: int,
        action_dim: int,
        lr: float = 0.001,
        gamma: float = 0.99,
        epsilon: float = 1.0,
        epsilon_decay: float = 0.995,
        epsilon_min: float = 0.01
    ):
        self.state_dim = state_dim
        self.action_dim = action_dim
        self.gamma = gamma
        self.epsilon = epsilon
        self.epsilon_decay = epsilon_decay
        self.epsilon_min = epsilon_min

        # Q-network and target network
        self.q_network = DQN(state_dim, action_dim)
        self.target_network = DQN(state_dim, action_dim)
        self.target_network.load_state_dict(self.q_network.state_dict())

        self.optimizer = optim.Adam(self.q_network.parameters(), lr=lr)
        self.criterion = nn.MSELoss()

        # Experience replay
        self.memory = deque(maxlen=10000)
        self.batch_size = 64

    def select_action(self, state: np.ndarray) -> int:
        """
        Epsilon-greedy action selection.
        """
        if random.random() < self.epsilon:
            return random.randint(0, self.action_dim - 1)
        else:
            state_tensor = torch.FloatTensor(state).unsqueeze(0)
            with torch.no_grad():
                q_values = self.q_network(state_tensor)
            return q_values.argmax().item()

    def store_transition(
        self,
        state: np.ndarray,
        action: int,
        reward: float,
        next_state: np.ndarray,
        done: bool
    ):
        """Store experience in replay buffer."""
        self.memory.append((state, action, reward, next_state, done))

    def train(self):
        """Train Q-network on batch from replay buffer."""
        if len(self.memory) < self.batch_size:
            return

        # Sample batch
        batch = random.sample(self.memory, self.batch_size)
        states, actions, rewards, next_states, dones = zip(*batch)

        states = torch.FloatTensor(np.array(states))
        actions = torch.LongTensor(actions)
        rewards = torch.FloatTensor(rewards)
        next_states = torch.FloatTensor(np.array(next_states))
        dones = torch.FloatTensor(dones)

        # Current Q-values
        q_values = self.q_network(states).gather(1, actions.unsqueeze(1)).squeeze()

        # Target Q-values
        with torch.no_grad():
            next_q_values = self.target_network(next_states).max(1)[0]
            target_q_values = rewards + self.gamma * next_q_values * (1 - dones)

        # Loss and optimization
        loss = self.criterion(q_values, target_q_values)

        self.optimizer.zero_grad()
        loss.backward()
        self.optimizer.step()

        # Decay epsilon
        self.epsilon = max(self.epsilon_min, self.epsilon * self.epsilon_decay)

    def update_target_network(self):
        """Copy weights from Q-network to target network."""
        self.target_network.load_state_dict(self.q_network.state_dict())


# Example training loop
if __name__ == "__main__":
    # Create environment (with dummy data)
    road_network = pd.DataFrame({
        'src': [0, 1, 2],
        'dst': [1, 2, 3],
        'distance_km': [10, 15, 20],
        'elevation_gain_m': [50, 100, 150],
        'speed_limit_kmh': [60, 80, 100]
    })

    charging_stations = pd.DataFrame({
        'lat': [37.7749],
        'lon': [-122.4194],
        'power_kw': [150],
        'cost_per_kwh': [0.25]
    })

    env = EVRoutingEnvironment(road_network, charging_stations)

    # Initialize agent
    agent = DQNAgent(
        state_dim=env.observation_space.shape[0],
        action_dim=len(road_network)  # Simplified
    )

    # Training
    episodes = 500
    for episode in range(episodes):
        state = env.reset()
        total_reward = 0

        for step in range(100):
            action = agent.select_action(state)
            # Simplified action (no charging decision for now)
            next_state, reward, done, info = env.step((action, 0))

            agent.store_transition(state, action, reward, next_state, done)
            agent.train()

            state = next_state
            total_reward += reward

            if done:
                break

        # Update target network every 10 episodes
        if episode % 10 == 0:
            agent.update_target_network()

        if episode % 50 == 0:
            print(f"Episode {episode}, Total Reward: {total_reward:.2f}, "
                  f"Epsilon: {agent.epsilon:.3f}")

Charging Strategy Optimization

python
import pulp
import numpy as np
import pandas as pd
from typing import List, Dict, Tuple

class ChargingStrategyOptimizer:
    """
    Optimize charging schedule for fleet to minimize cost and maximize availability.

    Objective: Minimize electricity cost while ensuring all vehicles charged for next day.

    Constraints:
    - Each vehicle must reach target SOC before departure
    - Charger capacity limits (power, number of chargers)
    - Time-of-use electricity pricing
    """

    def __init__(
        self,
        n_vehicles: int,
        n_chargers: int,
        charger_power_kw: float,
        time_slots: int = 24,
        battery_capacity_kwh: float = 75.0
    ):
        """
        Args:
            n_vehicles: Number of vehicles in fleet
            n_chargers: Number of available chargers
            charger_power_kw: Charging power per charger [kW]
            time_slots: Planning horizon (hours)
            battery_capacity_kwh: Vehicle battery capacity
        """
        self.n_vehicles = n_vehicles
        self.n_chargers = n_chargers
        self.charger_power = charger_power_kw
        self.time_slots = time_slots
        self.battery_capacity = battery_capacity_kwh

    def optimize_schedule(
        self,
        arrival_soc: np.ndarray,
        target_soc: np.ndarray,
        arrival_times: np.ndarray,
        departure_times: np.ndarray,
        electricity_prices: np.ndarray
    ) -> Dict:
        """
        Optimize charging schedule using linear programming.

        Args:
            arrival_soc: Current SOC for each vehicle [%] (n_vehicles,)
            target_soc: Desired SOC for each vehicle [%] (n_vehicles,)
            arrival_times: Arrival time slot for each vehicle (n_vehicles,)
            departure_times: Departure time slot for each vehicle (n_vehicles,)
            electricity_prices: Price per kWh for each time slot [$/kWh] (time_slots,)

        Returns:
            Dictionary with optimal charging schedule
        """
        # Create optimization problem
        prob = pulp.LpProblem("Fleet_Charging_Optimization", pulp.LpMinimize)

        # Decision variables: charging[vehicle, time] = kWh charged in time slot
        charging = {}
        for v in range(self.n_vehicles):
            for t in range(self.time_slots):
                charging[v, t] = pulp.LpVariable(
                    f"charge_v{v}_t{t}",
                    lowBound=0,
                    upBound=self.charger_power  # Max charge rate
                )

        # Binary variable: is_charging[vehicle, time] = 1 if vehicle is charging
        is_charging = {}
        for v in range(self.n_vehicles):
            for t in range(self.time_slots):
                is_charging[v, t] = pulp.LpVariable(
                    f"is_charging_v{v}_t{t}",
                    cat='Binary'
                )

        # Objective: Minimize total electricity cost
        total_cost = pulp.lpSum([
            charging[v, t] * electricity_prices[t]
            for v in range(self.n_vehicles)
            for t in range(self.time_slots)
        ])
        prob += total_cost

        # Constraint 1: Each vehicle reaches target SOC
        for v in range(self.n_vehicles):
            energy_needed = (
                (target_soc[v] - arrival_soc[v]) / 100 * self.battery_capacity
            )

            total_charged = pulp.lpSum([
                charging[v, t]
                for t in range(int(arrival_times[v]), int(departure_times[v]))
            ])

            prob += total_charged >= energy_needed, f"target_soc_v{v}"

        # Constraint 2: Charger capacity (max n_chargers in use at once)
        for t in range(self.time_slots):
            prob += (
                pulp.lpSum([is_charging[v, t] for v in range(self.n_vehicles)]) <=
                self.n_chargers
            ), f"charger_capacity_t{t}"

        # Constraint 3: Link charging amount to is_charging binary
        for v in range(self.n_vehicles):
            for t in range(self.time_slots):
                prob += charging[v, t] <= self.charger_power * is_charging[v, t]

        # Constraint 4: Only charge when vehicle is available
        for v in range(self.n_vehicles):
            for t in range(self.time_slots):
                if t < arrival_times[v] or t >= departure_times[v]:
                    prob += charging[v, t] == 0

        # Solve
        prob.solve(pulp.PULP_CBC_CMD(msg=0))

        # Extract solution
        schedule = np.zeros((self.n_vehicles, self.time_slots))
        for v in range(self.n_vehicles):
            for t in range(self.time_slots):
                schedule[v, t] = pulp.value(charging[v, t])

        total_cost_value = pulp.value(prob.objective)
        total_energy = schedule.sum()

        return {
            'schedule': schedule,
            'total_cost_usd': total_cost_value,
            'total_energy_kwh': total_energy,
            'status': pulp.LpStatus[prob.status]
        }

    def visualize_schedule(self, schedule: np.ndarray, electricity_prices: np.ndarray):
        """
        Visualize charging schedule.

        Args:
            schedule: Charging schedule matrix (n_vehicles, time_slots)
            electricity_prices: Price per time slot
        """
        import matplotlib.pyplot as plt

        fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8), sharex=True)

        # Heatmap of charging schedule
        im = ax1.imshow(schedule, aspect='auto', cmap='YlOrRd')
        ax1.set_ylabel('Vehicle ID')
        ax1.set_title('Charging Schedule (kWh per hour)')
        plt.colorbar(im, ax=ax1, label='Energy (kWh)')

        # Electricity prices
        ax2.bar(range(self.time_slots), electricity_prices, color='steelblue')
        ax2.set_xlabel('Hour of Day')
        ax2.set_ylabel('Price ($/kWh)')
        ax2.set_title('Time-of-Use Electricity Pricing')

        plt.tight_layout()
        plt.show()


# Example usage
if __name__ == "__main__":
    # Fleet parameters
    n_vehicles = 10
    n_chargers = 5

    # Vehicle data
    np.random.seed(42)
    arrival_soc = np.random.uniform(20, 50, n_vehicles)  # Return with 20-50% SOC
    target_soc = np.full(n_vehicles, 90.0)  # Charge to 90%
    arrival_times = np.random.randint(18, 20, n_vehicles)  # Return 6-8 PM
    departure_times = np.random.randint(6, 8, n_vehicles)  # Depart 6-8 AM next day

    # Time-of-use pricing (higher during peak hours)
    electricity_prices = np.array([
        0.08, 0.08, 0.08, 0.08, 0.08, 0.08,  # 12AM-6AM: Off-peak
        0.15, 0.15, 0.20, 0.20, 0.15, 0.15,  # 6AM-12PM: Morning peak
        0.12, 0.12, 0.12, 0.12, 0.20, 0.20,  # 12PM-6PM: Afternoon peak
        0.18, 0.18, 0.12, 0.12, 0.08, 0.08   # 6PM-12AM: Evening/off-peak
    ])

    # Optimize
    optimizer = ChargingStrategyOptimizer(
        n_vehicles=n_vehicles,
        n_chargers=n_chargers,
        charger_power_kw=50.0
    )

    result = optimizer.optimize_schedule(
        arrival_soc=arrival_soc,
        target_soc=target_soc,
        arrival_times=arrival_times,
        departure_times=departure_times,
        electricity_prices=electricity_prices
    )

    print(f"Optimization Status: {result['status']}")
    print(f"Total Cost: ${result['total_cost_usd']:.2f}")
    print(f"Total Energy: {result['total_energy_kwh']:.2f} kWh")

    # Visualize
    optimizer.visualize_schedule(result['schedule'], electricity_prices)

Eco-Routing with Multi-Objective Optimization

python
import numpy as np
import pandas as pd
from typing import List, Tuple
from scipy.optimize import minimize

class EcoRouter:
    """
    Multi-objective route optimization: minimize energy AND time.

    Uses Pareto optimization to explore trade-off frontier.
    """

    def __init__(
        self,
        road_network: pd.DataFrame,
        battery_capacity_kwh: float = 75.0
    ):
        """
        Args:
            road_network: Graph edges with distance, elevation, speed_limit
            battery_capacity_kwh: Vehicle battery capacity
        """
        self.road_network = road_network
        self.battery_capacity = battery_capacity_kwh

    def compute_route_metrics(
        self,
        route: List[int],
        preference_weight: float = 0.5
    ) -> Tuple[float, float, float]:
        """
        Compute energy, time, and weighted score for route.

        Args:
            route: List of node IDs
            preference_weight: Weight for energy vs time (0=all time, 1=all energy)

        Returns:
            energy_kwh, time_hours, weighted_score
        """
        total_energy = 0
        total_time = 0

        for i in range(len(route) - 1):
            src, dst = route[i], route[i + 1]

            # Find edge
            edge = self.road_network[
                (self.road_network['src'] == src) &
                (self.road_network['dst'] == dst)
            ]

            if edge.empty:
                return float('inf'), float('inf'), float('inf')

            edge = edge.iloc[0]

            # Energy (simplified model)
            base_consumption = 18.0  # kWh per 100 km
            elevation_energy = max(0, edge['elevation_gain_m'] / 100)
            segment_energy = (
                base_consumption * edge['distance_km'] / 100 + elevation_energy
            )

            # Time
            segment_time = edge['distance_km'] / edge['speed_limit_kmh']

            total_energy += segment_energy
            total_time += segment_time

        # Weighted score (normalized)
        # Assuming typical trip: 50 km, 0.5 hours, 10 kWh
        energy_normalized = total_energy / 10
        time_normalized = total_time / 0.5

        weighted_score = (
            preference_weight * energy_normalized +
            (1 - preference_weight) * time_normalized
        )

        return total_energy, total_time, weighted_score

    def find_pareto_routes(
        self,
        start: int,
        end: int,
        candidate_routes: List[List[int]],
        n_pareto: int = 5
    ) -> pd.DataFrame:
        """
        Find Pareto-optimal routes (non-dominated solutions).

        Args:
            start: Start node
            end: End node
            candidate_routes: List of possible routes
            n_pareto: Number of Pareto solutions to return

        Returns:
            DataFrame with Pareto-optimal routes
        """
        results = []

        for route in candidate_routes:
            if route[0] != start or route[-1] != end:
                continue

            energy, time, _ = self.compute_route_metrics(route, preference_weight=0.5)

            if energy < float('inf'):
                results.append({
                    'route': route,
                    'energy_kwh': energy,
                    'time_hours': time
                })

        df = pd.DataFrame(results)

        # Find Pareto frontier (non-dominated solutions)
        is_pareto = np.ones(len(df), dtype=bool)

        for i in range(len(df)):
            for j in range(len(df)):
                if i == j:
                    continue

                # j dominates i if j is better in both objectives
                if (df.loc[j, 'energy_kwh'] <= df.loc[i, 'energy_kwh'] and
                    df.loc[j, 'time_hours'] <= df.loc[i, 'time_hours'] and
                    (df.loc[j, 'energy_kwh'] < df.loc[i, 'energy_kwh'] or
                     df.loc[j, 'time_hours'] < df.loc[i, 'time_hours'])):
                    is_pareto[i] = False
                    break

        pareto_df = df[is_pareto].sort_values('energy_kwh').head(n_pareto)

        return pareto_df


# Example usage
if __name__ == "__main__":
    # Sample road network
    road_network = pd.DataFrame({
        'src': [0, 0, 1, 1, 2],
        'dst': [1, 2, 3, 2, 3],
        'distance_km': [30, 40, 20, 15, 25],
        'elevation_gain_m': [100, 50, 200, 300, 100],
        'speed_limit_kmh': [80, 100, 60, 50, 80]
    })

    # Candidate routes (could be generated by routing algorithm)
    candidate_routes = [
        [0, 1, 3],  # Route A
        [0, 2, 3],  # Route B
        [0, 1, 2, 3]  # Route C
    ]

    router = EcoRouter(road_network)

    # Find Pareto-optimal routes
    pareto_routes = router.find_pareto_routes(
        start=0,
        end=3,
        candidate_routes=candidate_routes
    )

    print("Pareto-Optimal Routes:")
    print(pareto_routes[['energy_kwh', 'time_hours']])

Production Deployment

yaml
# Energy optimization pipeline
architecture:
  route_service:
    framework: FastAPI
    endpoints:
      - POST /api/v1/optimize-route (request: origin, destination, preferences)
      - GET /api/v1/charging-stations (params: location, radius)
    latency: <500ms p95

  optimization_engine:
    algorithm: DQN (pre-trained) + heuristic fallback
    model_storage: S3 / Azure Blob
    inference: CPU (sufficient for routing)

  charging_optimizer:
    algorithm: Linear programming (PuLP)
    trigger: Nightly batch (10 PM)
    runtime: <5 minutes for 100 vehicles

  data_sources:
    - Road network: OpenStreetMap + HERE Maps
    - Traffic: Real-time API (Google, TomTom)
    - Weather: OpenWeatherMap API
    - Charging stations: PlugShare, ChargePoint API
    - Electricity prices: Utility API

monitoring:
  metrics:
    - Route quality (energy savings vs baseline)
    - Prediction accuracy (estimated vs actual energy)
    - API latency
    - Model drift (consumption patterns change)

  alerts:
    - Prediction error > 20%
    - API latency > 1s
    - Charging optimizer failure

Production Checklist

• ] Road network data updated (quarterly refresh)
• ] Traffic API integrated with fallback to historical averages
• ] Energy model calibrated on real vehicle data (R² > 0.85)
• ] Pareto routes validated by drivers (user study)
• ] Charging optimizer tested for feasibility (always finds solution)
• ] API rate limits configured (avoid external API overages)
• ] Model versioning in place (A/B test new algorithms)
• ] User preference storage (remember eco vs fast preference)
• ] Privacy compliance (location data anonymization)
• ] Offline mode (cache recent routes, fallback to heuristics)

Fleet Analytics

Fleet Analytics for Connected Vehicles

Build comprehensive analytics dashboards and KPI tracking for vehicle fleets. Focus on operational efficiency, cost optimization, and performance monitoring.

Key Performance Indicators (KPIs)

Vehicle Health

• Average SOH: Fleet-wide battery health
• Maintenance Compliance: % vehicles on schedule
• Fault Rate: Faults per 1000 km
• Downtime: Hours unavailable per vehicle per month

Energy & Efficiency

• Energy Efficiency: kWh per 100 km (fleet average, by vehicle type)
• Charging Efficiency: % of energy delivered vs drawn from grid
• Idle Time: % time parked with systems active
• Regenerative Braking: % energy recovered

Utilization

• Fleet Utilization: % time vehicles in use
• Distance per Day: Average km per vehicle per day
• Trip Count: Number of trips per vehicle per week
• Occupancy Rate: % trips with passengers (if applicable)

Cost Metrics

• Total Cost of Ownership (TCO): Per vehicle per year
• Energy Cost: $ per kWh charged
• Maintenance Cost: $ per km
• Insurance Cost: $ per vehicle per year

Safety & Compliance

• Incident Rate: Incidents per million km
• Driver Safety Score: 0-100 composite score
• Compliance Violations: Speeding, harsh braking events
• Recall Compliance: % vehicles with open recalls addressed

Dashboard Architecture

python
import pandas as pd
import numpy as np
import plotly.graph_objects as go
import plotly.express as px
from plotly.subplots import make_subplots
from typing import Dict, List
import streamlit as st
from datetime import datetime, timedelta

class FleetAnalyticsDashboard:
    """
    Comprehensive fleet analytics dashboard with real-time KPIs.

    Data sources:
    - Vehicle telemetry (PostgreSQL/TimescaleDB)
    - Maintenance records (ERP system)
    - Charging sessions (Charging management platform)
    - Driver profiles (Driver management system)
    """

    def __init__(self, db_connection):
        """
        Args:
            db_connection: Database connection object
        """
        self.db = db_connection

    def compute_fleet_kpis(
        self,
        start_date: datetime,
        end_date: datetime,
        vehicle_filter: List[str] = None
    ) -> Dict:
        """
        Compute all fleet-level KPIs for specified time range.

        Args:
            start_date: Start of analysis period
            end_date: End of analysis period
            vehicle_filter: Optional list of vehicle IDs to analyze

        Returns:
            Dictionary of KPIs
        """
        # Query vehicle data
        query = f"""
        SELECT
            vehicle_id,
            timestamp,
            odometer_km,
            soc,
            soh,
            energy_consumed_kwh,
            fault_codes,
            is_driving,
            driver_id
        FROM vehicle_telemetry
        WHERE timestamp BETWEEN '{start_date}' AND '{end_date}'
        """

        if vehicle_filter:
            query += f" AND vehicle_id IN ({','.join(map(repr, vehicle_filter))})"

        df = pd.read_sql(query, self.db)

        kpis = {}

        # Vehicle Health KPIs
        kpis['avg_soh'] = df.groupby('vehicle_id')['soh'].last().mean()
        kpis['soh_std'] = df.groupby('vehicle_id')['soh'].last().std()

        # Fault rate
        total_km = (
            df.groupby('vehicle_id')['odometer_km'].max() -
            df.groupby('vehicle_id')['odometer_km'].min()
        ).sum()

        fault_count = df['fault_codes'].notna().sum()
        kpis['fault_rate_per_1000km'] = (fault_count / total_km) * 1000 if total_km > 0 else 0

        # Energy efficiency
        total_energy = df['energy_consumed_kwh'].sum()
        kpis['fleet_efficiency_kwh_per_100km'] = (total_energy / total_km) * 100 if total_km > 0 else 0

        # Utilization
        total_hours = (end_date - start_date).total_seconds() / 3600
        driving_hours = df[df['is_driving']]['timestamp'].count() * (1/3600)  # Assuming 1 Hz data
        kpis['fleet_utilization_pct'] = (driving_hours / (total_hours * df['vehicle_id'].nunique())) * 100

        # Distance metrics
        n_vehicles = df['vehicle_id'].nunique()
        n_days = (end_date - start_date).days
        kpis['avg_km_per_vehicle_per_day'] = total_km / (n_vehicles * n_days) if n_days > 0 else 0

        # Cost metrics (assuming cost models)
        kpis['total_energy_cost_usd'] = total_energy * 0.12  # $0.12 per kWh
        kpis['energy_cost_per_km_usd'] = kpis['total_energy_cost_usd'] / total_km if total_km > 0 else 0

        return kpis

    def plot_soh_distribution(self, vehicle_data: pd.DataFrame) -> go.Figure:
        """
        Plot SOH distribution across fleet.

        Args:
            vehicle_data: DataFrame with columns: vehicle_id, soh
        """
        fig = go.Figure()

        # Histogram
        fig.add_trace(go.Histogram(
            x=vehicle_data['soh'],
            nbinsx=30,
            name='SOH Distribution',
            marker_color='steelblue'
        ))

        # Add threshold lines
        fig.add_vline(x=80, line_dash="dash", line_color="red",
                      annotation_text="EOL Threshold (80%)")
        fig.add_vline(x=90, line_dash="dash", line_color="orange",
                      annotation_text="Maintenance Alert (90%)")

        fig.update_layout(
            title='Fleet Battery SOH Distribution',
            xaxis_title='State of Health (%)',
            yaxis_title='Number of Vehicles',
            showlegend=False
        )

        return fig

    def plot_energy_efficiency_by_vehicle_type(
        self,
        vehicle_data: pd.DataFrame
    ) -> go.Figure:
        """
        Box plot of energy efficiency by vehicle type.

        Args:
            vehicle_data: Columns: vehicle_type, efficiency_kwh_per_100km
        """
        fig = px.box(
            vehicle_data,
            x='vehicle_type',
            y='efficiency_kwh_per_100km',
            color='vehicle_type',
            title='Energy Efficiency by Vehicle Type',
            labels={
                'efficiency_kwh_per_100km': 'Energy Efficiency (kWh/100km)',
                'vehicle_type': 'Vehicle Type'
            }
        )

        fig.update_layout(showlegend=False)
        return fig

    def plot_utilization_heatmap(
        self,
        usage_data: pd.DataFrame
    ) -> go.Figure:
        """
        Heatmap of fleet utilization by day and hour.

        Args:
            usage_data: Columns: timestamp, vehicle_id, is_driving
        """
        # Aggregate to hourly utilization
        usage_data['hour'] = usage_data['timestamp'].dt.hour
        usage_data['day_of_week'] = usage_data['timestamp'].dt.day_name()

        utilization = (
            usage_data.groupby(['day_of_week', 'hour'])['is_driving']
            .mean() * 100
        ).reset_index()

        # Pivot for heatmap
        pivot = utilization.pivot(index='day_of_week', columns='hour', values='is_driving')

        # Reorder days
        day_order = ['Monday', 'Tuesday', 'Wednesday', 'Thursday', 'Friday', 'Saturday', 'Sunday']
        pivot = pivot.reindex(day_order)

        fig = go.Figure(data=go.Heatmap(
            z=pivot.values,
            x=pivot.columns,
            y=pivot.index,
            colorscale='YlOrRd',
            colorbar=dict(title='Utilization %')
        ))

        fig.update_layout(
            title='Fleet Utilization by Day and Hour',
            xaxis_title='Hour of Day',
            yaxis_title='Day of Week'
        )

        return fig

    def plot_maintenance_cost_trend(
        self,
        maintenance_data: pd.DataFrame
    ) -> go.Figure:
        """
        Line plot of monthly maintenance costs.

        Args:
            maintenance_data: Columns: date, cost_usd, vehicle_id
        """
        # Aggregate to monthly
        monthly = (
            maintenance_data
            .groupby(pd.Grouper(key='date', freq='M'))['cost_usd']
            .agg(['sum', 'mean', 'count'])
            .reset_index()
        )

        fig = make_subplots(specs=[[{"secondary_y": True}]])

        # Total cost
        fig.add_trace(
            go.Bar(
                x=monthly['date'],
                y=monthly['sum'],
                name='Total Cost',
                marker_color='steelblue'
            ),
            secondary_y=False
        )

        # Average cost per vehicle
        fig.add_trace(
            go.Scatter(
                x=monthly['date'],
                y=monthly['mean'],
                name='Avg Cost per Vehicle',
                mode='lines+markers',
                marker_color='darkorange',
                line=dict(width=3)
            ),
            secondary_y=True
        )

        fig.update_layout(
            title='Monthly Maintenance Costs',
            xaxis_title='Month'
        )

        fig.update_yaxes(title_text='Total Cost ($)', secondary_y=False)
        fig.update_yaxes(title_text='Avg Cost per Vehicle ($)', secondary_y=True)

        return fig

    def plot_driver_safety_scores(
        self,
        driver_data: pd.DataFrame
    ) -> go.Figure:
        """
        Bar chart of driver safety scores with risk categories.

        Args:
            driver_data: Columns: driver_id, safety_score, risk_category
        """
        # Sort by safety score
        driver_data = driver_data.sort_values('safety_score', ascending=False)

        # Color mapping
        color_map = {
            'Low Risk': 'green',
            'Medium Risk': 'orange',
            'High Risk': 'red'
        }

        colors = driver_data['risk_category'].map(color_map)

        fig = go.Figure()

        fig.add_trace(go.Bar(
            x=driver_data['driver_id'],
            y=driver_data['safety_score'],
            marker_color=colors,
            text=driver_data['safety_score'].round(1),
            textposition='outside'
        ))

        fig.add_hline(y=70, line_dash="dash", line_color="red",
                      annotation_text="Retraining Threshold")

        fig.update_layout(
            title='Driver Safety Scores',
            xaxis_title='Driver ID',
            yaxis_title='Safety Score (0-100)',
            showlegend=False
        )

        return fig

    def generate_executive_summary(self, kpis: Dict) -> str:
        """
        Generate text summary of key findings.

        Args:
            kpis: Dictionary of computed KPIs

        Returns:
            Markdown-formatted summary
        """
        summary = f"""
## Fleet Executive Summary

### Vehicle Health
- **Average SOH**: {kpis['avg_soh']:.1f}% ± {kpis['soh_std']:.1f}%
- **Fault Rate**: {kpis['fault_rate_per_1000km']:.2f} faults per 1000 km

### Energy & Efficiency
- **Fleet Efficiency**: {kpis['fleet_efficiency_kwh_per_100km']:.2f} kWh per 100 km
- **Total Energy Cost**: ${kpis['total_energy_cost_usd']:.2f}
- **Cost per km**: ${kpis['energy_cost_per_km_usd']:.4f}

### Utilization
- **Fleet Utilization**: {kpis['fleet_utilization_pct']:.1f}%
- **Avg Distance**: {kpis['avg_km_per_vehicle_per_day']:.1f} km per vehicle per day

### Recommendations
"""

        # Add recommendations based on KPIs
        if kpis['avg_soh'] < 85:
            summary += "- **ACTION REQUIRED**: Fleet average SOH below 85%. Schedule battery assessments.\n"

        if kpis['fleet_utilization_pct'] < 50:
            summary += "- **OPTIMIZATION**: Fleet utilization below 50%. Consider fleet size optimization.\n"

        if kpis['fleet_efficiency_kwh_per_100km'] > 25:
            summary += "- **EFFICIENCY**: Energy efficiency above industry average. Investigate driver training needs.\n"

        return summary


# Streamlit Dashboard Application
def main():
    st.set_page_config(page_title="Fleet Analytics Dashboard", layout="wide")

    st.title("🚗 Fleet Analytics Dashboard")

    # Sidebar filters
    st.sidebar.header("Filters")
    date_range = st.sidebar.date_input(
        "Date Range",
        value=(datetime.now() - timedelta(days=30), datetime.now())
    )

    vehicle_types = st.sidebar.multiselect(
        "Vehicle Types",
        options=['Sedan', 'SUV', 'Van', 'Truck'],
        default=['Sedan', 'SUV']
    )

    # Initialize dashboard
    # db_connection = create_database_connection()  # Implement as needed
    # dashboard = FleetAnalyticsDashboard(db_connection)

    # Placeholder data for demo
    st.header("Key Performance Indicators")

    col1, col2, col3, col4 = st.columns(4)
    with col1:
        st.metric("Avg SOH", "92.3%", delta="-1.2%")
    with col2:
        st.metric("Fleet Efficiency", "18.5 kWh/100km", delta="+0.3")
    with col3:
        st.metric("Utilization", "68%", delta="+5%")
    with col4:
        st.metric("Fault Rate", "2.1/1000km", delta="-0.5")

    # Charts
    st.header("Detailed Analytics")

    tab1, tab2, tab3, tab4 = st.tabs(
        ["Vehicle Health", "Energy & Efficiency", "Utilization", "Costs"]
    )

    with tab1:
        st.subheader("Battery SOH Distribution")
        # Sample data
        soh_data = pd.DataFrame({
            'vehicle_id': [f'V{i:03d}' for i in range(100)],
            'soh': np.random.normal(92, 5, 100).clip(70, 100)
        })
        # fig = dashboard.plot_soh_distribution(soh_data)
        # st.plotly_chart(fig, use_container_width=True)

    with tab2:
        st.subheader("Energy Efficiency by Vehicle Type")
        # Charts would go here

    with tab3:
        st.subheader("Fleet Utilization Heatmap")
        # Charts would go here

    with tab4:
        st.subheader("Maintenance Cost Trends")
        # Charts would go here


if __name__ == "__main__":
    main()

Advanced Analytics: Clustering & Segmentation

python
from sklearn.cluster import KMeans, DBSCAN
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

class FleetSegmentation:
    """
    Segment fleet vehicles into clusters based on usage patterns.

    Use cases:
    - Identify maintenance risk groups
    - Optimize charging schedules by usage pattern
    - Personalize driver training
    """

    def __init__(self, n_clusters: int = 5):
        """
        Args:
            n_clusters: Number of vehicle segments
        """
        self.n_clusters = n_clusters
        self.scaler = StandardScaler()
        self.kmeans = KMeans(n_clusters=n_clusters, random_state=42)
        self.pca = PCA(n_components=2)

    def engineer_clustering_features(
        self,
        vehicle_data: pd.DataFrame
    ) -> pd.DataFrame:
        """
        Create features for clustering.

        Expected columns:
        - vehicle_id, total_km, avg_trip_distance_km, max_speed_kmh,
        - harsh_accel_rate, harsh_brake_rate, idle_time_pct,
        - energy_efficiency_kwh_per_100km, soh, age_days
        """
        features = pd.DataFrame()

        # Usage intensity
        features['total_km'] = vehicle_data['total_km']
        features['avg_trip_distance_km'] = vehicle_data['avg_trip_distance_km']

        # Driving style
        features['max_speed_kmh'] = vehicle_data['max_speed_kmh']
        features['harsh_accel_rate'] = vehicle_data['harsh_accel_rate']
        features['harsh_brake_rate'] = vehicle_data['harsh_brake_rate']
        features['idle_time_pct'] = vehicle_data['idle_time_pct']

        # Efficiency & health
        features['energy_efficiency'] = vehicle_data['energy_efficiency_kwh_per_100km']
        features['soh'] = vehicle_data['soh']

        # Age
        features['age_days'] = vehicle_data['age_days']

        return features

    def fit_predict(self, X: pd.DataFrame) -> np.ndarray:
        """
        Cluster vehicles and return cluster labels.

        Args:
            X: Feature matrix

        Returns:
            Cluster labels (0 to n_clusters-1)
        """
        # Normalize features
        X_scaled = self.scaler.fit_transform(X)

        # Cluster
        labels = self.kmeans.fit_predict(X_scaled)

        return labels

    def visualize_clusters(self, X: pd.DataFrame, labels: np.ndarray):
        """
        Visualize clusters in 2D using PCA.

        Args:
            X: Feature matrix
            labels: Cluster labels
        """
        # PCA for visualization
        X_scaled = self.scaler.transform(X)
        X_pca = self.pca.fit_transform(X_scaled)

        plt.figure(figsize=(10, 6))
        scatter = plt.scatter(
            X_pca[:, 0],
            X_pca[:, 1],
            c=labels,
            cmap='viridis',
            alpha=0.6,
            edgecolors='k'
        )
        plt.colorbar(scatter, label='Cluster')
        plt.xlabel(f'PC1 ({self.pca.explained_variance_ratio_[0]:.1%} variance)')
        plt.ylabel(f'PC2 ({self.pca.explained_variance_ratio_[1]:.1%} variance)')
        plt.title('Vehicle Fleet Segmentation (PCA Projection)')
        plt.tight_layout()
        plt.show()

    def describe_clusters(
        self,
        X: pd.DataFrame,
        labels: np.ndarray
    ) -> pd.DataFrame:
        """
        Describe each cluster with summary statistics.

        Returns:
            DataFrame with cluster profiles
        """
        X_with_labels = X.copy()
        X_with_labels['cluster'] = labels

        cluster_profiles = X_with_labels.groupby('cluster').agg([
            'mean', 'std', 'min', 'max', 'count'
        ])

        return cluster_profiles

    def assign_cluster_names(
        self,
        cluster_profiles: pd.DataFrame
    ) -> Dict[int, str]:
        """
        Assign interpretable names to clusters based on characteristics.

        Example logic:
        - High km + low SOH = "High Usage, Aging"
        - Low km + high SOH = "Light Usage, Healthy"
        - High harsh_accel + low efficiency = "Aggressive Drivers"
        """
        names = {}

        for cluster_id in cluster_profiles.index.get_level_values(0).unique():
            profile = cluster_profiles.loc[cluster_id]

            # Extract key metrics (using 'mean' aggregation)
            total_km = profile[('total_km', 'mean')]
            soh = profile[('soh', 'mean')]
            harsh_rate = profile[('harsh_accel_rate', 'mean')]
            efficiency = profile[('energy_efficiency', 'mean')]

            # Rule-based naming
            if total_km > 50000 and soh < 85:
                names[cluster_id] = "High Usage, Aging"
            elif total_km < 20000 and soh > 95:
                names[cluster_id] = "Light Usage, Healthy"
            elif harsh_rate > 5 and efficiency > 20:
                names[cluster_id] = "Aggressive Drivers"
            elif efficiency < 15 and soh > 90:
                names[cluster_id] = "Efficient, Well-Maintained"
            else:
                names[cluster_id] = f"Cluster {cluster_id}"

        return names


# Example usage
if __name__ == "__main__":
    # Load vehicle data
    vehicle_df = pd.read_csv('fleet_vehicle_profiles.csv')

    # Initialize segmentation
    segmenter = FleetSegmentation(n_clusters=5)

    # Engineer features
    X = segmenter.engineer_clustering_features(vehicle_df)

    # Cluster
    labels = segmenter.fit_predict(X)

    # Visualize
    segmenter.visualize_clusters(X, labels)

    # Describe clusters
    profiles = segmenter.describe_clusters(X, labels)
    print("\nCluster Profiles:")
    print(profiles)

    # Assign names
    names = segmenter.assign_cluster_names(profiles)
    print("\nCluster Names:")
    for cluster_id, name in names.items():
        print(f"  Cluster {cluster_id}: {name}")

Real-Time Streaming Analytics

python
from kafka import KafkaConsumer
from influxdb_client import InfluxDBClient, Point
from influxdb_client.client.write_api import SYNCHRONOUS
import json
from typing import Dict
import logging

logger = logging.getLogger(__name__)

class RealTimeFleetAnalytics:
    """
    Stream processing for real-time fleet analytics.

    Architecture:
    - Kafka: Ingest vehicle telemetry
    - Processing: Compute rolling KPIs
    - InfluxDB: Store time-series metrics
    - Grafana: Visualization
    """

    def __init__(
        self,
        kafka_brokers: list,
        kafka_topic: str,
        influx_url: str,
        influx_token: str,
        influx_org: str,
        influx_bucket: str
    ):
        """
        Args:
            kafka_brokers: List of Kafka broker addresses
            kafka_topic: Topic name for vehicle telemetry
            influx_url: InfluxDB server URL
            influx_token: InfluxDB authentication token
            influx_org: InfluxDB organization
            influx_bucket: InfluxDB bucket name
        """
        # Kafka consumer
        self.consumer = KafkaConsumer(
            kafka_topic,
            bootstrap_servers=kafka_brokers,
            value_deserializer=lambda m: json.loads(m.decode('utf-8')),
            auto_offset_reset='latest',
            enable_auto_commit=True
        )

        # InfluxDB client
        self.influx_client = InfluxDBClient(
            url=influx_url,
            token=influx_token,
            org=influx_org
        )
        self.write_api = self.influx_client.write_api(write_options=SYNCHRONOUS)
        self.bucket = influx_bucket
        self.org = influx_org

        # State management (for rolling metrics)
        self.vehicle_states = {}

    def process_telemetry_message(self, message: Dict):
        """
        Process single telemetry message and compute metrics.

        Expected message format:
        {
            "vehicle_id": "V001",
            "timestamp": "2024-03-19T10:00:00Z",
            "odometer_km": 12345.6,
            "soc": 75,
            "soh": 92,
            "voltage": 400.5,
            "current": 150.2,
            "temperature": 35.0,
            "speed_kmh": 80,
            "is_driving": true
        }
        """
        vehicle_id = message['vehicle_id']
        timestamp = message['timestamp']

        # Update vehicle state
        if vehicle_id not in self.vehicle_states:
            self.vehicle_states[vehicle_id] = {
                'last_odometer': message['odometer_km'],
                'energy_consumed_kwh': 0,
                'trip_count': 0
            }

        state = self.vehicle_states[vehicle_id]

        # Compute derived metrics
        distance_delta = message['odometer_km'] - state['last_odometer']

        # Power consumption (simplified)
        power_kw = (message['voltage'] * message['current']) / 1000
        energy_delta_kwh = power_kw * (1 / 3600)  # Assuming 1 second interval

        state['energy_consumed_kwh'] += energy_delta_kwh
        state['last_odometer'] = message['odometer_km']

        # Trip detection
        if message['is_driving'] and not state.get('was_driving', False):
            state['trip_count'] += 1

        state['was_driving'] = message['is_driving']

        # Write metrics to InfluxDB
        self.write_vehicle_metrics(message, distance_delta, energy_delta_kwh)

    def write_vehicle_metrics(
        self,
        message: Dict,
        distance_delta: float,
        energy_delta: float
    ):
        """
        Write processed metrics to InfluxDB.
        """
        vehicle_id = message['vehicle_id']
        timestamp = message['timestamp']

        # Core telemetry
        point = Point("vehicle_telemetry") \
            .tag("vehicle_id", vehicle_id) \
            .field("soc", message['soc']) \
            .field("soh", message['soh']) \
            .field("voltage", message['voltage']) \
            .field("current", message['current']) \
            .field("temperature", message['temperature']) \
            .field("speed_kmh", message['speed_kmh']) \
            .time(timestamp)

        self.write_api.write(bucket=self.bucket, org=self.org, record=point)

        # Derived metrics
        if distance_delta > 0:
            efficiency_kwh_per_100km = (energy_delta / distance_delta) * 100

            point_derived = Point("vehicle_metrics") \
                .tag("vehicle_id", vehicle_id) \
                .field("energy_consumed_kwh", energy_delta) \
                .field("distance_km", distance_delta) \
                .field("efficiency_kwh_per_100km", efficiency_kwh_per_100km) \
                .time(timestamp)

            self.write_api.write(bucket=self.bucket, org=self.org, record=point_derived)

    def run(self):
        """
        Start consuming messages from Kafka and process in real-time.
        """
        logger.info("Starting real-time fleet analytics processor...")

        try:
            for message in self.consumer:
                try:
                    self.process_telemetry_message(message.value)
                except Exception as e:
                    logger.error(f"Error processing message: {e}")

        except KeyboardInterrupt:
            logger.info("Shutting down...")
        finally:
            self.consumer.close()
            self.influx_client.close()


# Example usage
if __name__ == "__main__":
    processor = RealTimeFleetAnalytics(
        kafka_brokers=['localhost:9092'],
        kafka_topic='vehicle-telemetry',
        influx_url='http://localhost:8086',
        influx_token='your-influx-token',
        influx_org='your-org',
        influx_bucket='fleet-analytics'
    )

    processor.run()

Deployment Checklist

• ] Data pipeline validated (Kafka -> Processing -> Storage)
• ] Dashboard responsive (<2s load time)
• ] Real-time metrics update frequency configured (1-10s)
• ] Historical data retention policy defined
• ] Alert thresholds configured for critical KPIs
• ] User access controls implemented (RBAC)
• ] Dashboard mobile-responsive
• ] Export functionality (CSV, PDF reports)
• ] Scheduled reports automated (daily, weekly, monthly)
• ] Data quality monitoring (missing data, outliers)

Predictive Maintenance

Predictive Maintenance for Automotive Components

Predict component failures before they occur using supervised ML and survival analysis. Focus on battery SOH, tire wear, brake pad life, and electric motor degradation.

Use Cases

1. Battery State of Health (SOH): Predict remaining useful life, capacity fade, impedance rise
2. Tire Wear: Estimate tread depth, predict replacement needs
3. Brake System: Brake pad thickness, rotor wear, fluid degradation
4. Electric Motor: Bearing wear, insulation degradation, magnet demagnetization
5. Cooling System: Pump failure, radiator clogging, coolant degradation

Problem Formulation

Regression (Remaining Useful Life)

Predict time-to-failure or remaining capacity as continuous value.

Algorithms:

• Gradient Boosting (XGBoost, LightGBM, CatBoost)
• Random Forest
• LSTM/GRU for sequential degradation

Classification (Failure within Window)

Predict if failure will occur in next N days/cycles.

Algorithms:

• Gradient Boosting Classifiers
• Logistic Regression (baseline)
• Neural Networks

Survival Analysis (Time-to-Event)

Model probability of survival beyond time t, handling censored data.

Algorithms:

• Cox Proportional Hazards
• Random Survival Forests
• Deep survival models (DeepSurv)

Battery SOH Prediction

Feature Engineering

python
import numpy as np
import pandas as pd
from scipy.stats import skew, kurtosis
from typing import Dict

class BatterySOHFeatureEngineer:
    """
    Extract features from battery charge/discharge cycles for SOH prediction.

    Features capture:
    - Capacity fade
    - Impedance rise
    - Charge acceptance degradation
    - Thermal behavior changes
    """

    @staticmethod
    def extract_cycle_features(
        cycle_data: pd.DataFrame,
        metadata: Dict
    ) -> Dict[str, float]:
        """
        Extract features from a single charge/discharge cycle.

        Expected columns:
        - timestamp, voltage, current, temperature, soc
        """
        features = {}

        # Metadata
        features['cycle_number'] = metadata['cycle_number']
        features['total_kwh_throughput'] = metadata['total_kwh_throughput']
        features['age_days'] = metadata['age_days']
        features['ambient_temp_avg'] = metadata.get('ambient_temp_avg', 25.0)

        # Capacity features
        charge_data = cycle_data[cycle_data['current'] > 0]
        discharge_data = cycle_data[cycle_data['current'] < 0]

        if len(charge_data) > 0:
            features['charge_capacity_ah'] = (
                charge_data['current'].sum() *
                (charge_data['timestamp'].diff().dt.total_seconds().mean() / 3600)
            )
            features['charge_duration_min'] = (
                (charge_data['timestamp'].max() - charge_data['timestamp'].min())
                .total_seconds() / 60
            )
        else:
            features['charge_capacity_ah'] = 0
            features['charge_duration_min'] = 0

        if len(discharge_data) > 0:
            features['discharge_capacity_ah'] = abs(
                discharge_data['current'].sum() *
                (discharge_data['timestamp'].diff().dt.total_seconds().mean() / 3600)
            )
            features['discharge_duration_min'] = (
                (discharge_data['timestamp'].max() - discharge_data['timestamp'].min())
                .total_seconds() / 60
            )
        else:
            features['discharge_capacity_ah'] = 0
            features['discharge_duration_min'] = 0

        # Coulombic efficiency
        if features['charge_capacity_ah'] > 0:
            features['coulombic_efficiency'] = (
                features['discharge_capacity_ah'] / features['charge_capacity_ah']
            )
        else:
            features['coulombic_efficiency'] = 1.0

        # Voltage features
        features['voltage_mean'] = cycle_data['voltage'].mean()
        features['voltage_std'] = cycle_data['voltage'].std()
        features['voltage_min'] = cycle_data['voltage'].min()
        features['voltage_max'] = cycle_data['voltage'].max()

        # dV/dSOC (charge acceptance indicator)
        if len(charge_data) > 10:
            soc_bins = pd.cut(charge_data['soc'], bins=10)
            voltage_per_soc = charge_data.groupby(soc_bins)['voltage'].mean()
            features['dv_dsoc_slope'] = np.polyfit(
                range(len(voltage_per_soc)),
                voltage_per_soc.values,
                1
            )[0]
        else:
            features['dv_dsoc_slope'] = 0

        # Impedance proxy (voltage drop at constant current)
        if len(discharge_data) > 100:
            # Find steady-state discharge region (SOC 80% -> 20%)
            steady_discharge = discharge_data[
                (discharge_data['soc'] >= 20) & (discharge_data['soc'] <= 80)
            ]
            if len(steady_discharge) > 10:
                # Impedance ~ voltage_drop / current
                current_mean = steady_discharge['current'].mean()
                voltage_drop_per_amp = steady_discharge['voltage'].std() / abs(current_mean)
                features['impedance_proxy'] = voltage_drop_per_amp
            else:
                features['impedance_proxy'] = 0
        else:
            features['impedance_proxy'] = 0

        # Thermal features
        features['temp_mean'] = cycle_data['temperature'].mean()
        features['temp_max'] = cycle_data['temperature'].max()
        features['temp_std'] = cycle_data['temperature'].std()
        features['temp_rise'] = (
            cycle_data['temperature'].max() - cycle_data['temperature'].min()
        )

        # Temperature-power correlation (thermal management quality)
        if len(cycle_data) > 10:
            power = cycle_data['voltage'] * cycle_data['current']
            temp_power_corr = cycle_data['temperature'].corr(power)
            features['temp_power_correlation'] = temp_power_corr
        else:
            features['temp_power_correlation'] = 0

        # Resting voltage (open-circuit voltage after rest period)
        # Indicates equilibrium state
        rest_data = cycle_data[abs(cycle_data['current']) < 1.0]  # <1A = rest
        if len(rest_data) > 5:
            features['rest_voltage'] = rest_data['voltage'].mean()
        else:
            features['rest_voltage'] = cycle_data['voltage'].mean()

        # Statistical features
        features['voltage_skewness'] = skew(cycle_data['voltage'])
        features['voltage_kurtosis'] = kurtosis(cycle_data['voltage'])

        return features

    @staticmethod
    def create_degradation_trends(
        cycle_features_df: pd.DataFrame,
        window: int = 10
    ) -> pd.DataFrame:
        """
        Create trend features from cycle history.

        Args:
            cycle_features_df: Features per cycle (sorted by cycle_number)
            window: Rolling window size for trend calculation
        """
        trends = cycle_features_df.copy()

        # Capacity fade rate
        trends['capacity_fade_rate'] = (
            -trends['discharge_capacity_ah']
            .rolling(window=window, min_periods=3)
            .apply(lambda x: np.polyfit(range(len(x)), x, 1)[0])
        )

        # Impedance rise rate
        trends['impedance_rise_rate'] = (
            trends['impedance_proxy']
            .rolling(window=window, min_periods=3)
            .apply(lambda x: np.polyfit(range(len(x)), x, 1)[0])
        )

        # Efficiency degradation
        trends['efficiency_decline_rate'] = (
            -trends['coulombic_efficiency']
            .rolling(window=window, min_periods=3)
            .apply(lambda x: np.polyfit(range(len(x)), x, 1)[0])
        )

        # Temperature increase trend (thermal management degradation)
        trends['temp_increase_rate'] = (
            trends['temp_mean']
            .rolling(window=window, min_periods=3)
            .apply(lambda x: np.polyfit(range(len(x)), x, 1)[0])
        )

        # Volatility features (increased variance = instability)
        trends['capacity_volatility'] = (
            trends['discharge_capacity_ah'].rolling(window=window).std()
        )

        return trends.fillna(0)


# Example usage
if __name__ == "__main__":
    # Load cycle data
    cycle_df = pd.read_parquet('battery_cycle_001.parquet')
    metadata = {
        'cycle_number': 1,
        'total_kwh_throughput': 120.5,
        'age_days': 45
    }

    engineer = BatterySOHFeatureEngineer()
    features = engineer.extract_cycle_features(cycle_df, metadata)

    print("Extracted features:")
    for key, value in features.items():
        print(f"  {key}: {value:.4f}")

Model Implementation: Gradient Boosting

python
import lightgbm as lgb
import numpy as np
import pandas as pd
from sklearn.model_selection import TimeSeriesSplit
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
from typing import Tuple
import joblib
import logging

logger = logging.getLogger(__name__)

class BatterySOHPredictor:
    """
    LightGBM-based State of Health predictor.

    Predicts:
    - Current SOH (0-100%)
    - Remaining useful cycles (until 80% SOH)
    - Capacity at next cycle
    """

    def __init__(self, model_params: dict = None):
        """
        Args:
            model_params: LightGBM hyperparameters
        """
        default_params = {
            'objective': 'regression',
            'metric': 'rmse',
            'boosting_type': 'gbdt',
            'num_leaves': 31,
            'learning_rate': 0.05,
            'feature_fraction': 0.8,
            'bagging_fraction': 0.8,
            'bagging_freq': 5,
            'max_depth': 6,
            'min_child_samples': 20,
            'verbosity': -1,
            'n_jobs': -1
        }
        self.params = model_params or default_params
        self.model = None
        self.feature_names = None

    def train(
        self,
        X: pd.DataFrame,
        y: np.ndarray,
        n_splits: int = 5,
        n_estimators: int = 500,
        early_stopping_rounds: int = 50
    ) -> dict:
        """
        Train SOH predictor with time-series cross-validation.

        Args:
            X: Feature matrix (cycle features)
            y: Target (SOH percentage, 0-100)
            n_splits: Number of CV splits
            n_estimators: Number of boosting rounds
            early_stopping_rounds: Early stopping patience

        Returns:
            Training metrics
        """
        self.feature_names = X.columns.tolist()

        # Time-series cross-validation (preserve temporal order)
        tscv = TimeSeriesSplit(n_splits=n_splits)

        cv_scores = []
        for fold, (train_idx, val_idx) in enumerate(tscv.split(X)):
            X_train, X_val = X.iloc[train_idx], X.iloc[val_idx]
            y_train, y_val = y[train_idx], y[val_idx]

            train_data = lgb.Dataset(X_train, label=y_train)
            val_data = lgb.Dataset(X_val, label=y_val, reference=train_data)

            model = lgb.train(
                self.params,
                train_data,
                num_boost_round=n_estimators,
                valid_sets=[val_data],
                valid_names=['val'],
                callbacks=[
                    lgb.early_stopping(early_stopping_rounds),
                    lgb.log_evaluation(period=0)  # Silent
                ]
            )

            # Evaluate
            y_pred = model.predict(X_val)
            mae = mean_absolute_error(y_val, y_pred)
            rmse = np.sqrt(mean_squared_error(y_val, y_pred))
            r2 = r2_score(y_val, y_pred)

            cv_scores.append({'mae': mae, 'rmse': rmse, 'r2': r2})
            logger.info(f"Fold {fold+1}: MAE={mae:.3f}, RMSE={rmse:.3f}, R2={r2:.3f}")

        # Train final model on full data
        full_data = lgb.Dataset(X, label=y)
        self.model = lgb.train(
            self.params,
            full_data,
            num_boost_round=n_estimators
        )

        # Aggregate CV metrics
        avg_metrics = {
            'mae': np.mean([s['mae'] for s in cv_scores]),
            'rmse': np.mean([s['rmse'] for s in cv_scores]),
            'r2': np.mean([s['r2'] for s in cv_scores]),
            'cv_scores': cv_scores
        }

        logger.info(f"Average CV MAE: {avg_metrics['mae']:.3f}%")
        logger.info(f"Average CV RMSE: {avg_metrics['rmse']:.3f}%")
        logger.info(f"Average CV R2: {avg_metrics['r2']:.3f}")

        return avg_metrics

    def predict(self, X: pd.DataFrame) -> np.ndarray:
        """
        Predict SOH for new samples.

        Returns:
            Predicted SOH values (0-100%)
        """
        if self.model is None:
            raise ValueError("Model must be trained before prediction")

        return self.model.predict(X)

    def predict_with_uncertainty(
        self,
        X: pd.DataFrame,
        n_iterations: int = 100
    ) -> Tuple[np.ndarray, np.ndarray]:
        """
        Predict with uncertainty estimation using dropout simulation.

        Returns:
            mean_predictions: Mean predicted SOH
            std_predictions: Standard deviation (uncertainty)
        """
        # LightGBM doesn't have native dropout, so we use bootstrapping
        predictions = []

        for _ in range(n_iterations):
            # Random feature subsample
            feature_fraction = self.params.get('feature_fraction', 0.8)
            n_features = int(len(X.columns) * feature_fraction)
            random_features = np.random.choice(
                X.columns,
                size=n_features,
                replace=False
            )

            pred = self.model.predict(X[random_features])
            predictions.append(pred)

        predictions = np.array(predictions)
        mean_pred = predictions.mean(axis=0)
        std_pred = predictions.std(axis=0)

        return mean_pred, std_pred

    def predict_remaining_cycles(
        self,
        X: pd.DataFrame,
        current_soh: np.ndarray,
        eol_threshold: float = 80.0
    ) -> np.ndarray:
        """
        Predict remaining useful cycles until end-of-life.

        Args:
            X: Current feature matrix
            current_soh: Current SOH values
            eol_threshold: End-of-life SOH threshold (default 80%)

        Returns:
            Estimated remaining cycles
        """
        # Estimate degradation rate from recent trend
        if 'capacity_fade_rate' in X.columns:
            fade_rates = X['capacity_fade_rate'].values
            # Prevent division by zero
            fade_rates = np.where(np.abs(fade_rates) < 1e-6, -0.1, fade_rates)

            # RUL = (current_soh - eol_threshold) / abs(fade_rate)
            rul = (current_soh - eol_threshold) / np.abs(fade_rates)

            # Clip to reasonable range [0, 5000]
            rul = np.clip(rul, 0, 5000)
        else:
            # Fallback: assume constant 0.05% fade per cycle
            rul = (current_soh - eol_threshold) / 0.05

        return rul

    def feature_importance(self, plot: bool = False) -> pd.DataFrame:
        """
        Get feature importance scores.

        Args:
            plot: If True, display bar plot

        Returns:
            DataFrame with feature names and importance scores
        """
        if self.model is None:
            raise ValueError("Model must be trained before feature importance")

        importance = self.model.feature_importance(importance_type='gain')
        importance_df = pd.DataFrame({
            'feature': self.feature_names,
            'importance': importance
        }).sort_values('importance', ascending=False)

        if plot:
            import matplotlib.pyplot as plt
            plt.figure(figsize=(10, 6))
            plt.barh(importance_df['feature'][:15], importance_df['importance'][:15])
            plt.xlabel('Importance (Gain)')
            plt.title('Top 15 Features for SOH Prediction')
            plt.gca().invert_yaxis()
            plt.tight_layout()
            plt.show()

        return importance_df

    def save(self, path: str):
        """Save trained model."""
        if self.model is None:
            raise ValueError("Cannot save untrained model")

        joblib.dump({
            'model': self.model,
            'params': self.params,
            'feature_names': self.feature_names
        }, path)
        logger.info(f"Model saved to {path}")

    @classmethod
    def load(cls, path: str) -> 'BatterySOHPredictor':
        """Load trained model."""
        data = joblib.load(path)
        predictor = cls(model_params=data['params'])
        predictor.model = data['model']
        predictor.feature_names = data['feature_names']
        logger.info(f"Model loaded from {path}")
        return predictor


# Example usage
if __name__ == "__main__":
    # Load preprocessed cycle features
    df = pd.read_parquet('battery_cycle_features.parquet')

    # Features and target
    feature_cols = [col for col in df.columns if col not in ['soh', 'battery_id']]
    X = df[feature_cols]
    y = df['soh'].values

    # Train predictor
    predictor = BatterySOHPredictor()
    metrics = predictor.train(X, y, n_splits=5)

    # Feature importance
    importance_df = predictor.feature_importance()
    print("\nTop 10 Most Important Features:")
    print(importance_df.head(10))

    # Save model
    predictor.save('battery_soh_predictor.pkl')

    # Predict with uncertainty
    X_test = X.iloc[-100:]
    mean_pred, std_pred = predictor.predict_with_uncertainty(X_test)

    print("\nPredictions with uncertainty:")
    for i in range(5):
        print(f"  Sample {i}: {mean_pred[i]:.2f}% +/- {std_pred[i]:.2f}%")

    # Remaining useful life
    rul = predictor.predict_remaining_cycles(X_test, mean_pred)
    print(f"\nEstimated RUL: {rul[:5]} cycles")

LSTM for Sequential Degradation

python
import torch
import torch.nn as nn
from torch.utils.data import Dataset, DataLoader
import numpy as np
import pandas as pd

class BatterySequenceDataset(Dataset):
    """Dataset for multi-cycle battery sequences."""

    def __init__(
        self,
        sequences: np.ndarray,
        targets: np.ndarray,
        seq_len: int = 50
    ):
        """
        Args:
            sequences: (N, max_cycles, n_features) array
            targets: (N,) array of target SOH
            seq_len: Sequence length to use
        """
        self.sequences = torch.FloatTensor(sequences[:, -seq_len:, :])
        self.targets = torch.FloatTensor(targets)

    def __len__(self):
        return len(self.sequences)

    def __getitem__(self, idx):
        return self.sequences[idx], self.targets[idx]


class BatterySOHLSTM(nn.Module):
    """
    LSTM-based SOH predictor for sequential degradation modeling.

    Input: (batch, seq_len, n_features) - last N cycles
    Output: (batch,) - predicted SOH at next cycle
    """

    def __init__(
        self,
        n_features: int,
        hidden_dim: int = 64,
        n_layers: int = 2,
        dropout: float = 0.2
    ):
        super().__init__()

        self.lstm = nn.LSTM(
            input_size=n_features,
            hidden_size=hidden_dim,
            num_layers=n_layers,
            batch_first=True,
            dropout=dropout if n_layers > 1 else 0
        )

        self.fc = nn.Sequential(
            nn.Linear(hidden_dim, 32),
            nn.ReLU(),
            nn.Dropout(dropout),
            nn.Linear(32, 1)
        )

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        # x: (batch, seq_len, n_features)
        lstm_out, (hidden, _) = self.lstm(x)

        # Use last hidden state
        out = self.fc(hidden[-1])  # (batch, 1)
        return out.squeeze()  # (batch,)


def train_lstm_predictor(
    train_sequences: np.ndarray,
    train_targets: np.ndarray,
    val_sequences: np.ndarray,
    val_targets: np.ndarray,
    n_features: int,
    epochs: int = 100,
    batch_size: int = 32,
    lr: float = 0.001
) -> BatterySOHLSTM:
    """
    Train LSTM SOH predictor.

    Args:
        train_sequences: (N_train, seq_len, n_features)
        train_targets: (N_train,) - SOH values
        val_sequences: (N_val, seq_len, n_features)
        val_targets: (N_val,)
        n_features: Number of input features
        epochs: Training epochs
        batch_size: Batch size
        lr: Learning rate

    Returns:
        Trained model
    """
    device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')

    # Create datasets
    train_dataset = BatterySequenceDataset(train_sequences, train_targets)
    val_dataset = BatterySequenceDataset(val_sequences, val_targets)

    train_loader = DataLoader(train_dataset, batch_size=batch_size, shuffle=True)
    val_loader = DataLoader(val_dataset, batch_size=batch_size)

    # Model
    model = BatterySOHLSTM(n_features=n_features).to(device)
    optimizer = torch.optim.Adam(model.parameters(), lr=lr)
    criterion = nn.MSELoss()

    best_val_loss = float('inf')
    patience = 10
    patience_counter = 0

    for epoch in range(epochs):
        # Training
        model.train()
        train_loss = 0
        for sequences, targets in train_loader:
            sequences = sequences.to(device)
            targets = targets.to(device)

            optimizer.zero_grad()
            outputs = model(sequences)
            loss = criterion(outputs, targets)
            loss.backward()
            optimizer.step()

            train_loss += loss.item()

        train_loss /= len(train_loader)

        # Validation
        model.eval()
        val_loss = 0
        with torch.no_grad():
            for sequences, targets in val_loader:
                sequences = sequences.to(device)
                targets = targets.to(device)
                outputs = model(sequences)
                loss = criterion(outputs, targets)
                val_loss += loss.item()

        val_loss /= len(val_loader)

        if (epoch + 1) % 10 == 0:
            print(f"Epoch {epoch+1}/{epochs} - "
                  f"Train Loss: {train_loss:.4f}, Val Loss: {val_loss:.4f}")

        # Early stopping
        if val_loss < best_val_loss:
            best_val_loss = val_loss
            patience_counter = 0
            best_model_state = model.state_dict()
        else:
            patience_counter += 1
            if patience_counter >= patience:
                print(f"Early stopping at epoch {epoch+1}")
                break

    # Load best model
    model.load_state_dict(best_model_state)
    return model

Tire Wear Prediction

python
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split
import numpy as np
import pandas as pd

class TireWearPredictor:
    """
    Predict tire tread depth based on usage patterns.

    Features:
    - Mileage
    - Driving style (acceleration, braking intensity)
    - Road conditions (highway vs city, temperature)
    - Tire pressure history
    - Load distribution
    """

    def __init__(self):
        self.model = RandomForestRegressor(
            n_estimators=200,
            max_depth=15,
            min_samples_split=10,
            random_state=42,
            n_jobs=-1
        )
        self.initial_tread_depth_mm = 8.0  # New tire

    def extract_features(self, vehicle_data: pd.DataFrame) -> pd.DataFrame:
        """
        Extract tire wear features from vehicle telemetry.

        Expected columns:
        - mileage, speed, acceleration, brake_pressure,
        - tire_pressure_fl, tire_pressure_fr, tire_pressure_rl, tire_pressure_rr,
        - ambient_temp, road_surface
        """
        features = pd.DataFrame()

        # Usage features
        features['total_mileage_km'] = vehicle_data['mileage']
        features['avg_speed_kmh'] = vehicle_data.groupby('vehicle_id')['speed'].transform('mean')

        # Driving style
        features['harsh_accel_count'] = (
            (vehicle_data['acceleration'] > 2.5)
            .groupby(vehicle_data['vehicle_id']).transform('sum')
        )
        features['harsh_brake_count'] = (
            (vehicle_data['brake_pressure'] > 0.7)
            .groupby(vehicle_data['vehicle_id']).transform('sum')
        )

        # Tire pressure (underinflation accelerates wear)
        tire_cols = ['tire_pressure_fl', 'tire_pressure_fr',
                     'tire_pressure_rl', 'tire_pressure_rr']
        features['avg_tire_pressure_bar'] = vehicle_data[tire_cols].mean(axis=1)
        features['tire_pressure_variance'] = vehicle_data[tire_cols].var(axis=1)

        # Underinflation events
        optimal_pressure = 2.5  # bar
        features['underinflation_events'] = (
            (vehicle_data[tire_cols] < optimal_pressure * 0.9).sum(axis=1)
            .groupby(vehicle_data['vehicle_id']).transform('sum')
        )

        # Environmental
        features['avg_temp_c'] = vehicle_data.groupby('vehicle_id')['ambient_temp'].transform('mean')

        # Road surface distribution
        if 'road_surface' in vehicle_data.columns:
            road_dummies = pd.get_dummies(vehicle_data['road_surface'], prefix='road')
            for col in road_dummies.columns:
                features[col] = road_dummies[col].groupby(vehicle_data['vehicle_id']).transform('mean')

        return features

    def train(self, X: pd.DataFrame, y: np.ndarray):
        """
        Train tire wear predictor.

        Args:
            X: Feature matrix
            y: Measured tread depth [mm]
        """
        X_train, X_val, y_train, y_val = train_test_split(
            X, y, test_size=0.2, random_state=42
        )

        self.model.fit(X_train, y_train)

        # Evaluate
        train_score = self.model.score(X_train, y_train)
        val_score = self.model.score(X_val, y_val)
        print(f"Training R2: {train_score:.3f}")
        print(f"Validation R2: {val_score:.3f}")

    def predict_tread_depth(self, X: pd.DataFrame) -> np.ndarray:
        """Predict current tread depth [mm]."""
        return self.model.predict(X)

    def predict_replacement_mileage(
        self,
        X: pd.DataFrame,
        current_mileage: np.ndarray,
        min_tread_depth_mm: float = 1.6
    ) -> np.ndarray:
        """
        Predict mileage at which tire should be replaced.

        Args:
            X: Current feature matrix
            current_mileage: Current vehicle mileage [km]
            min_tread_depth_mm: Legal minimum tread depth

        Returns:
            Estimated mileage at replacement [km]
        """
        current_tread = self.predict_tread_depth(X)

        # Estimate wear rate (mm per 1000 km)
        wear_so_far = self.initial_tread_depth_mm - current_tread
        wear_rate = wear_so_far / (current_mileage / 1000)

        # Remaining mileage
        remaining_tread = current_tread - min_tread_depth_mm
        remaining_mileage = (remaining_tread / wear_rate) * 1000

        replacement_mileage = current_mileage + remaining_mileage

        return replacement_mileage

Deployment Architecture

yaml
# Predictive maintenance pipeline
pipeline:
  data_ingestion:
    source: Vehicle telemetry stream (Kafka)
    frequency: Real-time for critical, daily batch for non-critical

  feature_engineering:
    engine: Apache Spark / Pandas
    caching: Redis (recent features)

  inference:
    model_serving: TensorFlow Serving / FastAPI
    latency: <100ms for real-time, <1h for batch

  output:
    storage: PostgreSQL (predictions) + Grafana dashboards
    alerts: PagerDuty / Slack for critical predictions
    recommendations: Maintenance scheduling system

monitoring:
  metrics:
    - Prediction accuracy (MAE, RMSE)
    - Model drift (feature distribution shift)
    - Alert precision/recall
    - Business KPIs (maintenance cost reduction, downtime)

  retraining:
    trigger: Performance degradation OR new failure modes
    frequency: Monthly for battery SOH, quarterly for tires
    validation: Holdout test set + field validation

Production Checklist

• ] Training data spans full degradation lifecycle (0 -> EOL)
• ] Model validated on multiple battery/component types
• ] Feature engineering handles missing sensors gracefully
• ] Prediction confidence intervals calibrated
• ] Alert thresholds set based on cost-benefit analysis
• ] Integration with maintenance scheduling system tested
• ] Model versioning and rollback procedure in place
• ] Monitoring dashboards for prediction accuracy
• ] Retraining pipeline automated with data quality checks
• ] A/B testing framework for model updates

Time Series Forecasting

Time-Series Forecasting for Automotive Systems

Forecast battery degradation, energy consumption, and charging demand using specialized time-series models.

Use Cases

1. Battery Degradation Forecasting: Predict SOH trajectory over 5-10 years
2. Energy Consumption Prediction: Forecast energy usage for route planning
3. Charging Demand Forecasting: Predict station load for grid management
4. Range Estimation: Predict remaining range under varying conditions
5. Thermal Management: Forecast cooling/heating needs

Algorithm Selection

Prophet (Meta)

Best for: Business forecasting with seasonality and holidays

Pros:

• Handles missing data and outliers
• Automatic seasonality detection
• Interpretable components (trend, weekly, yearly)
• Fast training

Cons:

• Assumes additive model structure
• Limited for multivariate forecasting

Use cases: Charging station demand, fleet energy consumption

LSTM/GRU

Best for: Complex non-linear temporal dependencies

Pros:

• Captures long-term dependencies
• Handles multivariate inputs
• Flexible architecture

Cons:

• Requires large training data
• Slow to train
• Hyperparameter sensitive

Use cases: Battery SOH trajectory, multi-sensor degradation

ARIMA/SARIMAX

Best for: Stationary time-series with clear seasonality

Pros:

• Statistical rigor
• Confidence intervals
• Interpretable parameters

Cons:

• Requires stationarity
• Struggles with non-linear patterns
• Manual order selection

Use cases: Short-term energy consumption, charging demand

Temporal Fusion Transformer (TFT)

Best for: Multi-horizon forecasting with multiple covariates

Pros:

• State-of-the-art performance
• Attention mechanism for interpretability
• Handles static and dynamic features

Cons:

• Computationally expensive
• Requires significant data
• Complex implementation

Use cases: Long-term battery degradation, fleet-wide energy optimization

Battery SOH Forecasting with Prophet

python
from prophet import Prophet
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from typing import Tuple, Dict

class BatterySOHForecaster:
    """
    Forecast battery State of Health using Facebook Prophet.

    Models degradation as:
    - Trend: Long-term capacity fade
    - Seasonal: Usage patterns (daily, weekly)
    - Regressor: Temperature, charge cycles
    """

    def __init__(
        self,
        growth: str = 'linear',
        seasonality_mode: str = 'additive'
    ):
        """
        Args:
            growth: 'linear' or 'logistic' (for bounded degradation)
            seasonality_mode: 'additive' or 'multiplicative'
        """
        self.model = Prophet(
            growth=growth,
            seasonality_mode=seasonality_mode,
            yearly_seasonality=False,  # Not relevant for batteries
            weekly_seasonality=True,   # Usage patterns
            daily_seasonality=True,    # Charge cycles
            changepoint_prior_scale=0.05  # Flexibility of trend changes
        )
        self.trained = False

    def prepare_data(
        self,
        soh_history: pd.DataFrame,
        temperature_history: pd.DataFrame = None
    ) -> pd.DataFrame:
        """
        Prepare data for Prophet.

        Args:
            soh_history: Columns: timestamp, soh
            temperature_history: Columns: timestamp, avg_temp_c

        Returns:
            DataFrame with columns: ds (timestamp), y (SOH), temp (optional)
        """
        df = soh_history.rename(columns={'timestamp': 'ds', 'soh': 'y'})

        # Cap for logistic growth (max 100%, min 0%)
        df['cap'] = 100.0
        df['floor'] = 0.0

        # Add temperature as regressor
        if temperature_history is not None:
            temp_df = temperature_history.rename(columns={'timestamp': 'ds'})
            df = df.merge(temp_df, on='ds', how='left')
            df['avg_temp_c'] = df['avg_temp_c'].fillna(method='ffill')

        return df

    def fit(self, df: pd.DataFrame):
        """
        Train Prophet model.

        Args:
            df: Prepared DataFrame (ds, y, optional regressors)
        """
        # Add regressors if present
        if 'avg_temp_c' in df.columns:
            self.model.add_regressor('avg_temp_c', prior_scale=0.5)

        self.model.fit(df)
        self.trained = True
        print("Model trained successfully")

    def forecast(
        self,
        periods: int,
        freq: str = 'D',
        include_history: bool = True,
        temperature_forecast: pd.DataFrame = None
    ) -> pd.DataFrame:
        """
        Generate SOH forecast.

        Args:
            periods: Number of periods to forecast
            freq: Frequency ('D' for daily, 'H' for hourly)
            include_history: Include historical fit
            temperature_forecast: Future temperature values

        Returns:
            DataFrame with columns: ds, yhat (predicted SOH), yhat_lower, yhat_upper
        """
        if not self.trained:
            raise ValueError("Model must be trained before forecasting")

        # Create future dataframe
        future = self.model.make_future_dataframe(
            periods=periods,
            freq=freq,
            include_history=include_history
        )

        # Add cap/floor for logistic growth
        future['cap'] = 100.0
        future['floor'] = 0.0

        # Add temperature regressor
        if temperature_forecast is not None and 'avg_temp_c' in self.model.extra_regressors:
            temp_df = temperature_forecast.rename(columns={'timestamp': 'ds'})
            future = future.merge(temp_df, on='ds', how='left')
            # Forward fill missing values
            future['avg_temp_c'] = future['avg_temp_c'].fillna(method='ffill')

        # Forecast
        forecast = self.model.predict(future)

        return forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']]

    def plot_forecast(self, forecast: pd.DataFrame):
        """Plot forecast with components."""
        fig1 = self.model.plot(forecast)
        plt.title('Battery SOH Forecast')
        plt.ylabel('SOH (%)')
        plt.xlabel('Date')

        fig2 = self.model.plot_components(forecast)
        plt.show()

    def cross_validate(
        self,
        df: pd.DataFrame,
        initial: str = '365 days',
        period: str = '90 days',
        horizon: str = '180 days'
    ) -> pd.DataFrame:
        """
        Perform time-series cross-validation.

        Args:
            df: Training data
            initial: Initial training period
            period: Spacing between cutoff dates
            horizon: Forecast horizon

        Returns:
            DataFrame with CV results
        """
        from prophet.diagnostics import cross_validation, performance_metrics

        self.fit(df)

        df_cv = cross_validation(
            self.model,
            initial=initial,
            period=period,
            horizon=horizon
        )

        df_metrics = performance_metrics(df_cv)
        print("Cross-Validation Metrics:")
        print(df_metrics[['horizon', 'mse', 'rmse', 'mae', 'mape']].head())

        return df_cv


# Example usage
if __name__ == "__main__":
    # Load historical SOH data
    soh_df = pd.read_csv('battery_soh_history.csv', parse_dates=['timestamp'])

    # Optional: Load temperature data
    temp_df = pd.read_csv('temperature_history.csv', parse_dates=['timestamp'])

    # Initialize forecaster
    forecaster = BatterySOHForecaster(growth='logistic')

    # Prepare data
    train_df = forecaster.prepare_data(soh_df, temp_df)

    # Train
    forecaster.fit(train_df)

    # Forecast 365 days ahead
    forecast = forecaster.forecast(periods=365, freq='D')

    # Plot
    forecaster.plot_forecast(forecast)

    # Print key predictions
    future_forecast = forecast[forecast['ds'] > soh_df['timestamp'].max()]
    print("\nSOH Predictions (next 12 months):")
    print(future_forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']].head(12))

    # Estimate EOL (80% SOH)
    eol_date = future_forecast[future_forecast['yhat'] <= 80.0]['ds'].iloc[0]
    print(f"\nEstimated End-of-Life Date (80% SOH): {eol_date}")

LSTM for Multi-Variate Forecasting

python
import torch
import torch.nn as nn
from torch.utils.data import Dataset, DataLoader
import numpy as np
import pandas as pd
from typing import Tuple

class BatteryMultivariateDataset(Dataset):
    """
    Dataset for multivariate battery forecasting.

    Features: voltage, current, temperature, SOC
    Target: Future SOH
    """

    def __init__(
        self,
        data: np.ndarray,
        target: np.ndarray,
        seq_len: int = 100,
        pred_len: int = 10
    ):
        """
        Args:
            data: (N_samples, N_features) array
            target: (N_samples,) array (SOH)
            seq_len: Input sequence length
            pred_len: Prediction horizon
        """
        self.seq_len = seq_len
        self.pred_len = pred_len

        # Create sequences
        self.sequences = []
        self.targets = []

        for i in range(len(data) - seq_len - pred_len + 1):
            self.sequences.append(data[i:i+seq_len])
            self.targets.append(target[i+seq_len+pred_len-1])

        self.sequences = torch.FloatTensor(np.array(self.sequences))
        self.targets = torch.FloatTensor(np.array(self.targets))

    def __len__(self):
        return len(self.sequences)

    def __getitem__(self, idx):
        return self.sequences[idx], self.targets[idx]


class BatteryLSTMForecaster(nn.Module):
    """
    Multi-layer LSTM for battery SOH forecasting.

    Architecture:
    - Encoder LSTM: Process input sequence
    - Attention: Weight important time steps
    - Decoder: Generate future predictions
    """

    def __init__(
        self,
        n_features: int,
        hidden_dim: int = 128,
        n_layers: int = 3,
        dropout: float = 0.3
    ):
        super().__init__()

        self.n_features = n_features
        self.hidden_dim = hidden_dim

        # LSTM layers
        self.lstm = nn.LSTM(
            input_size=n_features,
            hidden_size=hidden_dim,
            num_layers=n_layers,
            batch_first=True,
            dropout=dropout
        )

        # Attention mechanism
        self.attention = nn.Sequential(
            nn.Linear(hidden_dim, hidden_dim),
            nn.Tanh(),
            nn.Linear(hidden_dim, 1)
        )

        # Output layer
        self.fc = nn.Sequential(
            nn.Linear(hidden_dim, 64),
            nn.ReLU(),
            nn.Dropout(dropout),
            nn.Linear(64, 1)
        )

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        # x: (batch, seq_len, n_features)

        # LSTM encoding
        lstm_out, _ = self.lstm(x)
        # lstm_out: (batch, seq_len, hidden_dim)

        # Attention weights
        attention_weights = self.attention(lstm_out)
        # attention_weights: (batch, seq_len, 1)
        attention_weights = torch.softmax(attention_weights, dim=1)

        # Weighted sum
        context = torch.sum(lstm_out * attention_weights, dim=1)
        # context: (batch, hidden_dim)

        # Prediction
        out = self.fc(context)
        return out.squeeze()


def train_lstm_forecaster(
    train_data: np.ndarray,
    train_target: np.ndarray,
    val_data: np.ndarray,
    val_target: np.ndarray,
    n_features: int,
    seq_len: int = 100,
    pred_len: int = 10,
    epochs: int = 100,
    batch_size: int = 64,
    lr: float = 0.001
) -> BatteryLSTMForecaster:
    """
    Train LSTM forecaster.

    Args:
        train_data: Training features (N, n_features)
        train_target: Training SOH values (N,)
        val_data: Validation features
        val_target: Validation SOH values
        n_features: Number of input features
        seq_len: Input sequence length
        pred_len: Prediction horizon (time steps ahead)
        epochs: Training epochs
        batch_size: Batch size
        lr: Learning rate

    Returns:
        Trained model
    """
    device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')

    # Create datasets
    train_dataset = BatteryMultivariateDataset(
        train_data, train_target, seq_len, pred_len
    )
    val_dataset = BatteryMultivariateDataset(
        val_data, val_target, seq_len, pred_len
    )

    train_loader = DataLoader(train_dataset, batch_size=batch_size, shuffle=True)
    val_loader = DataLoader(val_dataset, batch_size=batch_size)

    # Model
    model = BatteryLSTMForecaster(n_features=n_features).to(device)
    optimizer = torch.optim.Adam(model.parameters(), lr=lr)
    criterion = nn.MSELoss()

    best_val_loss = float('inf')
    patience = 15
    patience_counter = 0

    for epoch in range(epochs):
        # Training
        model.train()
        train_loss = 0
        for sequences, targets in train_loader:
            sequences = sequences.to(device)
            targets = targets.to(device)

            optimizer.zero_grad()
            outputs = model(sequences)
            loss = criterion(outputs, targets)
            loss.backward()

            # Gradient clipping
            torch.nn.utils.clip_grad_norm_(model.parameters(), max_norm=1.0)

            optimizer.step()
            train_loss += loss.item()

        train_loss /= len(train_loader)

        # Validation
        model.eval()
        val_loss = 0
        with torch.no_grad():
            for sequences, targets in val_loader:
                sequences = sequences.to(device)
                targets = targets.to(device)
                outputs = model(sequences)
                loss = criterion(outputs, targets)
                val_loss += loss.item()

        val_loss /= len(val_loader)

        if (epoch + 1) % 10 == 0:
            print(f"Epoch {epoch+1}/{epochs} - "
                  f"Train Loss: {train_loss:.4f}, Val Loss: {val_loss:.4f}")

        # Early stopping
        if val_loss < best_val_loss:
            best_val_loss = val_loss
            patience_counter = 0
            best_model_state = model.state_dict()
        else:
            patience_counter += 1
            if patience_counter >= patience:
                print(f"Early stopping at epoch {epoch+1}")
                break

    # Load best model
    model.load_state_dict(best_model_state)
    print(f"Best validation loss: {best_val_loss:.4f}")

    return model


def forecast_multi_step(
    model: BatteryLSTMForecaster,
    initial_sequence: np.ndarray,
    n_steps: int,
    device: str = 'cpu'
) -> np.ndarray:
    """
    Multi-step ahead forecasting with recursive prediction.

    Args:
        model: Trained LSTM model
        initial_sequence: (seq_len, n_features) - starting sequence
        n_steps: Number of steps to forecast ahead
        device: 'cpu' or 'cuda'

    Returns:
        Forecasted SOH values (n_steps,)
    """
    model.eval()
    model.to(device)

    predictions = []
    current_seq = torch.FloatTensor(initial_sequence).unsqueeze(0).to(device)
    # current_seq: (1, seq_len, n_features)

    with torch.no_grad():
        for _ in range(n_steps):
            # Predict next value
            pred = model(current_seq)
            predictions.append(pred.item())

            # Update sequence (shift left, add prediction)
            # Note: This assumes last feature is SOH. Adjust as needed.
            new_step = current_seq[0, -1, :].clone()
            new_step[-1] = pred  # Update SOH feature

            # Shift sequence
            current_seq = torch.cat([
                current_seq[:, 1:, :],
                new_step.unsqueeze(0).unsqueeze(0)
            ], dim=1)

    return np.array(predictions)

Energy Consumption Forecasting

python
from sklearn.ensemble import GradientBoostingRegressor
import numpy as np
import pandas as pd
from typing import Dict

class EnergyConsumptionForecaster:
    """
    Forecast vehicle energy consumption based on route and conditions.

    Features:
    - Route profile (distance, elevation change, road type)
    - Weather (temperature, wind, precipitation)
    - Traffic (congestion level, stop frequency)
    - Vehicle state (SOC, weight, tire pressure)
    - Driver behavior (historical consumption)
    """

    def __init__(self):
        self.model = GradientBoostingRegressor(
            n_estimators=200,
            learning_rate=0.05,
            max_depth=6,
            min_samples_split=20,
            random_state=42
        )

    def engineer_route_features(self, route_data: pd.DataFrame) -> pd.DataFrame:
        """
        Extract features from route data.

        Expected columns:
        - distance_km, elevation_gain_m, avg_speed_kmh,
        - road_type, ambient_temp_c, wind_speed_kmh,
        - traffic_level, stop_count
        """
        features = pd.DataFrame()

        # Route characteristics
        features['distance_km'] = route_data['distance_km']
        features['elevation_gain_m'] = route_data['elevation_gain_m']
        features['avg_speed_kmh'] = route_data['avg_speed_kmh']

        # Energy-relevant derived features
        features['elevation_per_km'] = (
            route_data['elevation_gain_m'] / (route_data['distance_km'] + 1e-3)
        )

        # Speed efficiency (most efficient around 60 km/h)
        features['speed_efficiency'] = np.exp(-((route_data['avg_speed_kmh'] - 60) / 30) ** 2)

        # Road type encoding
        if 'road_type' in route_data.columns:
            road_dummies = pd.get_dummies(route_data['road_type'], prefix='road')
            features = pd.concat([features, road_dummies], axis=1)

        # Weather impact
        features['ambient_temp_c'] = route_data['ambient_temp_c']
        features['temp_deviation_from_20'] = abs(route_data['ambient_temp_c'] - 20)
        features['wind_speed_kmh'] = route_data.get('wind_speed_kmh', 0)

        # Traffic impact
        features['traffic_level'] = route_data.get('traffic_level', 0)
        features['stops_per_km'] = (
            route_data.get('stop_count', 0) / (route_data['distance_km'] + 1e-3)
        )

        return features

    def train(self, X: pd.DataFrame, y: np.ndarray):
        """
        Train energy consumption model.

        Args:
            X: Feature matrix
            y: Energy consumed [kWh]
        """
        self.model.fit(X, y)
        score = self.model.score(X, y)
        print(f"Training R2: {score:.3f}")

    def predict_consumption(self, X: pd.DataFrame) -> np.ndarray:
        """
        Predict energy consumption for routes.

        Returns:
            Energy consumption [kWh]
        """
        return self.model.predict(X)

    def predict_range(
        self,
        route_features: pd.DataFrame,
        current_soc: float,
        battery_capacity_kwh: float
    ) -> Dict[str, float]:
        """
        Predict remaining range based on route profile.

        Args:
            route_features: Features for candidate routes
            current_soc: Current state of charge (0-100%)
            battery_capacity_kwh: Total battery capacity

        Returns:
            Dictionary with range estimates
        """
        # Available energy
        available_energy_kwh = (current_soc / 100) * battery_capacity_kwh

        # Predict consumption per route
        consumption_per_km = self.predict_consumption(route_features)

        # Range estimates
        estimated_range_km = available_energy_kwh / (consumption_per_km + 1e-6)

        return {
            'available_energy_kwh': available_energy_kwh,
            'estimated_consumption_kwh_per_km': consumption_per_km.mean(),
            'estimated_range_km': estimated_range_km.mean(),
            'range_uncertainty_km': estimated_range_km.std()
        }


# Example usage
if __name__ == "__main__":
    # Load historical trip data
    trips_df = pd.read_csv('historical_trips.csv')

    # Feature engineering
    forecaster = EnergyConsumptionForecaster()
    X = forecaster.engineer_route_features(trips_df)
    y = trips_df['energy_consumed_kwh'].values

    # Train
    forecaster.train(X, y)

    # Predict for new route
    new_route = pd.DataFrame([{
        'distance_km': 150,
        'elevation_gain_m': 500,
        'avg_speed_kmh': 80,
        'road_type': 'highway',
        'ambient_temp_c': 15,
        'wind_speed_kmh': 20,
        'traffic_level': 2,
        'stop_count': 5
    }])

    route_features = forecaster.engineer_route_features(new_route)
    range_estimate = forecaster.predict_range(
        route_features,
        current_soc=70,
        battery_capacity_kwh=75
    )

    print("\nRange Estimate:")
    for key, value in range_estimate.items():
        print(f"  {key}: {value:.2f}")

Charging Demand Forecasting (SARIMAX)

python
from statsmodels.tsa.statespace.sarimax import SARIMAX
import pandas as pd
import numpy as np
from typing import Tuple

class ChargingDemandForecaster:
    """
    Forecast charging station demand using SARIMAX.

    Captures:
    - Hourly seasonality (peak hours)
    - Weekly seasonality (weekday vs weekend)
    - Exogenous factors (weather, events)
    """

    def __init__(
        self,
        order: Tuple[int, int, int] = (1, 1, 1),
        seasonal_order: Tuple[int, int, int, int] = (1, 1, 1, 24)
    ):
        """
        Args:
            order: (p, d, q) for ARIMA
            seasonal_order: (P, D, Q, s) for seasonal component
                s=24 for hourly data with daily seasonality
        """
        self.order = order
        self.seasonal_order = seasonal_order
        self.model = None

    def fit(
        self,
        y: pd.Series,
        exog: pd.DataFrame = None
    ):
        """
        Fit SARIMAX model.

        Args:
            y: Time-series of charging demand (indexed by datetime)
            exog: Exogenous variables (temperature, events, etc.)
        """
        self.model = SARIMAX(
            y,
            exog=exog,
            order=self.order,
            seasonal_order=self.seasonal_order,
            enforce_stationarity=False,
            enforce_invertibility=False
        )

        self.results = self.model.fit(disp=False)
        print("Model fitted successfully")
        print(self.results.summary())

    def forecast(
        self,
        steps: int,
        exog: pd.DataFrame = None
    ) -> pd.DataFrame:
        """
        Forecast future demand.

        Args:
            steps: Number of steps ahead
            exog: Future exogenous variables

        Returns:
            DataFrame with forecast and confidence intervals
        """
        if self.model is None:
            raise ValueError("Model must be fitted before forecasting")

        forecast = self.results.get_forecast(steps=steps, exog=exog)
        forecast_df = forecast.summary_frame()

        return forecast_df.rename(columns={
            'mean': 'demand_forecast',
            'mean_ci_lower': 'lower_bound',
            'mean_ci_upper': 'upper_bound'
        })


# Example usage
if __name__ == "__main__":
    # Load hourly charging demand data
    demand_df = pd.read_csv('charging_demand.csv', parse_dates=['timestamp'])
    demand_df = demand_df.set_index('timestamp')

    # Exogenous variables
    exog_df = demand_df[['temperature_c', 'is_weekend', 'is_holiday']]
    y = demand_df['num_active_chargers']

    # Fit model
    forecaster = ChargingDemandForecaster(
        order=(1, 1, 1),
        seasonal_order=(1, 1, 1, 24)  # Daily seasonality
    )
    forecaster.fit(y, exog=exog_df)

    # Forecast next 48 hours
    future_exog = pd.DataFrame({
        'temperature_c': np.random.normal(20, 5, 48),
        'is_weekend': [0] * 24 + [1] * 24,
        'is_holiday': [0] * 48
    })

    forecast = forecaster.forecast(steps=48, exog=future_exog)
    print("\nDemand Forecast (next 48 hours):")
    print(forecast.head(12))

Deployment Strategy

yaml
# Time-series forecasting pipeline
pipeline:
  data_preparation:
    resampling: Hourly / Daily aggregation
    missing_data: Forward fill + interpolation
    outlier_handling: Median filter + IQR clipping

  feature_engineering:
    lag_features: [1, 7, 30] days
    rolling_statistics: Mean, std over [7, 30] day windows
    calendar_features: Hour, day_of_week, month, is_holiday

  model_training:
    framework: Prophet (quick) / LSTM (complex)
    validation: Time-series cross-validation
    hyperparameter_tuning: Bayesian optimization

  inference:
    frequency: Daily batch for long-term, real-time for short-term
    latency: <1s for API requests

  monitoring:
    metrics:
      - MAPE (Mean Absolute Percentage Error)
      - Coverage (% actual within confidence interval)
      - Forecast bias (systematic over/under prediction)

  retraining:
    trigger: MAPE > 10% OR weekly schedule
    validation: Backtest on last 90 days

Production Checklist

• ] Historical data spans multiple seasons/years (if seasonal)
• ] Missing data strategy validated (forward fill, interpolation)
• ] Outliers identified and handled (capping, removal)
• ] Stationarity tested (ADF test) and differencing applied if needed
• ] Cross-validation on time-series splits (not random)
• ] Confidence intervals calibrated (90% coverage target)
• ] Model performance monitored (MAPE, coverage)
• ] Forecast horizon validated against business needs
• ] Retraining pipeline automated with data quality checks
• ] A/B testing for model updates (gradual rollout)