1) Problem Statement¶
1.1) Business Context¶
Renewable energy sources play an increasingly important role in the global energy mix, as the effort to reduce the environmental impact of energy production increases.
Out of all the renewable energy alternatives, wind energy is one of the most developed technologies worldwide. The U.S Department of Energy has put together a guide to achieving operational efficiency using predictive maintenance practices.
Predictive maintenance uses sensor information and analysis methods to measure and predict degradation and future component capability. The idea behind predictive maintenance is that failure patterns are predictable and if component failure can be predicted accurately and the component is replaced before it fails, the costs of operation and maintenance will be much lower.
The sensors fitted across different machines involved in the process of energy generation collect data related to various environmental factors (temperature, humidity, wind speed, etc.) and additional features related to various parts of the wind turbine (gearbox, tower, blades, break, etc.).
1.2) Objective¶
“ReneWind” is a company working on improving the machinery/processes involved in the production of wind energy using machine learning and has collected data of generator failure of wind turbines using sensors. They have shared a ciphered version of the data, as the data collected through sensors is confidential (the type of data collected varies with companies). Data has 40 predictors, 20000 observations in the training set and 5000 in the test set.
The objective is to build various classification models, tune them, and find the best one that will help identify failures so that the generators could be repaired before failing/breaking to reduce the overall maintenance cost. The nature of predictions made by the classification model will translate as follows:
- True positives (TP) are failures correctly predicted by the model. These will result in repairing costs.
- False negatives (FN) are real failures where there is no detection by the model. These will result in replacement costs.
- False positives (FP) are detections where there is no failure. These will result in inspection costs.
It is given that the cost of repairing a generator is much less than the cost of replacing it, and the cost of inspection is less than the cost of repair.
“1” in the target variables should be considered as “failure” and “0” represents “No failure”.
1.3) Data Description¶
The data provided is a transformed version of the original data which was collected using sensors.
- Train.csv - To be used for training and tuning of models.
- Test.csv - To be used only for testing the performance of the final best model.
Both the datasets consist of 40 predictor variables and 1 target variable.
2) Project Setup¶
2.1) Install Libraries¶
In [ ]:
# 2.1) Install Libraries
!pip install --quiet \
numpy pandas scikit-learn \
tensorflow keras keras-tuner \
imbalanced-learn lime \
xgboost lightgbm catboost \
matplotlib seaborn
print("All required libraries have been successfully installed.")
2.2) Import Libraries¶
In [ ]:
# 2.2) Import Libraries
# Set a global random seed
RANDOM_STATE = 48
# Core & Utility imports
import logging
import os
import random
import sys
import warnings
from collections import Counter
from pathlib import Path
from google.colab import drive
from pathlib import Path
# Data Handling & Visualization
import joblib
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from IPython.display import display
# Preprocessing, modeling, and evaluation
import sklearn
from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.compose import ColumnTransformer
from sklearn.decomposition import PCA
from sklearn.ensemble import RandomForestClassifier
from sklearn.feature_selection import SelectKBest, mutual_info_classif
from sklearn.impute import SimpleImputer
from sklearn.linear_model import LogisticRegression
# metrics import
from sklearn.metrics import (ConfusionMatrixDisplay, accuracy_score,
classification_report, confusion_matrix, f1_score,
make_scorer, precision_recall_curve,
precision_score, recall_score, roc_auc_score,
roc_curve)
# model_selection import
from sklearn.model_selection import (GridSearchCV, StratifiedKFold,
cross_validate, train_test_split)
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import FunctionTransformer, StandardScaler
from sklearn.utils import class_weight
# For handling class imbalance and powerful gradient boosting models
from imblearn.over_sampling import SMOTE
from imblearn.pipeline import Pipeline as ImbPipeline
import catboost
import lightgbm as lgb
import xgboost as xgb
# For explaining model predictions
import shap
from lime.lime_tabular import LimeTabularExplainer
# For building neural networks
import keras_tuner as kt
import tensorflow as tf
from tensorflow import keras
from tensorflow.keras import backend as K
from tensorflow.keras import layers, regularizers
from tensorflow.keras.callbacks import EarlyStopping, ReduceLROnPlateau
# Ensure reproducibility
os.environ['TF_DETERMINISTIC_OPS'] = '1'
random.seed(RANDOM_STATE)
np.random.seed(RANDOM_STATE)
tf.random.set_seed(RANDOM_STATE)
keras.utils.set_random_seed(RANDOM_STATE)
tf.config.experimental.enable_op_determinism()
# Suppress warnings (use with caution)
warnings.filterwarnings('ignore', category=FutureWarning)
warnings.filterwarnings('ignore', category=UserWarning)
# Pandas display options
pd.set_option('display.max_columns', None)
pd.set_option('display.width', 1000)
# Logging configuration
logging.basicConfig(
level=logging.INFO,
format="%(asctime)s [%(levelname)s] %(message)s",
handlers=[logging.StreamHandler(sys.stdout)]
)
# Set Seaborn style after all imports
sns.set_style('whitegrid')
print("All libraries successfully imported.")
2.3) Setup Google Drive¶
In [ ]:
# 2.3) Setup Google Drive
# Mount Google Drive
drive.mount('/content/drive')
# Define a constant for your project's root directory in Drive
# This makes it easy to update if you move your project
PROJECT_ROOT = Path('/content/drive/MyDrive/Colab Notebooks/')
# Create specific paths for your data files
TRAIN_FILE = PROJECT_ROOT / 'Train.csv'
TEST_FILE = PROJECT_ROOT / 'Test.csv'
# Now, you can use these variables directly when you need them
# Example: df = pd.read_csv(TRAIN_FILE)
# Verify that the training file exists
if TRAIN_FILE.is_file():
print(f"Train.csv found at: {TRAIN_FILE}")
else:
print("Train.csv not found. Check the PROJECT_ROOT path.")
3) Data Loading & Overview¶
3.1) Data Loading¶
In [ ]:
# 3.1) Data Loading
# Load the original data
raw_train = pd.read_csv(TRAIN_FILE)
raw_test = pd.read_csv(TEST_FILE)
# Create copies to work with
train_df = raw_train.copy()
test_df = raw_test.copy()
# Display dataset shapes
print(f"Training data shape: {train_df.shape}")
print(f"Test data shape: {test_df.shape}")
3.2) Data Overview¶
In [ ]:
# 3.2) Data Overview
# Preview first few rows of training data
display(train_df.head())
# Basic info and data types
print("\n-----------------------------------------------------------------\n")
print("Training Data Info:")
train_df.info()
# Summary statistics for numerical features
print("\n-----------------------------------------------------------------\n")
print("Training Data Summary Statistics:")
display(train_df.describe())
print("\n-----------------------------------------------------------------\n")
# Check for missing values in train and test sets
def missing_values_summary(train_df, test_df):
"""Creates a summary of missing values for train and test sets."""
missing_train = train_df.isnull().sum()
missing_test = test_df.isnull().sum()
missing_summary = pd.DataFrame({
'Train Missing Count': missing_train[missing_train > 0],
'Train Missing %': (train_df.isnull().mean() * 100)[missing_train > 0],
'Test Missing Count': missing_test[missing_test > 0],
'Test Missing %': (test_df.isnull().mean() * 100)[missing_test > 0]
}).sort_values(by='Train Missing %', ascending=False)
# Fill NaN for display purposes if a column has missing values in one set but not the other
return missing_summary.fillna(0)
# Generate and display the summary
missing_summary = missing_values_summary(train_df, test_df)
print("Missing Values Summary:")
if missing_summary.empty:
print("No missing values found in either dataset. ✅")
else:
display(missing_summary)
print("\n-----------------------------------------------------------------\n")
# Visualize target class distribution
plt.figure(figsize=(8, 5))
# Create the plot
ax = sns.countplot(x='Target', data=train_df, palette='viridis')
# Add percentage labels
total = len(train_df)
for p in ax.patches:
percentage = f'{100 * p.get_height() / total:.1f}%'
x = p.get_x() + p.get_width() / 2
y = p.get_height()
ax.annotate(percentage, (x, y), ha='center', va='bottom', fontsize=12, weight='bold')
# Customize titles and labels
plt.title('Target Class Distribution (0: No Failure, 1: Failure)', fontsize=14, weight='bold')
plt.xlabel('Class', fontsize=12)
plt.ylabel('Count', fontsize=12)
plt.xticks(ticks=[0, 1], labels=['No Failure', 'Failure'])
plt.tight_layout()
plt.show()
print("\n-----------------------------------------------------------------\n")
# Print the exact value counts and ratio
print("Target Variable Value Counts:")
print(train_df['Target'].value_counts())
print("\nRatio of No Failure to Failure:")
print(f"1:{train_df['Target'].value_counts()[0] / train_df['Target'].value_counts()[1]:.1f}")
OBSERVATIONS:
- 20K training rows and 5K test rows, each with 40 continuous predictors and a binary failure target.
- Only features V1 and V2 have missing values (18 each); all other predictors are fully observed.
- Failures (Target = 1) account for just 1 110 of 20 000 cases (~5.6%), so I will need to address imbalance.
- All predictors are floating-point measurements, ready for scaling or normalization.
4) Exploratory Data Analysis¶
4.1) Univariate EDA¶
In [ ]:
# 4.1) Univariate Analysis
# Define predictor variables (all columns except 'Target')
predictors = train_df.columns.drop('Target')
# Set up the plotting area
fig, axes = plt.subplots(8, 5, figsize=(20, 28))
axes = axes.flatten() # Flatten the 2D array of axes for easy iteration
# Iterate through each predictor and plot its distribution colored by the Target
for i, col in enumerate(predictors):
sns.histplot(
data=train_df,
x=col,
hue='Target', # Color bars based on the 'Target' value
multiple="stack", # Stack the bars for each class
kde=True, # Overlay a Kernel Density Estimate curve
ax=axes[i],
palette='viridis', # Use a color-blind friendly palette
legend=False # Hide individual legends to avoid clutter
)
axes[i].set_title(f'Distribution of {col}', fontsize=12)
axes[i].set_xlabel('') # Remove x-axis labels for a cleaner look
axes[i].set_ylabel('') # Remove y-axis labels
# Add a clear, overarching title
plt.suptitle('Feature Distributions by Target Class (0 vs 1)', fontsize=22, y=1.01)
# Create a single legend for the entire figure
handles, labels = axes[0].get_legend_handles_labels()
fig.legend(handles, ['No Failure', 'Failure'], loc='upper right', title='Class')
# Adjust layout to prevent titles and plots from overlapping
plt.tight_layout(rect=[0, 0, 0.95, 1])
plt.show()
OBSERVATIONS:
- The data from each sensor has a unique shape and pattern, with some showing a wider range of values than others.
- The sensors record data on completely different scales, which confirms we need to normalize them before feeding them into our models.
- For several key sensors, the readings for failing generators are noticeably different from healthy ones. These are our most promising clues for predicting failure.
- Features like V10, V14, and V17 show a strong link to failures when their values are high, making them very important for our analysis.
4.2) Bivariate EDA¶
In [ ]:
# 4.2) Bivariate EDA
# Calculate the correlation of all features with the 'Target'
target_corr = train_df.corr()['Target'].drop('Target').sort_values(ascending=False)
# Plotting the correlations
plt.figure(figsize=(10, 8))
sns.barplot(x=target_corr.values, y=target_corr.index, palette='coolwarm')
plt.title('Feature Correlation with Target', fontsize=16)
plt.xlabel('Pearson Correlation Coefficient', fontsize=12)
plt.ylabel('Features', fontsize=12)
plt.axvline(x=0, color='grey', linestyle='--') # Add a line at zero for reference
plt.tight_layout()
plt.show()
# Inter-Feature Correlation Heatmap
print("\n-----------------------------------------------------------------\n")
print("--- Full Correlation Matrix Heatmap ---")
print("\n-----------------------------------------------------------------\n")
# Calculate the full correlation matrix
corr_matrix = train_df.drop('Target', axis=1).corr()
# Create a mask to hide the upper triangle (since the matrix is symmetric)
mask = np.triu(np.ones_like(corr_matrix, dtype=bool))
# Set up the matplotlib figure
plt.figure(figsize=(10, 8))
sns.heatmap(
corr_matrix,
mask=mask,
annot=False,
cmap='coolwarm',
vmin=-1,
vmax=1
)
plt.title('Inter-Feature Correlation Heatmap', fontsize=18)
plt.show()
OBSERVATIONS:
- A few sensors are strong red flags. High readings on V17, V10, and V38, and low readings on V11 and V21, are the clearest signals of a potential generator failure.
- Some Sensors Move Together: The heatmap shows that certain groups of sensors provide similar information, meaning they are likely monitoring related parts of the turbine.
- Confident Direction: These findings confirm our earlier visual analysis and give us a reliable, ranked list of the most important features for predicting failure.
5) Preprocessing & Feature Engineering¶
- Split our data into training and validation sets.
- Use stratification during the split to guarantee both the training and validation sets have the same percentage of failures.
- Build a Pipeline that bundles our two preprocessing steps: filling missing values and scaling features.
- The pipeline learns how to process the data using only the training set. It then applies these learned rules to the validation set.
5.1) Data Splitting¶
In [ ]:
# 5.1) Data Splitting
# Separate predictors (X) from the target variable (y)
X = train_df.drop('Target', axis=1)
y = train_df['Target']
# Split the data into training and validation sets
# We use stratify=y to ensure that the proportion of failures is the same in both sets.
# This is crucial for imbalanced datasets.
X_train, X_val, y_train, y_val = train_test_split(
X, y,
test_size=0.2, # Using 20% of the data for validation
random_state=RANDOM_STATE, # Ensures reproducibility
stratify=y # <<< Critical for maintaining class distribution
)
# --- Verification ---
# Print the shapes of the resulting datasets to confirm the split
print("Data splitting complete.")
print(f"Shape of X_train: {X_train.shape}")
print(f"Shape of X_val: {X_val.shape}")
print(f"Shape of y_train: {y_train.shape}")
print(f"Shape of y_val: {y_val.shape}")
# Verify the class distribution in the training and validation sets
print("\n--- Class Distribution Verification ---")
print(f"Failure percentage in training set: {y_train.mean() * 100:.2f}%")
print(f"Failure percentage in validation set: {y_val.mean() * 100:.2f}%")
5.2) Create Preprocessing Pipeline¶
In [ ]:
# 5.2) Create Preprocessing Pipeline
# Define the steps for the preprocessing pipeline
preprocessing_pipeline = Pipeline(steps=[
('imputer', SimpleImputer(strategy='median')),
('scaler', StandardScaler())
])
# This shows the structure of our preprocessing steps.
print("Preprocessing pipeline created successfully.")
display(preprocessing_pipeline)
5.3) Apply Pipeline to Data¶
In [ ]:
# 5.3) Apply Pipeline to Data
# Fit the pipeline on the training data and transform it
X_train_processed = preprocessing_pipeline.fit_transform(X_train)
# Use the already-fitted pipeline to transform the validation data
X_val_processed = preprocessing_pipeline.transform(X_val)
# Convert the processed numpy arrays back to DataFrames for easier inspection
X_train_processed = pd.DataFrame(X_train_processed, columns=X_train.columns)
X_val_processed = pd.DataFrame(X_val_processed, columns=X_val.columns)
print("Preprocessing complete. X_train and X_val have been transformed.")
print("\nFirst 5 rows of the processed training data:")
display(X_train_processed.head())
# Verify that the processed data has a mean near 0 and a standard deviation near 1
print("\nSummary statistics of the processed training data:")
display(X_train_processed.describe().round(2))
6) Model Building & Tuning¶
- I will use a custom Total Business Cost function as the primary metric to find the most economically valuable model.
- This function calculates a total cost by assigning a cost for each type of prediction error and a cost for correct repairs.
- This metric directly reflects the project's goal of minimizing maintenance costs, making it more useful than technical metrics like recall or accuracy alone. We will use it to compare all models and guide our hyperparameter tuning.
- Disclaimer: The costs used are informed estimates, as real-world figures are not available. Our framework is flexible, allowing for these values to be easily updated with precise accounting data for any future analysis.
6.1) Helper Functions¶
- calculate_business_cost: This function computes a single "total cost" number for a model's predictions. It takes the true and predicted labels, counts the TPs, FPs, and FNs, and multiplies them by their assigned business cost (FN=100, TP=5, FP=1).
- evaluate_model: For any model, it prints a full classification report, displays the total business cost calculated by the first function, and generates a clear confusion matrix plot to visualize the model's performance.
In [ ]:
# 6.1) Helper Functions
# This function calculates the total financial cost based on our custom cost matrix.
# This will be our primary metric for comparing models.
def calculate_business_cost(y_true, y_pred, cost_matrix):
"""Calculates the total business cost of a model's predictions."""
# Get the confusion matrix components
tn, fp, fn, tp = confusion_matrix(y_true, y_pred).ravel()
# Calculate total cost
total_cost = (
(tn * cost_matrix['TN']) +
(fp * cost_matrix['FP']) +
(fn * cost_matrix['FN']) +
(tp * cost_matrix['TP'])
)
return total_cost
# This function will print a full classification report, the total business cost,
# and plot a visually clear confusion matrix.
def evaluate_model(y_true, y_pred, cost_matrix, model_name="Model"):
"""
Prints a classification report, business cost, and plots a confusion matrix.
"""
# Print Classification Report
print(f"--- Evaluation Report for {model_name} ---")
print(classification_report(y_true, y_pred, target_names=['No Failure', 'Failure']))
# Calculate and Print Business Cost
cost = calculate_business_cost(y_true, y_pred, cost_matrix)
print(f"Total Business Cost: ${cost:,.2f}")
print("--------------------------------------------------\n")
# Plot Confusion Matrix
cm = confusion_matrix(y_true, y_pred)
disp = ConfusionMatrixDisplay(
confusion_matrix=cm,
display_labels=['No Failure', 'Failure']
)
disp.plot(cmap=plt.cm.Blues)
plt.title(f'Confusion Matrix - {model_name}')
plt.show()
# Define the cost matrix
COST_MATRIX = {
'TN': 0, # No cost for true negatives
'FP': 1, # Cost for a false alarm inspection
'TP': 5, # Cost for a correctly identified repair
'FN': 100 # Cost for a missed failure replacement
}
print("Helper functions and COST_MATRIX are defined and ready.")
6.2) Baseline Models¶
- Establish a performance benchmark by running 5 off-the-shelf baseline models (without hyperparameter tuning) to compare with neural network models.
In [ ]:
# 6.2) Baseline Model Training & Evaluation
# This dictionary holds the models that work well with the 'class_weight' parameter.
models_with_class_weight = {
'Logistic Regression': LogisticRegression(random_state=RANDOM_STATE, class_weight='balanced', solver='liblinear'),
'Random Forest': RandomForestClassifier(random_state=RANDOM_STATE, class_weight='balanced', n_jobs=-1)
}
# This dictionary holds the gradient boosting models.
boosting_models = {
'XGBoost': xgb.XGBClassifier(random_state=RANDOM_STATE, objective='binary:logistic', eval_metric='logloss'),
'LightGBM': lgb.LGBMClassifier(random_state=RANDOM_STATE),
'CatBoost': catboost.CatBoostClassifier(random_state=RANDOM_STATE, verbose=0) # verbose=0 keeps the output clean
}
# Store results for final comparison
model_performance = []
# Train and evaluate models using class_weight
for name, model in models_with_class_weight.items():
model.fit(X_train_processed, y_train)
y_pred = model.predict(X_val_processed)
cost = calculate_business_cost(y_val, y_pred, COST_MATRIX)
recall = recall_score(y_val, y_pred)
f1 = f1_score(y_val, y_pred)
model_performance.append({'Model': name, 'Total Business Cost': cost, 'Recall': recall, 'F1-Score': f1})
# Manually apply SMOTE and then train boosting models
smote = SMOTE(random_state=RANDOM_STATE)
X_train_smote, y_train_smote = smote.fit_resample(X_train_processed, y_train)
for name, model in boosting_models.items():
model.fit(X_train_smote, y_train_smote)
y_pred = model.predict(X_val_processed)
cost = calculate_business_cost(y_val, y_pred, COST_MATRIX)
recall = recall_score(y_val, y_pred)
f1 = f1_score(y_val, y_pred)
model_performance.append({'Model': name, 'Total Business Cost': cost, 'Recall': recall, 'F1-Score': f1})
# Final Comparison
print("--- Baseline Model Performance Summary ---")
performance_df = pd.DataFrame(model_performance).sort_values(by='Total Business Cost', ascending=True).reset_index(drop=True)
display(performance_df)
OBSERVATIONS:
- With a total cost of $3,614, LightGBM is the best-performing baseline model. It successfully identified over 88% of the actual failures, which is a solid result.
- The three gradient boosting models clearly outperformed the others. Overall, this sets a fairly high bar.
6.3) Neural Net Models¶
6.3.1) NN Model 1 - Shallow, regularized MLP w/ 2 hidden layers¶
In [ ]:
# 6.3.1) NN Model 1
# Define and Compile the Model
simple_nn = keras.Sequential([
layers.Input(shape=(X_train_processed.shape[1],)),
layers.Dense(64, activation='relu', kernel_regularizer=regularizers.l2(0.001)),
layers.Dropout(0.3),
layers.Dense(32, activation='relu', kernel_regularizer=regularizers.l2(0.001)),
layers.Dropout(0.3),
layers.Dense(1, activation='sigmoid')
])
simple_nn.compile(
optimizer=keras.optimizers.Adam(learning_rate=0.001),
loss='binary_crossentropy',
metrics=['accuracy']
)
print("--- Architecture of Baseline NN ---")
simple_nn.summary()
print("\n" + "="*50 + "\n")
# Calculate class weights to handle the imbalanced dataset
class_weights = class_weight.compute_class_weight(
'balanced',
classes=np.unique(y_train),
y=y_train
)
class_weights_dict = dict(enumerate(class_weights))
print("--- Starting Model Training ---")
history = simple_nn.fit(
X_train_processed,
y_train,
epochs=50,
batch_size=32,
validation_data=(X_val_processed, y_val),
class_weight=class_weights_dict,
verbose=0 # Set to 0 to show cleaner output
)
print("\n" + "="*50 + "\n")
# Visualize Training History
history_df = pd.DataFrame(history.history)
plt.figure(figsize=(10, 6))
plt.plot(history_df['loss'], label='Training Loss')
plt.plot(history_df['val_loss'], label='Validation Loss')
plt.title('Model Loss Over Epochs', fontsize=16)
plt.ylabel('Loss')
plt.xlabel('Epoch')
plt.legend(loc='upper right')
plt.grid(True)
plt.show()
print("\n" + "="*50 + "\n")
# Predict on the validation set and convert probabilities to classes
y_pred_proba = simple_nn.predict(X_val_processed).flatten()
y_pred_class = (y_pred_proba > 0.5).astype(int)
# Use our helper function for a full evaluation
evaluate_model(y_val, y_pred_class, COST_MATRIX, model_name="Simple Baseline NN")
# Append results for final comparison table
cost = calculate_business_cost(y_val, y_pred_class, COST_MATRIX)
recall = recall_score(y_val, y_pred_class)
f1 = f1_score(y_val, y_pred_class)
# Check if 'Simple NN' is already in the performance dataframe to avoid duplicates
if not any(d['Model'] == 'Simple NN' for d in model_performance):
model_performance.append({'Model': 'Simple NN', 'Total Business Cost': cost, 'Recall': recall, 'F1-Score': f1})
OBSERVATIONS:
- The model achieved a recall of 0.91 for the failure class, meaning it correctly identified 91% of the actual generator failures.
- Training history plot shows that both the training and validation loss decreased steadily and began to flatten out. This is a sign of a healthy, stable training process without significant overfitting.
- Most importantly, the model achieved a total cost of $3,025.
6.3.2) NN Model 2 - Deeper, more heavily regularized MLP w/ 3 hidden layers¶
In [ ]:
# 6.3.2) NN Model 2
# Define and Compile the Model
deeper_nn = keras.Sequential([
layers.Input(shape=(X_train_processed.shape[1],)),
layers.Dense(128, activation='relu', kernel_regularizer=regularizers.l2(0.001)),
layers.Dropout(0.4),
layers.Dense(64, activation='relu', kernel_regularizer=regularizers.l2(0.001)),
layers.Dropout(0.4),
layers.Dense(32, activation='relu', kernel_regularizer=regularizers.l2(0.001)),
layers.Dropout(0.4),
layers.Dense(1, activation='sigmoid')
])
deeper_nn.compile(
optimizer=keras.optimizers.Adam(learning_rate=0.001),
loss='binary_crossentropy',
metrics=['accuracy']
)
print("--- Architecture of Deeper NN ---")
deeper_nn.summary()
print("\n" + "="*50 + "\n")
# use the same class weights as before to handle imbalance
print("--- Starting Model Training ---")
history_deep = deeper_nn.fit(
X_train_processed,
y_train,
epochs=50,
batch_size=32,
validation_data=(X_val_processed, y_val),
class_weight=class_weights_dict,
verbose=0 # Set to 0 for a cleaner output
)
print("\n" + "="*50 + "\n")
# Visualize Training History
history_df_deep = pd.DataFrame(history_deep.history)
plt.figure(figsize=(10, 6))
plt.plot(history_df_deep['loss'], label='Training Loss')
plt.plot(history_df_deep['val_loss'], label='Validation Loss')
plt.title('Deeper NN: Model Loss Over Epochs', fontsize=16)
plt.ylabel('Loss')
plt.xlabel('Epoch')
plt.legend(loc='upper right')
plt.grid(True)
plt.show()
print("\n" + "="*50 + "\n")
# Evaluate the Model
y_pred_proba_deep = deeper_nn.predict(X_val_processed).flatten()
y_pred_class_deep = (y_pred_proba_deep > 0.5).astype(int)
evaluate_model(y_val, y_pred_class_deep, COST_MATRIX, model_name="Deeper NN")
# Append results for final comparison table
cost = calculate_business_cost(y_val, y_pred_class_deep, COST_MATRIX)
recall = recall_score(y_val, y_pred_class_deep)
f1 = f1_score(y_val, y_pred_class_deep)
model_performance.append({'Model': 'Deeper NN', 'Total Business Cost': cost, 'Recall': recall, 'F1-Score': f1})
OBSERVATIONS:
- New best total cost of $2,940.
- The model maintained an strong recall of 0.91 while only slightly decreasing precision to 0.89. This means we are still catching the vast majority of failures without a significant increase in false alarms.
- The history plot shows another healthy learning curve.
6.4) Hyperparameter Tuning with KerasTuner¶
In [ ]:
# 6.4) Hyperparameter Tuning with KerasTuner
# Define the Model-Building Function
def build_model(hp):
model = keras.Sequential()
model.add(layers.Input(shape=(X_train_processed.shape[1],)))
hp_units_1 = hp.Int('units_1', min_value=64, max_value=256, step=32)
model.add(layers.Dense(units=hp_units_1, activation='relu', kernel_regularizer=regularizers.l2(0.001)))
model.add(layers.Dropout(0.4))
hp_units_2 = hp.Int('units_2', min_value=32, max_value=128, step=32)
model.add(layers.Dense(units=hp_units_2, activation='relu', kernel_regularizer=regularizers.l2(0.001)))
model.add(layers.Dropout(0.4))
model.add(layers.Dense(32, activation='relu', kernel_regularizer=regularizers.l2(0.001)))
model.add(layers.Dropout(0.4))
model.add(layers.Dense(1, activation='sigmoid'))
hp_learning_rate = hp.Choice('learning_rate', values=[1e-2, 1e-3, 1e-4])
model.compile(
optimizer=keras.optimizers.Adam(learning_rate=hp_learning_rate),
loss='binary_crossentropy',
metrics=['accuracy']
)
return model
# Instantiate and Run the Tuner
tuner = kt.Hyperband(
build_model,
objective='val_loss',
max_epochs=25,
factor=3,
directory='keras_tuner_dir',
project_name='rene_wind_tuning',
overwrite=True
)
stop_early = EarlyStopping(monitor='val_loss', patience=5)
print("--- Starting Hyperparameter Search ---")
tuner.search(
X_train_processed,
y_train,
epochs=50,
validation_data=(X_val_processed, y_val),
class_weight=class_weights_dict,
callbacks=[stop_early],
verbose=1
)
print("--- Hyperparameter Search Complete ---\n")
print("\n" + "="*50 + "\n")
# Get the optimal hyperparameters
best_hps = tuner.get_best_hyperparameters(num_trials=1)[0]
print("--- Optimal Hyperparameters Found ---")
print(f"Optimal units in first layer: {best_hps.get('units_1')}")
print(f"Optimal units in second layer: {best_hps.get('units_2')}")
print(f"Optimal learning rate: {best_hps.get('learning_rate')}\n")
# Build the model with the best HPs
tuned_nn = tuner.hypermodel.build(best_hps)
# Retrain the best model
print("--- Training the Final Tuned Model ---")
history_tuned = tuned_nn.fit(
X_train_processed, y_train,
epochs=50,
validation_data=(X_val_processed, y_val),
class_weight=class_weights_dict,
callbacks=[stop_early],
verbose=0
)
print("--- Final Tuned Model Training Complete ---\n")
# Evaluate the Final Tuned Model
print("--- Evaluating the Final Tuned Model ---")
y_pred_proba_tuned = tuned_nn.predict(X_val_processed).flatten()
y_pred_class_tuned = (y_pred_proba_tuned > 0.5).astype(int)
evaluate_model(y_val, y_pred_class_tuned, COST_MATRIX, model_name="Tuned NN")
# Append results for final comparison table
cost = calculate_business_cost(y_val, y_pred_class_tuned, COST_MATRIX)
recall = recall_score(y_val, y_pred_class_tuned)
f1 = f1_score(y_val, y_pred_class_tuned)
model_performance.append({'Model': 'Tuned NN', 'Total Business Cost': cost, 'Recall': recall, 'F1-Score': f1})
# Build and Train the Best Model
# Get the optimal hyperparameters found by the tuner
best_hps = tuner.get_best_hyperparameters(num_trials=1)[0]
print(f"""
--- Optimal Hyperparameters Found ---
Optimal units in first layer: {best_hps.get('units_1')}
Optimal units in second layer: {best_hps.get('units_2')}
Optimal learning rate: {best_hps.get('learning_rate')}
""")
# Build the model with the optimal hyperparameters
tuned_nn = tuner.hypermodel.build(best_hps)
# Retrain this single best model on the data
print("\n--- Training the Final Tuned Model ---")
history_tuned = tuned_nn.fit(
X_train_processed,
y_train,
epochs=50, # Train for a sufficient number of epochs
validation_data=(X_val_processed, y_val),
class_weight=class_weights_dict,
callbacks=[stop_early], # Use early stopping to prevent overfitting
verbose=0
)
OBSERVATIONS:
- Tuning brought some success, and brought total cost down to $2,939, a minor improvement.
- Tuning shows the best-performing architecture used 160 neurons in the first hidden layer and 32 neurons in the second, with a learning rate of 0.001.
- This tuned model maintained an excellent recall of 0.91. This ensures it continues to catch the vast majority of generator failures, which is the primary goal.
6.5) Threshold Optimization¶
In [ ]:
# 6.5) Threshold Optimization
# Get Predictions from the Tuned Model
y_pred_proba = tuned_nn.predict(X_val_processed).flatten()
# Find the Optimal Threshold
thresholds = np.arange(0.01, 1.0, 0.01)
costs = []
for t in thresholds:
y_pred_class = (y_pred_proba > t).astype(int)
costs.append(calculate_business_cost(y_val, y_pred_class, COST_MATRIX))
optimal_idx = np.argmin(costs)
optimal_threshold = thresholds[optimal_idx]
min_cost = costs[optimal_idx]
print(f"Optimal Threshold for Tuned Model: {optimal_threshold:.2f}")
print(f"Minimum possible business cost: ${min_cost:,.2f}")
print("\n" + "="*50 + "\n")
# Evaluate the Tuned Model with the Optimal Threshold
y_pred_optimal = (y_pred_proba > optimal_threshold).astype(int)
evaluate_model(y_val, y_pred_optimal, COST_MATRIX, model_name="Tuned NN with Optimal Threshold")
# Append results for final comparison table
cost = calculate_business_cost(y_val, y_pred_optimal, COST_MATRIX)
recall = recall_score(y_val, y_pred_optimal)
f1 = f1_score(y_val, y_pred_optimal)
model_performance.append({'Model': 'Tuned NN (Optimal Threshold)', 'Total Business Cost': cost, 'Recall': recall, 'F1-Score': f1})
OBSERVATIONS:
- Threshold optimization was effective. We brought cost down to $2,816 with an optimal decision threshold of 0.24.
- Recall of 0.92 is strong.
7) Final Model Training & Evaluation¶
7.1) Train Final Model on Full Training Set¶
In [ ]:
# 7.1) Train Final Model on Full Training Set
# Prepare the Full Training Dat
X_full_train = raw_train.drop('Target', axis=1)
y_full_train = raw_train['Target']
# Use the same pipeline to process this full dataset.
X_full_train_processed = preprocessing_pipeline.fit_transform(X_full_train)
# Build the Final Model Architecture
print("--- Building the final champion model ---")
final_model = keras.Sequential([
layers.Input(shape=(X_full_train_processed.shape[1],)),
layers.Dense(160, activation='relu', kernel_regularizer=regularizers.l2(0.001)),
layers.Dropout(0.4),
layers.Dense(32, activation='relu', kernel_regularizer=regularizers.l2(0.001)),
layers.Dropout(0.4),
layers.Dense(32, activation='relu', kernel_regularizer=regularizers.l2(0.001)),
layers.Dropout(0.4),
layers.Dense(1, activation='sigmoid')
])
final_model.compile(
optimizer=keras.optimizers.Adam(learning_rate=0.001),
loss='binary_crossentropy',
metrics=['accuracy']
)
# Train the Model on All Data
print("--- Training the final model on all 20,000 training samples... ---")
final_model.fit(
X_full_train_processed,
y_full_train,
epochs=40,
batch_size=32,
class_weight=class_weights_dict,
verbose=0
)
print("--- Final Model is Trained and Ready for Evaluation ---")
7.2) Evaluate Final Model on Test Set¶
In [ ]:
# 7.2) Evaluate Final Model on Test Set
# Prepare the Test Data
X_test = raw_test.drop('Target', axis=1)
y_test = raw_test['Target']
# Use the pipeline that was already fitted on the full training data
X_test_processed = preprocessing_pipeline.transform(X_test)
# Make Predictions on the Test Set
print("--- Evaluating Final Model on Unseen Test Data ---")
y_pred_proba_final = final_model.predict(X_test_processed).flatten()
# Apply Optimal Threshold and Evaluate
optimal_threshold = 0.24
y_pred_final = (y_pred_proba_final > optimal_threshold).astype(int)
# Use our helper function for the final report
evaluate_model(y_test, y_pred_final, COST_MATRIX, model_name="Final Champion Model (on Test Set)")
# Display Final Model Comparison
print("\n--- Final Model Performance Comparison ---")
# This dataframe summarizes the performance of all models we tested on the validation set
performance_df = pd.DataFrame(model_performance).sort_values(by='Total Business Cost', ascending=True).reset_index(drop=True)
display(performance_df)
OBSERVATIONS
- The final model performed well on the test data, successfully identifying 88% of all generator failures.
- As expected, the model's performance was slightly lower on the test set than on the validation set. The final business cost of $4,639 reflects this.
- This evaluation provides the most honest and reliable measure of how the model would perform in a real-world production environment. The results confirm that the model is a highly effective and valuable tool.
8) Feature Importance & Effects¶
8.1) Global Feature Importance¶
In [ ]:
# 8.1) Feature Importance
# Get the Baseline Score
y_pred_baseline = (final_model.predict(X_test_processed, verbose=0) > 0.29).astype(int)
baseline_accuracy = accuracy_score(y_test, y_pred_baseline)
print(f"Baseline Model Accuracy: {baseline_accuracy:.4f}")
# Manually Calculate Permutation Importance
print("\n--- Calculating Permutation Importance... ---")
importances = {}
for col in X_test.columns:
# Create a copy of the processed test data
X_test_permuted = X_test_processed.copy()
# Get the index of the current column
col_idx = X_test.columns.get_loc(col)
# Shuffle the values in that single column
np.random.shuffle(X_test_permuted[:, col_idx])
# Get the new score with the shuffled data
y_pred_permuted = (final_model.predict(X_test_permuted, verbose=0) > 0.29).astype(int)
permuted_accuracy = accuracy_score(y_test, y_pred_permuted)
# The importance is the drop in score
feature_importance = baseline_accuracy - permuted_accuracy
importances[col] = feature_importance
print("--- Calculation Complete ---")
# Organize and Plot the Results
perm_importance = pd.Series(importances).sort_values(ascending=False)
print("\n--- Model Feature Importance ---")
plt.figure(figsize=(10, 8))
sns.barplot(x=perm_importance.values, y=perm_importance.index, palette='viridis')
plt.title('Feature Importance', fontsize=16)
plt.xlabel('Performance Drop (Decrease in Accuracy Score)')
plt.ylabel('Feature')
plt.tight_layout()
plt.show()
OBSERVATIONS:
- Top features identified here (V17, V10, V11, and V21) were identified as highly predictive during initial EDA, a strong sign that the model has learned the correct patterns from the data.
- There is a distinct group of about 7-8 highly impactful features at the top. After that, there is a noticeable drop-off in importance. This suggests that the model's performance is driven by a core set of reliable sensor signals.
- The model isn't relying on just one or two features. It has learned to weigh the input from several different sensors to make a decision, which makes it more robust and less likely to be thrown off by a single faulty reading.
- The features at the bottom of the list have an importance score near zero. This confirms that they contribute very little to the model's predictive power, and we could likely build a simpler, faster model by excluding them without a significant loss in performance.
8.2) Visualizing Feature Effects¶
In [ ]:
# 8.2) Visualizing Feature Effects
# Select Top Features
top_features = perm_importance.head(6).index.tolist()
print("--- Manually Calculating and Plotting Partial Dependence... ---")
# Create a figure to hold all the subplots
fig, axes = plt.subplots(2, 3, figsize=(18, 10))
axes = axes.flatten()
# Loop Through Each Top Feature
for i, feature_name in enumerate(top_features):
# Get the integer index of the current feature
feature_idx = X_test.columns.get_loc(feature_name)
# Create a grid of 10 values for this feature, from its min to max in the test set
feature_grid = np.linspace(
X_test_processed[:, feature_idx].min(),
X_test_processed[:, feature_idx].max(),
num=10
)
# This will store the average predictions
avg_preds = []
# Calculate Average Prediction for Each Grid Value
for grid_val in feature_grid:
# Create a copy of the processed test data
X_temp = X_test_processed.copy()
# Replace the entire column for the current feature with the grid value
X_temp[:, feature_idx] = grid_val
# Get the model's predictions on this modified data
preds = final_model.predict(X_temp, verbose=0)
# Calculate the average prediction and store it
avg_preds.append(preds.mean())
# Plot the Results
axes[i].plot(feature_grid, avg_preds, marker='o', linestyle='-')
axes[i].set_title(f'Partial Dependence of {feature_name}')
axes[i].set_xlabel('Feature Value')
axes[i].set_ylabel('Avg. Predicted Probability')
axes[i].grid(True)
# Finalize the Figure
fig.suptitle('Partial Dependence Plots for Top 6 Features', fontsize=20)
plt.tight_layout(rect=[0, 0, 1, 0.97])
plt.show()
OBSERVATIONS:
- Plots reveal that the model has learned complex, non-linear patterns. For features like V12, V9, V22, and V31, the probability of failure is highest at low and high values, and lowest near the middle of their range. This suggests the model has learned that any significant deviation from a normal operating value for these sensors, in either direction, increases the risk of failure.
- For features V10 and V14, the relationship is more straightforward. As the value of these features increases, the predicted probability of failure also increases steadily. This confirms they are strong positive predictors, which aligns perfectly with our earlier EDA and feature importance findings.
9) Conclusions & Business Recommendations¶
9.1) Project Summary & Final Model Performance¶
- This project aimed to develop a predictive maintenance model to minimize operational costs for ReneWind by accurately identifying generator failures. Detailed EDA revealed a class imbalance and the necessity for feature scaling. To directly address the business goal, a custom cost function was defined that heavily penalized missed failures (False Negatives).
- Strong benchmarks were set using five standard machine learning models, with LightGBM proving to be the best baseline. A series of neural network architectures were then developed, finding that a deep, regularized MLP as the most effective design. This champion architecture was further improved through systematic hyperparameter tuning and decision threshold optimization.
- The final and best-performing model is the Tuned Neural Network with an optimized classification threshold of 0.29.
- The model's strong performance, particularly its 88% recall on unseen data, confirms that it is a robust and valuable tool capable of generalizing well to new scenarios.
9.2) Recommendations for ReneWind¶
Based on the successful development and rigorous evaluation of the predictive maintenance model, we recommend the following actions to leverage this tool and maximize its value for ReneWind.
- Models perform best when they are kept up-to-date. We recommend establishing a cycle to:
- Monitor: Continuously track the model's performance on new data to ensure its accuracy does not degrade over time.
- Retrain: Periodically retrain the model on a monthly or quarterly basis, incorporating new sensor data as it is collected. This will allow the model to adapt to any changes in equipment behavior or environmental conditions.
- The model's performance is heavily influenced by the quality and quantity of the data it's trained on. While the current model is strong, its accuracy can be improved even further by gathering more examples of failure events. We recommend prioritizing the collection and labeling of data from any future generator failures to enhance the dataset for future retraining cycles.
The explainability plots are not just diagnostic tools; they are valuable operational assets. We recommend that when the model flags a turbine for inspection, the report should include an explanation showing which specific sensor readings drove the prediction. This will:
- Build Trust: Help engineers and technicians trust the model's recommendations.
- Speed Up Diagnostics: Guide maintenance teams on where to focus their initial inspection efforts.
By implementing these recommendations, ReneWind can fully capitalize on the predictive power of this model, leading to significant and sustained reductions in maintenance costs and operational downtime.
Export to HTML¶
In [ ]:
# Exporting to HTML
!jupyter nbconvert --to html "/content/drive/MyDrive/Colab Notebooks/Burns_Greg_ReneWind_Main_Project.ipynb"