Hybrid Machine Learning-Weibull for Remaining Useful Life Prediction

Mean Weibull Error (MWE)

Mean Weibull Error evaluates the cumulative distance of the predicted scale parameter \(\lambda\) from the actual failure point, effectively mapping standard deviation equivalents over survival profiles.

def weibull_error_dataset(x_test, x_pred, k):
    sum_error = 0
    for x, lmbda in zip(x_test, x_pred):
        sum_error += weibull_error(x, lmbda, k)
    return sum_error / len(x_test)

Mean Quadratic Weibull Error (MQWE)

MQWE enforces high quadratic penalties for extreme anomalies, evaluating the general reliability curve deviations on outliers. It is critical for safety-critical hardware where catastrophic failures must be avoided.

def quadratic_weibull_error_dataset(x_test, x_pred, k):
    sum_error = 0
    for x, lmbda in zip(x_test, x_pred):
        sum_error += weibull_error(x, lmbda, k) ** 2
    return sum_error / len(x_test)

Mean Directional Errors (MRWE & MLWE)

Mean Left Weibull Error (MLWE) and Mean Right Weibull Error (MRWE) are used to isolate whether a model is aggressively over-predicting or under-predicting the failure threshold.

def left_weibull_error(x, lmbda, k):
    lmbda = lmbda if lmbda > 0 else 1e-12
    s = 1.0 / k
    z = (x / lmbda) ** k
    gamma_s = gamma(s)
    l_gamma = gammaincc(s, z) * gamma_s
    l_gamma_0 = gammaincc(s, 0) * gamma_s
    return x + (lmbda * l_gamma / k) - (lmbda * l_gamma_0 / k)
 
def right_weibull_error(x, lmbda, k):
    lmbda = lmbda if lmbda > 0 else 1e-12
    s = 1.0 / k
    z = (x / lmbda) ** k
    gamma_s = gamma(s)
    l_gamma = gammaincc(s, z) * gamma_s
    return lmbda * l_gamma / k

Scipy Integration & Core Math

Calculations use Scipy's special incomplete gamma functions to evaluate the probability densities dynamically on every telemetry step.

from scipy.special import gamma, gammaincc
 
def weibull_error(x, lmbda, k):
    lmbda = lmbda if lmbda > 0 else 1e-12
    s = 1.0 / k
    z = (x / lmbda) ** k
    gamma_s = gamma(s)
    l_gamma = gammaincc(s, z) * gamma_s
    l_gamma_0 = gammaincc(s, 0) * gamma_s
    return x + (2 * lmbda * l_gamma / k) - (lmbda * l_gamma_0 / k)

Interactive Weibull Error Explorer

Adjust the telemetry failure target and the model's predicted parameters to see how errors are distributed.

IoT Edge Gateway & Telemetry

Data is captured on the vehicle using an ESP32-S3 microcontroller acting as an edge gateway. Live sensor packages collect TPMS data (tyre pressures and temperatures), OBD-II diagnostics (engine load and current speeds), and spatial GPS logs to track high-accuracy mileage.

Real-World Data Engineering & Preprocessing

Real-world telemetry gathered from IoT microcontrollers is highly asynchronous, noisy, and prone to packet drops. A comprehensive mathematical and data preprocessing pipeline was engineered to clean, smooth, and transform the raw sensor data:

• 1. Spatial Outlier Detection (IQR Bounds): Erratic coordinate jumps due to temporary GPS signal loss or multipath propagation are detected and filtered using the Interquartile Range (IQR) method:
  \[IQR = Q_3 - Q_1 \quad \implies \quad \text{Valid Range} = [Q_1 - 1.5 \times IQR, \; Q_3 + 1.5 \times IQR]\]

  def delete_outlier_gps(df):
      df_copy = df.copy()
      for col_name in ['latitude', 'longitude']:
          Q1 = df.loc[:, col_name].quantile(0.25)
          Q3 = df.loc[:, col_name].quantile(0.75)
          IQR = Q3 - Q1
          lower_bound = Q1 - 1.5 * IQR
          upper_bound = Q3 + 1.5 * IQR
          df_copy = df_copy[(df_copy[col_name] >= lower_bound) & (df_copy[col_name] <= upper_bound)]
          df_copy.reset_index(drop=True, inplace=True)
      return df_copy

