Chapter 4: Predictive Maintenance and Troubleshooting

📖 Chapter Overview

In food manufacturing processes, sudden equipment failures and quality issues are serious problems that can lead to production shutdowns and product waste. In this chapter, you will learn how to predict equipment failures in advance and achieve planned maintenance using AI-powered Predictive Maintenance technology. We will also explain root cause analysis and decision support systems for troubleshooting when anomalies occur.

🎯 Learning Objectives

• Machine learning modeling techniques for equipment failure prediction
• Remaining Useful Life (RUL) estimation algorithms
• Anomaly detection and early warning system construction
• Data analysis techniques for Root Cause Analysis (RCA)
• Troubleshooting decision support systems
• Maintenance planning optimization methods

🔧 4.1 Fundamentals of Equipment Failure Prediction

Three Approaches to Predictive Maintenance

Predictive maintenance has three main approaches:

• Time-Based Maintenance (TBM): Maintenance on a fixed schedule (traditional approach)
• Condition-Based Maintenance (CBM): Threshold management based on sensor data
• Predictive Maintenance (PdM): Failure prediction using machine learning and optimal maintenance timing determination

Definition of Remaining Useful Life (RUL)

💻 Code Example 4.1: Equipment Failure Prediction with Random Forest

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import classification_report, confusion_matrix, roc_curve, auc
import warnings
warnings.filterwarnings('ignore')

plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# Generate equipment sensor data (1000 samples)
np.random.seed(42)
n_samples = 1000

# Normal operation data (70%)
n_normal = int(n_samples * 0.7)
temp_normal = np.random.normal(75, 5, n_normal)  # Temperature (°C)
vibration_normal = np.random.normal(0.3, 0.05, n_normal)  # Vibration (mm/s)
pressure_normal = np.random.normal(2.0, 0.2, n_normal)  # Pressure (MPa)
current_normal = np.random.normal(15, 2, n_normal)  # Current (A)
runtime_normal = np.random.uniform(0, 5000, n_normal)  # Runtime (h)

# Pre-failure data (30%)
n_failure = n_samples - n_normal
temp_failure = np.random.normal(90, 8, n_failure)  # High temperature trend
vibration_failure = np.random.normal(0.6, 0.1, n_failure)  # Increased vibration
pressure_failure = np.random.normal(2.5, 0.3, n_failure)  # Pressure increase
current_failure = np.random.normal(20, 3, n_failure)  # Current increase
runtime_failure = np.random.uniform(4000, 8000, n_failure)  # Long operation time

# Create dataframe
data = pd.DataFrame({
    'temperature': np.concatenate([temp_normal, temp_failure]),
    'vibration': np.concatenate([vibration_normal, vibration_failure]),
    'pressure': np.concatenate([pressure_normal, pressure_failure]),
    'current': np.concatenate([current_normal, current_failure]),
    'runtime_hours': np.concatenate([runtime_normal, runtime_failure]),
    'failure_risk': np.concatenate([np.zeros(n_normal), np.ones(n_failure)])
})

# Create derived features
data['temp_vibration_ratio'] = data['temperature'] / (data['vibration'] * 100)
data['power_consumption'] = data['current'] * data['pressure']  # Virtual power consumption
data['runtime_category'] = pd.cut(data['runtime_hours'], bins=5, labels=['Very Low', 'Low', 'Medium', 'High', 'Very High'])
data['runtime_category'] = data['runtime_category'].cat.codes

# Split training and test data
X = data.drop('failure_risk', axis=1)
y = data['failure_risk']
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42, stratify=y)

# Train Random Forest model
rf_model = RandomForestClassifier(
    n_estimators=100,
    max_depth=10,
    min_samples_split=10,
    min_samples_leaf=5,
    random_state=42,
    class_weight='balanced'
)
rf_model.fit(X_train, y_train)

# Prediction and evaluation
y_pred = rf_model.predict(X_test)
y_pred_proba = rf_model.predict_proba(X_test)[:, 1]

print("=" * 60)
print("Equipment Failure Prediction Model Evaluation")
print("=" * 60)
print("\nClassification Report:")
print(classification_report(y_test, y_pred, target_names=['Normal', 'Failure Risk']))

# Visualization
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Confusion matrix
cm = confusion_matrix(y_test, y_pred)
sns.heatmap(cm, annot=True, fmt='d', cmap='Greens', ax=axes[0, 0],
            xticklabels=['Normal', 'Failure Risk'], yticklabels=['Normal', 'Failure Risk'])
axes[0, 0].set_title('Confusion Matrix', fontsize=12, fontweight='bold')
axes[0, 0].set_ylabel('True Label')
axes[0, 0].set_xlabel('Predicted Label')

# ROC curve
fpr, tpr, thresholds = roc_curve(y_test, y_pred_proba)
roc_auc = auc(fpr, tpr)
axes[0, 1].plot(fpr, tpr, color='#11998e', lw=2, label=f'ROC Curve (AUC = {roc_auc:.3f})')
axes[0, 1].plot([0, 1], [0, 1], color='gray', lw=1, linestyle='--', label='Random')
axes[0, 1].set_xlim([0.0, 1.0])
axes[0, 1].set_ylim([0.0, 1.05])
axes[0, 1].set_xlabel('False Positive Rate')
axes[0, 1].set_ylabel('True Positive Rate')
axes[0, 1].set_title('ROC Curve', fontsize=12, fontweight='bold')
axes[0, 1].legend(loc='lower right')
axes[0, 1].grid(alpha=0.3)

# Feature importance
feature_importance = pd.DataFrame({
    'feature': X.columns,
    'importance': rf_model.feature_importances_
}).sort_values('importance', ascending=True)

axes[1, 0].barh(feature_importance['feature'], feature_importance['importance'], color='#38ef7d')
axes[1, 0].set_xlabel('Importance')
axes[1, 0].set_title('Feature Importance', fontsize=12, fontweight='bold')
axes[1, 0].grid(axis='x', alpha=0.3)

# Failure risk score distribution
axes[1, 1].hist(y_pred_proba[y_test == 0], bins=30, alpha=0.6, label='Normal', color='green')
axes[1, 1].hist(y_pred_proba[y_test == 1], bins=30, alpha=0.6, label='Failure Risk', color='red')
axes[1, 1].axvline(x=0.5, color='black', linestyle='--', linewidth=1, label='Threshold (0.5)')
axes[1, 1].set_xlabel('Failure Risk Score')
axes[1, 1].set_ylabel('Frequency')
axes[1, 1].set_title('Failure Risk Score Distribution', fontsize=12, fontweight='bold')
axes[1, 1].legend()
axes[1, 1].grid(alpha=0.3)

plt.tight_layout()
plt.savefig('failure_prediction_analysis.png', dpi=300, bbox_inches='tight')
plt.show()

print(f"\nROC-AUC Score: {roc_auc:.4f}")
print("\nTop 3 Feature Importance:")
print(feature_importance.tail(3).to_string(index=False))

Implementation Points:

• Prediction model integrating sensor data (temperature, vibration, pressure, current)
• Utilization of historical information such as runtime as features
• Accuracy improvement through derived features (temperature-vibration ratio, power consumption)
• Addressing class imbalance (class_weight='balanced')
• ROC-AUC evaluation provides guidance for threshold adjustment

⏱️ 4.2 Remaining Useful Life (RUL) Estimation

Importance of RUL Estimation

If we can predict "when" a failure will occur, we can plan optimal maintenance timing. RUL estimation is a technology that estimates the remaining time until failure from the current equipment condition.

Degradation Curve Model

💻 Code Example 4.2: Remaining Useful Life Estimation with LSTM

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_absolute_error, mean_squared_error
import warnings
warnings.filterwarnings('ignore')

# TensorFlow/Keras import (required for actual implementation)
# from tensorflow.keras.models import Sequential
# from tensorflow.keras.layers import LSTM, Dense, Dropout

plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# Generate equipment degradation data (simulated RUL data)
np.random.seed(42)
n_cycles = 200  # Number of life cycles
time_steps = 100  # Time steps for each cycle

# Generate degradation patterns for 3 equipment units
def generate_degradation_data(n_cycles, time_steps, failure_threshold=20):
    """Generate equipment degradation data"""
    data_list = []

    for cycle_id in range(n_cycles):
        # Random degradation rate
        tau = np.random.uniform(60, 120)  # Degradation time constant
        noise_level = np.random.uniform(2, 5)

        # Health index calculation (exponential degradation + noise)
        t = np.arange(time_steps)
        health_index = 100 * np.exp(-t / tau) + np.random.normal(0, noise_level, time_steps)
        health_index = np.clip(health_index, 0, 100)

        # RUL calculation (time until health index falls below threshold)
        failure_time = np.where(health_index < failure_threshold)[0]
        if len(failure_time) > 0:
            failure_step = failure_time[0]
        else:
            failure_step = time_steps

        rul = np.maximum(0, failure_step - t)

        # Sensor data (changes according to health index)
        temperature = 60 + (100 - health_index) * 0.5 + np.random.normal(0, 2, time_steps)
        vibration = 0.2 + (100 - health_index) * 0.008 + np.random.normal(0, 0.05, time_steps)
        pressure = 1.8 + (100 - health_index) * 0.01 + np.random.normal(0, 0.1, time_steps)

        for step in range(time_steps):
            data_list.append({
                'cycle_id': cycle_id,
                'time_step': step,
                'temperature': temperature[step],
                'vibration': vibration[step],
                'pressure': pressure[step],
                'health_index': health_index[step],
                'RUL': rul[step]
            })

    return pd.DataFrame(data_list)

# Generate data
degradation_data = generate_degradation_data(n_cycles, time_steps)

# Data preprocessing
feature_columns = ['temperature', 'vibration', 'pressure', 'health_index']
scaler = MinMaxScaler()
degradation_data[feature_columns] = scaler.fit_transform(degradation_data[feature_columns])

# RUL normalization (0-1 range)
max_rul = degradation_data['RUL'].max()
degradation_data['RUL_normalized'] = degradation_data['RUL'] / max_rul

# Split training and test data
train_cycles = int(n_cycles * 0.8)
train_data = degradation_data[degradation_data['cycle_id'] < train_cycles]
test_data = degradation_data[degradation_data['cycle_id'] >= train_cycles]

print("=" * 60)
print("Remaining Useful Life (RUL) Estimation System")
print("=" * 60)
print(f"\nTotal cycles: {n_cycles}")
print(f"Training cycles: {train_cycles}")
print(f"Test cycles: {n_cycles - train_cycles}")
print(f"Time steps per cycle: {time_steps}")
print(f"\nRUL range: 0 ~ {max_rul:.0f} steps")

# Simple prediction (using moving average-based prediction instead of actual LSTM)
def simple_rul_prediction(data, window=10):
    """Simple RUL prediction based on moving average"""
    predictions = []
    actuals = []

    for cycle_id in data['cycle_id'].unique():
        cycle_data = data[data['cycle_id'] == cycle_id].copy()

        for i in range(window, len(cycle_data)):
            # Predict from average decline rate of health index over past window steps
            recent_health = cycle_data.iloc[i-window:i]['health_index'].values
            health_decline_rate = (recent_health[0] - recent_health[-1]) / window

            current_health = cycle_data.iloc[i]['health_index']

            # Estimate steps until health index reaches 20%
            if health_decline_rate > 0.0001:
                estimated_rul = max(0, (current_health - 0.2) / health_decline_rate)
            else:
                estimated_rul = max_rul  # If no degradation observed

            estimated_rul = min(estimated_rul, max_rul)  # Upper limit constraint

            predictions.append(estimated_rul / max_rul)  # Normalize
            actuals.append(cycle_data.iloc[i]['RUL_normalized'])

    return np.array(predictions), np.array(actuals)

# Predict on test data
y_pred, y_true = simple_rul_prediction(test_data, window=10)

# Evaluation metrics
mae = mean_absolute_error(y_true * max_rul, y_pred * max_rul)
rmse = np.sqrt(mean_squared_error(y_true * max_rul, y_pred * max_rul))

print(f"\nPrediction Performance:")
print(f"MAE (Mean Absolute Error): {mae:.2f} steps")
print(f"RMSE (Root Mean Square Error): {rmse:.2f} steps")

# Visualization
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Degradation curves for sample cycles
sample_cycles = [0, 1, 2]
for i, cycle_id in enumerate(sample_cycles):
    cycle_data = degradation_data[degradation_data['cycle_id'] == cycle_id]
    axes[0, 0].plot(cycle_data['time_step'], cycle_data['health_index'],
                    label=f'Cycle {cycle_id}', alpha=0.7)

axes[0, 0].axhline(y=0.2, color='red', linestyle='--', linewidth=1, label='Failure Threshold (20%)')
axes[0, 0].set_xlabel('Time Step')
axes[0, 0].set_ylabel('Health Index (Normalized)')
axes[0, 0].set_title('Equipment Degradation Curves (Sample)', fontsize=12, fontweight='bold')
axes[0, 0].legend()
axes[0, 0].grid(alpha=0.3)

# RUL prediction vs actual
axes[0, 1].scatter(y_true * max_rul, y_pred * max_rul, alpha=0.3, s=10, color='#11998e')
axes[0, 1].plot([0, max_rul], [0, max_rul], 'r--', lw=2, label='Ideal Prediction')
axes[0, 1].set_xlabel('Actual RUL (Steps)')
axes[0, 1].set_ylabel('Predicted RUL (Steps)')
axes[0, 1].set_title('RUL Prediction Accuracy', fontsize=12, fontweight='bold')
axes[0, 1].legend()
axes[0, 1].grid(alpha=0.3)

# Prediction error distribution
errors = (y_pred - y_true) * max_rul
axes[1, 0].hist(errors, bins=50, color='#38ef7d', alpha=0.7, edgecolor='black')
axes[1, 0].axvline(x=0, color='red', linestyle='--', linewidth=2, label='Zero Error')
axes[1, 0].set_xlabel('Prediction Error (Steps)')
axes[1, 0].set_ylabel('Frequency')
axes[1, 0].set_title('RUL Prediction Error Distribution', fontsize=12, fontweight='bold')
axes[1, 0].legend()
axes[1, 0].grid(alpha=0.3)

# RUL time series for specific cycle
test_cycle = test_data['cycle_id'].unique()[0]
test_cycle_data = test_data[test_data['cycle_id'] == test_cycle].iloc[10:]  # Skip window size
cycle_pred, cycle_true = simple_rul_prediction(
    test_data[test_data['cycle_id'] == test_cycle], window=10
)

axes[1, 1].plot(test_cycle_data['time_step'].values, cycle_true * max_rul,
                label='Actual RUL', color='blue', linewidth=2)
axes[1, 1].plot(test_cycle_data['time_step'].values, cycle_pred * max_rul,
                label='Predicted RUL', color='orange', linewidth=2, linestyle='--')
axes[1, 1].fill_between(test_cycle_data['time_step'].values,
                         cycle_pred * max_rul - 10, cycle_pred * max_rul + 10,
                         alpha=0.2, color='orange', label='Prediction Error Range (±10)')
axes[1, 1].set_xlabel('Time Step')
axes[1, 1].set_ylabel('RUL (Steps)')
axes[1, 1].set_title(f'RUL Time Series Prediction (Cycle {test_cycle})', fontsize=12, fontweight='bold')
axes[1, 1].legend()
axes[1, 1].grid(alpha=0.3)

plt.tight_layout()
plt.savefig('rul_estimation_analysis.png', dpi=300, bbox_inches='tight')
plt.show()

Implementation Points:

• Estimation of future failure timing from time series data (temperature, vibration, pressure, health index)
• In actual implementation, use RNN such as LSTM (Long Short-Term Memory)
• Define health index and learn degradation patterns
• Quantify maintenance plan reliability through prediction error evaluation
• Support decision-making of maintenance personnel through remaining time visualization

🚨 4.3 Anomaly Detection and Early Warning System

Two Approaches to Anomaly Detection

1. Supervised Learning: Learn from past failure cases (requires labeled data)
2. Unsupervised Learning: Detect deviations from normal patterns (no anomaly data required)

💻 Code Example 4.3: Anomaly Detection with Isolation Forest

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.ensemble import IsolationForest
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
import warnings
warnings.filterwarnings('ignore')

plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# Generate normal operation data
np.random.seed(42)
n_normal = 800
n_anomaly = 50

# Normal data (multivariate normal distribution)
mean_normal = [70, 0.25, 2.0, 14, 50]
cov_normal = [[25, 0, 0, 0, 0],
              [0, 0.0025, 0, 0, 0],
              [0, 0, 0.04, 0, 0],
              [0, 0, 0, 4, 0],
              [0, 0, 0, 0, 100]]

normal_data = np.random.multivariate_normal(mean_normal, cov_normal, n_normal)

# Anomaly data (various anomaly patterns)
anomaly_patterns = []

# Pattern 1: High temperature anomaly
high_temp = np.random.multivariate_normal([95, 0.25, 2.0, 14, 50], cov_normal, n_anomaly // 5)
anomaly_patterns.append(high_temp)

# Pattern 2: Vibration anomaly
high_vibration = np.random.multivariate_normal([70, 0.7, 2.0, 14, 50], cov_normal, n_anomaly // 5)
anomaly_patterns.append(high_vibration)

# Pattern 3: Pressure anomaly
high_pressure = np.random.multivariate_normal([70, 0.25, 3.5, 14, 50], cov_normal, n_anomaly // 5)
anomaly_patterns.append(high_pressure)

# Pattern 4: Current anomaly
high_current = np.random.multivariate_normal([70, 0.25, 2.0, 25, 50], cov_normal, n_anomaly // 5)
anomaly_patterns.append(high_current)

# Pattern 5: Combined anomaly
combined = np.random.multivariate_normal([90, 0.6, 2.8, 20, 80], cov_normal, n_anomaly // 5)
anomaly_patterns.append(combined)

anomaly_data = np.vstack(anomaly_patterns)

# Create dataframe
columns = ['temperature', 'vibration', 'pressure', 'current', 'humidity']
df_normal = pd.DataFrame(normal_data, columns=columns)
df_normal['label'] = 0  # Normal

df_anomaly = pd.DataFrame(anomaly_data, columns=columns)
df_anomaly['label'] = 1  # Anomaly

df = pd.concat([df_normal, df_anomaly], ignore_index=True)
df = df.sample(frac=1, random_state=42).reset_index(drop=True)  # Shuffle

# Data standardization
X = df[columns].values
y_true = df['label'].values

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# Isolation Forest model
iso_forest = IsolationForest(
    contamination=0.1,  # Estimated proportion of anomaly data
    random_state=42,
    n_estimators=100,
    max_samples='auto'
)

# Calculate anomaly scores
y_pred = iso_forest.fit_predict(X_scaled)
anomaly_scores = iso_forest.score_samples(X_scaled)

# Convert prediction labels to 0/1 (-1 → 1 (anomaly), 1 → 0 (normal))
y_pred_binary = np.where(y_pred == -1, 1, 0)

# Calculate evaluation metrics
from sklearn.metrics import precision_score, recall_score, f1_score, confusion_matrix

precision = precision_score(y_true, y_pred_binary)
recall = recall_score(y_true, y_pred_binary)
f1 = f1_score(y_true, y_pred_binary)
cm = confusion_matrix(y_true, y_pred_binary)

print("=" * 60)
print("Anomaly Detection System Evaluation (Isolation Forest)")
print("=" * 60)
print(f"\nTotal samples: {len(df)}")
print(f"Normal data: {n_normal} ({n_normal/len(df)*100:.1f}%)")
print(f"Anomaly data: {n_anomaly} ({n_anomaly/len(df)*100:.1f}%)")
print(f"\nDetection Performance:")
print(f"Precision: {precision:.3f}")
print(f"Recall: {recall:.3f}")
print(f"F1 Score: {f1:.3f}")
print(f"\nConfusion Matrix:")
print(f"True Negative: {cm[0,0]}, False Positive: {cm[0,1]}")
print(f"False Negative: {cm[1,0]}, True Positive: {cm[1,1]}")

# Dimensionality reduction to 2D with PCA (for visualization)
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X_scaled)

# Visualization
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Anomaly detection results in PCA space
scatter1 = axes[0, 0].scatter(X_pca[y_pred_binary==0, 0], X_pca[y_pred_binary==0, 1],
                               c='green', s=20, alpha=0.6, label='Classified as Normal')
scatter2 = axes[0, 0].scatter(X_pca[y_pred_binary==1, 0], X_pca[y_pred_binary==1, 1],
                               c='red', s=40, alpha=0.8, marker='X', label='Classified as Anomaly')
axes[0, 0].set_xlabel(f'PC1 (Variance: {pca.explained_variance_ratio_[0]:.1%})')
axes[0, 0].set_ylabel(f'PC2 (Variance: {pca.explained_variance_ratio_[1]:.1%})')
axes[0, 0].set_title('Anomaly Detection Results (PCA Space)', fontsize=12, fontweight='bold')
axes[0, 0].legend()
axes[0, 0].grid(alpha=0.3)

# Anomaly score distribution
axes[0, 1].hist(anomaly_scores[y_true==0], bins=50, alpha=0.6, label='Normal', color='green')
axes[0, 1].hist(anomaly_scores[y_true==1], bins=50, alpha=0.6, label='Anomaly', color='red')
axes[0, 1].set_xlabel('Anomaly Score (Lower = More Anomalous)')
axes[0, 1].set_ylabel('Frequency')
axes[0, 1].set_title('Anomaly Score Distribution', fontsize=12, fontweight='bold')
axes[0, 1].legend()
axes[0, 1].grid(alpha=0.3)

# Confusion matrix heatmap
import seaborn as sns
sns.heatmap(cm, annot=True, fmt='d', cmap='RdYlGn', ax=axes[1, 0],
            xticklabels=['Normal', 'Anomaly'], yticklabels=['Normal', 'Anomaly'])
axes[1, 0].set_title('Confusion Matrix', fontsize=12, fontweight='bold')
axes[1, 0].set_ylabel('True Label')
axes[1, 0].set_xlabel('Predicted Label')

# Time series anomaly detection (simulated real-time monitoring)
time_series_length = 200
time_indices = np.arange(time_series_length)

# Generate time series data (normal → anomaly → normal)
ts_data = []
for i in range(time_series_length):
    if 80 <= i <= 120:  # Anomaly period
        sample = np.random.multivariate_normal([90, 0.6, 2.8, 20, 80], cov_normal)
    else:
        sample = np.random.multivariate_normal(mean_normal, cov_normal)
    ts_data.append(sample)

ts_data = np.array(ts_data)
ts_scaled = scaler.transform(ts_data)
ts_scores = iso_forest.score_samples(ts_scaled)
ts_pred = iso_forest.predict(ts_scaled)

# Calculate anomaly threshold (10th percentile of anomaly scores)
threshold = np.percentile(anomaly_scores, 10)

axes[1, 1].plot(time_indices, ts_scores, color='#11998e', linewidth=1.5, label='Anomaly Score')
axes[1, 1].axhline(y=threshold, color='red', linestyle='--', linewidth=2, label=f'Anomaly Threshold ({threshold:.3f})')
axes[1, 1].fill_between(time_indices, threshold, ts_scores.min(),
                         where=(ts_scores < threshold), alpha=0.3, color='red', label='Anomaly Detected')
axes[1, 1].set_xlabel('Time')
axes[1, 1].set_ylabel('Anomaly Score')
axes[1, 1].set_title('Real-time Anomaly Detection (Simulated)', fontsize=12, fontweight='bold')
axes[1, 1].legend()
axes[1, 1].grid(alpha=0.3)

plt.tight_layout()
plt.savefig('anomaly_detection_analysis.png', dpi=300, bbox_inches='tight')
plt.show()

# Generate warning messages when anomalies are detected
detected_anomalies = np.where(ts_pred == -1)[0]
if len(detected_anomalies) > 0:
    print(f"\n⚠️ Warning: {len(detected_anomalies)} anomalies detected")
    print(f"Anomaly detection times: {detected_anomalies[:10]}..." if len(detected_anomalies) > 10 else f"Anomaly detection times: {detected_anomalies}")

Implementation Points:

• Isolation Forest can be trained with normal data only (unsupervised learning)
• Automatically detect anomaly patterns from multi-dimensional sensor data
• Easy integration into real-time monitoring systems
• Warning sensitivity can be tuned through threshold adjustment
• Achieve preventive maintenance before failure through early warning

🔍 4.4 Root Cause Analysis (RCA)

Analysis Flow When Trouble Occurs

1. Symptom Identification: What happened (quality defect, equipment shutdown, etc.)
2. Direct Cause Investigation: Confirm immediately preceding changes/events
3. Root Cause Identification: Why it happened (5 Whys)
4. Countermeasure Planning: Consider recurrence prevention measures

Data-Driven RCA

To identify factors that contributed to quality anomalies or equipment failures from large amounts of process data, statistical methods and machine learning are effective. Visualize parameters with high impact using correlation analysis, decision trees, SHAP values, etc.

💻 Code Example 4.4: Root Cause Analysis with Decision Trees and SHAP Values

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.tree import DecisionTreeClassifier, plot_tree
from sklearn.model_selection import train_test_split
import warnings
warnings.filterwarnings('ignore')

plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# Generate quality defect data (simulation data with clear root causes)
np.random.seed(42)
n_samples = 500

# Normal products (60%)
n_normal = int(n_samples * 0.6)
temp_normal = np.random.normal(75, 3, n_normal)
ph_normal = np.random.normal(6.5, 0.2, n_normal)
humidity_normal = np.random.normal(50, 5, n_normal)
mixing_time_normal = np.random.normal(120, 10, n_normal)
additive_amount_normal = np.random.normal(2.0, 0.1, n_normal)

# Defective pattern 1: High temperature cause (20%)
n_defect1 = int(n_samples * 0.2)
temp_defect1 = np.random.normal(90, 5, n_defect1)
ph_defect1 = np.random.normal(6.5, 0.2, n_defect1)
humidity_defect1 = np.random.normal(50, 5, n_defect1)
mixing_time_defect1 = np.random.normal(120, 10, n_defect1)
additive_amount_defect1 = np.random.normal(2.0, 0.1, n_defect1)

# Defective pattern 2: pH anomaly cause (10%)
n_defect2 = int(n_samples * 0.1)
temp_defect2 = np.random.normal(75, 3, n_defect2)
ph_defect2 = np.random.normal(7.5, 0.3, n_defect2)
humidity_defect2 = np.random.normal(50, 5, n_defect2)
mixing_time_defect2 = np.random.normal(120, 10, n_defect2)
additive_amount_defect2 = np.random.normal(2.0, 0.1, n_defect2)

# Defective pattern 3: Insufficient additive amount cause (10%)
n_defect3 = n_samples - n_normal - n_defect1 - n_defect2
temp_defect3 = np.random.normal(75, 3, n_defect3)
ph_defect3 = np.random.normal(6.5, 0.2, n_defect3)
humidity_defect3 = np.random.normal(50, 5, n_defect3)
mixing_time_defect3 = np.random.normal(120, 10, n_defect3)
additive_amount_defect3 = np.random.normal(1.5, 0.15, n_defect3)

# Create dataframe
data = pd.DataFrame({
    'temperature': np.concatenate([temp_normal, temp_defect1, temp_defect2, temp_defect3]),
    'pH': np.concatenate([ph_normal, ph_defect1, ph_defect2, ph_defect3]),
    'humidity': np.concatenate([humidity_normal, humidity_defect1, humidity_defect2, humidity_defect3]),
    'mixing_time': np.concatenate([mixing_time_normal, mixing_time_defect1, mixing_time_defect2, mixing_time_defect3]),
    'additive_amount': np.concatenate([additive_amount_normal, additive_amount_defect1, additive_amount_defect2, additive_amount_defect3]),
    'quality': np.concatenate([
        np.zeros(n_normal),
        np.ones(n_defect1),
        np.ones(n_defect2),
        np.ones(n_defect3)
    ])
})

# Shuffle data
data = data.sample(frac=1, random_state=42).reset_index(drop=True)

# Split training and test data
X = data.drop('quality', axis=1)
y = data['quality']
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42, stratify=y)

# Train decision tree model (limit depth to prioritize interpretability)
dt_model = DecisionTreeClassifier(
    max_depth=4,
    min_samples_split=20,
    min_samples_leaf=10,
    random_state=42
)
dt_model.fit(X_train, y_train)

# Prediction and evaluation
y_pred = dt_model.predict(X_test)
from sklearn.metrics import classification_report, accuracy_score

accuracy = accuracy_score(y_test, y_pred)

print("=" * 60)
print("Root Cause Analysis (Decision Tree)")
print("=" * 60)
print(f"\nModel Accuracy: {accuracy:.3f}")
print("\nClassification Report:")
print(classification_report(y_test, y_pred, target_names=['Normal', 'Defective']))

# Feature importance
feature_importance = pd.DataFrame({
    'feature': X.columns,
    'importance': dt_model.feature_importances_
}).sort_values('importance', ascending=False)

print("\nFeature Importance (Contribution to Quality Defects):")
print(feature_importance.to_string(index=False))

# Visualization
fig = plt.figure(figsize=(16, 10))

# Decision tree visualization
ax1 = plt.subplot(2, 2, (1, 2))
plot_tree(dt_model, feature_names=X.columns, class_names=['Normal', 'Defective'],
          filled=True, rounded=True, fontsize=9, ax=ax1)
ax1.set_title('Root Cause Analysis by Decision Tree', fontsize=14, fontweight='bold')

# Feature importance
ax2 = plt.subplot(2, 2, 3)
bars = ax2.barh(feature_importance['feature'], feature_importance['importance'], color='#38ef7d')
ax2.set_xlabel('Importance Score')
ax2.set_title('Contribution to Quality Defects', fontsize=12, fontweight='bold')
ax2.grid(axis='x', alpha=0.3)

# Important feature distribution comparison
ax3 = plt.subplot(2, 2, 4)
top_feature = feature_importance.iloc[0]['feature']

normal_values = data[data['quality'] == 0][top_feature]
defect_values = data[data['quality'] == 1][top_feature]

ax3.hist(normal_values, bins=30, alpha=0.6, label='Normal', color='green')
ax3.hist(defect_values, bins=30, alpha=0.6, label='Defective', color='red')
ax3.set_xlabel(f'{top_feature}')
ax3.set_ylabel('Frequency')
ax3.set_title(f'Distribution Comparison of Most Important Factor: {top_feature}', fontsize=12, fontweight='bold')
ax3.legend()
ax3.grid(alpha=0.3)

plt.tight_layout()
plt.savefig('root_cause_analysis.png', dpi=300, bbox_inches='tight')
plt.show()

# Generate RCA report
print("\n" + "=" * 60)
print("Root Cause Analysis Report")
print("=" * 60)

# Extract decision tree rules (top nodes only)
from sklearn.tree import _tree

def extract_rules(tree, feature_names, node=0, depth=0, max_depth=2):
    """Extract rules from decision tree"""
    if depth > max_depth:
        return

    feature = tree.feature[node]
    threshold = tree.threshold[node]

    if feature != _tree.TREE_UNDEFINED:
        name = feature_names[feature]
        print(f"{'  ' * depth}├─ If {name} <= {threshold:.2f}:")
        extract_rules(tree, feature_names, tree.children_left[node], depth + 1, max_depth)
        print(f"{'  ' * depth}└─ If {name} > {threshold:.2f}:")
        extract_rules(tree, feature_names, tree.children_right[node], depth + 1, max_depth)
    else:
        value = tree.value[node][0]
        total = value.sum()
        defect_rate = value[1] / total if total > 0 else 0
        print(f"{'  ' * depth}   → Defect Rate: {defect_rate:.1%} (Samples: {int(total)})")

print("\nKey Decision Rules:")
extract_rules(dt_model.tree_, X.columns)

print("\nRecommended Countermeasures:")
if feature_importance.iloc[0]['feature'] == 'temperature':
    print("- Strengthen temperature management: Control within 75±5°C range")
    print("- Inspect and improve cooling system")
elif feature_importance.iloc[0]['feature'] == 'pH':
    print("- Review pH adjustment process")
    print("- Increase pH meter calibration frequency")
elif feature_importance.iloc[0]['feature'] == 'additive_amount':
    print("- Stricter management of additive dosage")
    print("- Consider introducing automated dosing system")

Implementation Points:

• Visualize decision rules for defect occurrence using decision trees
• Identify parameters with high impact through feature importance
• Individual factors for each sample can be analyzed using SHAP values (SHAP library required in implementation example)
• Extract human-understandable rules to support countermeasure planning
• Data-driven RCA enables analysis not relying solely on empirical rules

📚 Summary

In this chapter, we learned about AI-powered predictive maintenance and troubleshooting technologies.

Key Points

• Building equipment failure prediction models using machine learning
• Optimal maintenance planning through Remaining Useful Life (RUL) estimation
• Anomaly detection using unsupervised learning (Isolation Forest)
• Root cause analysis using decision trees and SHAP values
• Improving maintenance efficiency through data-driven approaches

Practical Applications

• Continuous collection of sensor data and periodic model retraining
• Planned maintenance scheduling through early warning systems
• Equipment improvement and process optimization utilizing RCA results
• Optimizing the balance between maintenance costs and equipment uptime