• 2. Geodesic Distance Tracking (Haversine Formula): Spherical distance (meters) between consecutive spatial coordinates is computed using the Haversine formula (Earth's radius \(R \approx 6,371,000\text{ m}\)):
  \[a = \sin^2\left(\frac{\Delta\phi}{2}\right) + \cos(\phi_1)\cos(\phi_2)\sin^2\left(\frac{\Delta\theta}{2}\right)\] \[d = 2R \cdot \operatorname{arctan2}\left(\sqrt{a}, \sqrt{1-a}\right)\]
  Packets with anomalous jumps (\(d_i > 1000\text{ m}\)) are discarded. The cumulative distance traveled is calculated using a rolling prefix sum (\(\sum d_i\)).

  def change_to_distance(df):
      df_copy = df.copy()
      df_copy.reset_index(drop=True, inplace=True)
      for i in range(1, len(df_copy)):
          lat1, lon1 = df_copy.loc[i-1, 'latitude'], df_copy.loc[i-1, 'longitude']
          lat2, lon2 = df_copy.loc[i, 'latitude'], df_copy.loc[i, 'longitude']
          distance_km = haversine((lat1, lon1), (lat2, lon2), unit=Unit.METERS)
          df_copy.loc[i, 'distance_traveled_m'] = distance_km
      df_copy.at[0, 'distance_traveled_m'] = 0
      df_copy['sum_distance_traveled_m'] = df_copy['distance_traveled_m'].cumsum()
      return df_copy

• 3. Telemetry Smearing & Packet Drop Compensation: Loss of sensor telemetry packets (which register as 0 values in TPMS outputs) are converted to NaN and then filled sequentially. We apply a backward fill (bfill) to resolve initial zero readings, followed by a forward fill (ffill) to propagate the last known telemetry:
  \[\mathbf{x}_{\text{smeared}} = \operatorname{ffill}\Big(\operatorname{bfill}\big(\{x_i \mid x_i \neq 0, \text{ else } \text{NaN}\}\big)\Big)\]

  def smear_temp_pressure(df):
      df_copy = df.copy()
      for col_prefix in ['fl', 'fr', 'rl', 'rr']:
          for suffix in ['temp', 'pressure']:
              col_name = f'{col_prefix}_{suffix}'
              df_copy.loc[df_copy[col_name] == 0, col_name] = np.nan
              df_copy[col_name] = df_copy[col_name].bfill().ffill()
      return df_copy

• 4. Vertical Load Stressor (Elevation Delta): Vertical elevation adjustments between snapshots are computed to capture grade changes and additional mechanical loading on tires:
  \[\Delta E_i = H_i - H_{i-1} \quad (H = \text{altitude in meters})\]

  def change_to_elevation(df):
      df_copy = df.copy()
      df_copy.reset_index(drop=True, inplace=True)
      for i in range(1, len(df_copy)):
          alt_1 = df_copy.loc[i-1, 'altitude']
          alt_2 = df_copy.loc[i, 'altitude']
          df_copy.loc[i, 'elevation_traveled_m'] = alt_2 - alt_1
      df_copy.at[0, 'elevation_traveled_m'] = 0
      return df_copy

Model Optimization & Deployment

Models are trained to predict the dynamic \(\lambda\) parameter. LightGBM achieved a Mean Absolute Percentage Error (MAPE) of 0.0025 and Mean Absolute Scaled Error (MASE) of 0.2168. Models are then serialized for cloud server deployment.

from sklearn.model_selection import GridSearchCV
import lightgbm as lgb
import pickle

lightgbm_hyperparam = {
    'num_leaves': [20, 31, 50, 80],
    'min_child_samples': [20, 50, 100],
    'max_depth': [-1, 5, 8],
    'force_col_wise': [True]
}
lightgbm_grid = GridSearchCV(estimator=lgb.LGBMRegressor(), param_grid=lightgbm_hyperparam, cv=3)
lightgbm_grid.fit(X_train, y_train)
with open('lightgbm_grid.pkl', 'wb') as f:
    pickle.dump(lightgbm_grid.best_estimator_, f)

System Demonstration Video

Watch the real-time telemetry streaming from the ESP32-S3 edge gateway to the cloud dashboard and dynamic predictive degradation curve updates in action: