Chapter 3: Anomaly Detection and Fault Diagnosis

AI Terakoya Top›Process Informatics›Process Data Analysis›Chapter 3

🌐 EN | 🇯🇵 JP | Last sync: 2025-11-16

3.1 Introduction

To maintain the safety and productivity of chemical plants, early detection of anomalies and rapid diagnosis of failure causes are essential. Technologies that automatically detect anomalous patterns from sensor data and identify fault locations play a crucial role in advancing plant operations.

In this chapter, we will implement 8 practical anomaly detection and fault diagnosis techniques, ranging from statistical methods to machine learning algorithms and deep learning. You will understand the characteristics of each method and become able to select the appropriate approach according to process characteristics.

📊 What You Will Learn in This Chapter

• Statistical anomaly detection (Z-score, Modified Z-score)
• Machine learning-based anomaly detection (Isolation Forest, One-Class SVM)
• Deep learning-based anomaly detection (Autoencoder, LSTM)
• Fault classification and ensemble methods
• Root cause analysis (Granger causality test)

3.2 Statistical Anomaly Detection

Statistical methods are suitable for real-time monitoring because of their high interpretability and low computational cost. By combining Z-score-based methods assuming normal distribution with the robust Modified Z-score, we achieve high-precision detection.

Example 1: Statistical Anomaly Detection (Z-score & Modified Z-score)

We implement anomaly detection combining standard Z-score and the robust Modified Z-score.

# ===================================
# Example 1: Statistical Anomaly Detection
# ===================================
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats

# Generate process data (reactor temperature, including anomalies)
np.random.seed(42)
n_samples = 500
normal_data = np.random.normal(350, 2, n_samples)

# Add anomalous data (3 types)
anomaly_indices = [50, 150, 250, 350, 450]
normal_data[anomaly_indices] = [360, 340, 365, 335, 362]  # Anomalous values

df = pd.DataFrame({'temperature': normal_data})

def detect_anomalies_zscore(data, threshold=3.0):
    """Anomaly detection using standard Z-score method"""
    mean = data.mean()
    std = data.std()
    z_scores = np.abs((data - mean) / std)
    return z_scores > threshold, z_scores

def detect_anomalies_modified_zscore(data, threshold=3.5):
    """Modified Z-score method (median-based, robust to outliers)"""
    median = data.median()
    mad = np.median(np.abs(data - median))
    modified_z_scores = 0.6745 * (data - median) / mad
    return np.abs(modified_z_scores) > threshold, modified_z_scores

# Apply both methods
is_anomaly_z, z_scores = detect_anomalies_zscore(df['temperature'])
is_anomaly_mod, mod_z_scores = detect_anomalies_modified_zscore(df['temperature'])

# Compare results
print("Anomaly Detection Results:")
print(f"Z-score method: Detected {is_anomaly_z.sum()} anomalies")
print(f"Modified Z-score method: Detected {is_anomaly_mod.sum()} anomalies")

# Evaluate detection accuracy
true_anomalies = np.zeros(n_samples, dtype=bool)
true_anomalies[anomaly_indices] = True

tp_z = np.sum(is_anomaly_z & true_anomalies)
fp_z = np.sum(is_anomaly_z & ~true_anomalies)
precision_z = tp_z / is_anomaly_z.sum() if is_anomaly_z.sum() > 0 else 0
recall_z = tp_z / true_anomalies.sum()

tp_mod = np.sum(is_anomaly_mod & true_anomalies)
fp_mod = np.sum(is_anomaly_mod & ~true_anomalies)
precision_mod = tp_mod / is_anomaly_mod.sum() if is_anomaly_mod.sum() > 0 else 0
recall_mod = tp_mod / true_anomalies.sum()

print(f"\nZ-score method: Precision={precision_z:.2f}, Recall={recall_z:.2f}")
print(f"Modified Z-score method: Precision={precision_mod:.2f}, Recall={recall_mod:.2f}")

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

axes[0].plot(df['temperature'], 'b-', alpha=0.6, label='Data')
axes[0].scatter(np.where(is_anomaly_z)[0], df.loc[is_anomaly_z, 'temperature'],
                color='red', s=100, zorder=5, label=f'Anomalies (Z-score)')
axes[0].axhline(df['temperature'].mean() + 3*df['temperature'].std(),
                color='r', linestyle='--', alpha=0.5, label='±3σ threshold')
axes[0].axhline(df['temperature'].mean() - 3*df['temperature'].std(),
                color='r', linestyle='--', alpha=0.5)
axes[0].set_ylabel('Temperature (°C)')
axes[0].set_title('Z-score Method')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

axes[1].plot(df['temperature'], 'b-', alpha=0.6, label='Data')
axes[1].scatter(np.where(is_anomaly_mod)[0], df.loc[is_anomaly_mod, 'temperature'],
                color='red', s=100, zorder=5, label='Anomalies (Modified Z-score)')
axes[1].set_xlabel('Sample')
axes[1].set_ylabel('Temperature (°C)')
axes[1].set_title('Modified Z-score Method (Robust to outliers)')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('statistical_anomaly_detection.png', dpi=300)
print("\nResult: Modified Z-score is less affected by outliers and more accurate")

3.3 Machine Learning-Based Anomaly Detection

Machine learning algorithms excel at detecting high-dimensional data and complex nonlinear patterns. Unsupervised learning that can identify anomalies from normal data only is a major practical advantage.

Example 2: Anomaly Detection with Isolation Forest

We implement an efficient anomaly detection algorithm based on random forests.

# ===================================
# Example 2: Isolation Forest
# ===================================
from sklearn.ensemble import IsolationForest

# Generate multivariate process data
np.random.seed(42)
n_normal = 450
n_anomaly = 50

# Normal data (correlated temperature, pressure, and flow)
normal_temp = np.random.normal(350, 3, n_normal)
normal_pressure = 5 + 0.01 * normal_temp + np.random.normal(0, 0.2, n_normal)
normal_flow = 100 + 0.5 * normal_temp + np.random.normal(0, 5, n_normal)

# Anomalous data (correlation breaks down)
anomaly_temp = np.random.uniform(340, 370, n_anomaly)
anomaly_pressure = np.random.uniform(4, 7, n_anomaly)
anomaly_flow = np.random.uniform(80, 150, n_anomaly)

# Combine data
X = np.column_stack([
    np.concatenate([normal_temp, anomaly_temp]),
    np.concatenate([normal_pressure, anomaly_pressure]),
    np.concatenate([normal_flow, anomaly_flow])
])
y_true = np.concatenate([np.zeros(n_normal), np.ones(n_anomaly)])

df_process = pd.DataFrame(X, columns=['temperature', 'pressure', 'flow'])
df_process['is_anomaly_true'] = y_true

# Train Isolation Forest
iso_forest = IsolationForest(
    contamination=0.1,  # Expected anomaly proportion
    random_state=42,
    n_estimators=100
)
predictions = iso_forest.fit_predict(X)
anomaly_scores = iso_forest.score_samples(X)

df_process['anomaly_score'] = anomaly_scores
df_process['is_anomaly_pred'] = predictions == -1

# Evaluation
from sklearn.metrics import classification_report, confusion_matrix

print("Isolation Forest Results:")
print(classification_report(y_true, predictions == -1, target_names=['Normal', 'Anomaly']))
print("\nConfusion Matrix:")
print(confusion_matrix(y_true, predictions == -1))

# Visualization (3D scatter plot)
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure(figsize=(14, 6))
ax1 = fig.add_subplot(121, projection='3d')
ax2 = fig.add_subplot(122, projection='3d')

# True anomalies
scatter1 = ax1.scatter(df_process['temperature'], df_process['pressure'],
                       df_process['flow'], c=df_process['is_anomaly_true'],
                       cmap='RdYlGn_r', s=30, alpha=0.6)
ax1.set_xlabel('Temperature (°C)')
ax1.set_ylabel('Pressure (MPa)')
ax1.set_zlabel('Flow (m³/h)')
ax1.set_title('True Anomaly Labels')
plt.colorbar(scatter1, ax=ax1)

# Detection results
scatter2 = ax2.scatter(df_process['temperature'], df_process['pressure'],
                       df_process['flow'], c=df_process['anomaly_score'],
                       cmap='RdYlGn', s=30, alpha=0.6)
ax2.set_xlabel('Temperature (°C)')
ax2.set_ylabel('Pressure (MPa)')
ax2.set_zlabel('Flow (m³/h)')
ax2.set_title('Isolation Forest Anomaly Score')
plt.colorbar(scatter2, ax=ax2)

plt.tight_layout()
plt.savefig('isolation_forest.png', dpi=300)
print("\nResult: Efficiently detects anomalous patterns in high-dimensional space")

Example 3: Novelty Detection with One-Class SVM

We learn the boundary of normal data and detect unknown anomalous patterns.

# ===================================
# Example 3: One-Class SVM
# ===================================
from sklearn.svm import OneClassSVM
from sklearn.preprocessing import StandardScaler

# Standardize data (important for SVM)
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# Train One-Class SVM
oc_svm = OneClassSVM(
    kernel='rbf',
    gamma='auto',
    nu=0.1  # Upper bound on anomaly proportion
)
predictions_svm = oc_svm.fit_predict(X_scaled)
decision_scores = oc_svm.decision_function(X_scaled)

df_process['decision_score'] = decision_scores
df_process['is_anomaly_svm'] = predictions_svm == -1

# Evaluation
print("\nOne-Class SVM Results:")
print(classification_report(y_true, predictions_svm == -1, target_names=['Normal', 'Anomaly']))

# Visualize decision boundary (2D projection)
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Temperature-Pressure plane
xx, yy = np.meshgrid(np.linspace(X_scaled[:, 0].min()-1, X_scaled[:, 0].max()+1, 100),
                     np.linspace(X_scaled[:, 1].min()-1, X_scaled[:, 1].max()+1, 100))
Z = oc_svm.decision_function(np.c_[xx.ravel(), yy.ravel(), np.zeros(xx.ravel().shape[0])])
Z = Z.reshape(xx.shape)

axes[0].contourf(xx, yy, Z, levels=np.linspace(Z.min(), 0, 7), cmap='RdYlGn', alpha=0.3)
axes[0].contour(xx, yy, Z, levels=[0], linewidths=2, colors='red')
axes[0].scatter(X_scaled[y_true==0, 0], X_scaled[y_true==0, 1],
                c='green', s=20, alpha=0.6, label='Normal')
axes[0].scatter(X_scaled[y_true==1, 0], X_scaled[y_true==1, 1],
                c='red', s=50, alpha=0.8, marker='x', label='Anomaly')
axes[0].set_xlabel('Temperature (standardized)')
axes[0].set_ylabel('Pressure (standardized)')
axes[0].set_title('One-Class SVM Decision Boundary')
axes[0].legend()

# Distribution of decision scores
axes[1].hist(decision_scores[y_true==0], bins=30, alpha=0.6, label='Normal', color='green')
axes[1].hist(decision_scores[y_true==1], bins=30, alpha=0.6, label='Anomaly', color='red')
axes[1].axvline(0, color='black', linestyle='--', label='Decision threshold')
axes[1].set_xlabel('Decision Score')
axes[1].set_ylabel('Frequency')
axes[1].set_title('Decision Score Distribution')
axes[1].legend()

plt.tight_layout()
plt.savefig('one_class_svm.png', dpi=300)
print("\nResult: Rigorous modeling of normal data boundary, capable of handling novel anomalies")

3.4 Deep Learning-Based Anomaly Detection

Deep learning demonstrates its power in detecting complex spatiotemporal patterns and anomalies in high-dimensional data. We implement two approaches using Autoencoder and LSTM.

Example 4: Reconstruction Error-Based Anomaly Detection with Autoencoder

We compress and restore normal data, detecting anomalies from reconstruction errors.

# ===================================
# Example 4: Autoencoder Anomaly Detection
# ===================================
import tensorflow as tf
from tensorflow import keras
from sklearn.model_selection import train_test_split

# Train only on normal data (unsupervised learning)
X_train = X[y_true == 0]
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X)

# Build Autoencoder model
input_dim = X_train_scaled.shape[1]
encoding_dim = 2  # Compressed dimension

autoencoder = keras.Sequential([
    keras.layers.Dense(8, activation='relu', input_shape=(input_dim,)),
    keras.layers.Dense(encoding_dim, activation='relu', name='encoder'),
    keras.layers.Dense(8, activation='relu'),
    keras.layers.Dense(input_dim, activation='linear', name='decoder')
])

autoencoder.compile(optimizer='adam', loss='mse')

# Training
history = autoencoder.fit(
    X_train_scaled, X_train_scaled,
    epochs=50,
    batch_size=32,
    validation_split=0.2,
    verbose=0
)

# Calculate reconstruction errors
X_reconstructed = autoencoder.predict(X_test_scaled)
reconstruction_errors = np.mean(np.square(X_test_scaled - X_reconstructed), axis=1)

# Anomaly detection threshold (95th percentile of normal data)
threshold = np.percentile(reconstruction_errors[y_true==0], 95)
predictions_ae = reconstruction_errors > threshold

# Evaluation
print("\nAutoencoder Results:")
print(f"Threshold: {threshold:.4f}")
print(classification_report(y_true, predictions_ae, target_names=['Normal', 'Anomaly']))

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

# Learning curve
axes[0, 0].plot(history.history['loss'], label='Training loss')
axes[0, 0].plot(history.history['val_loss'], label='Validation loss')
axes[0, 0].set_xlabel('Epoch')
axes[0, 0].set_ylabel('MSE')
axes[0, 0].set_title('Autoencoder Learning Curve')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# Reconstruction error distribution
axes[0, 1].hist(reconstruction_errors[y_true==0], bins=30, alpha=0.6,
                label='Normal', color='green')
axes[0, 1].hist(reconstruction_errors[y_true==1], bins=30, alpha=0.6,
                label='Anomaly', color='red')
axes[0, 1].axvline(threshold, color='black', linestyle='--', label='Threshold')
axes[0, 1].set_xlabel('Reconstruction Error')
axes[0, 1].set_ylabel('Frequency')
axes[0, 1].set_title('Reconstruction Error Distribution')
axes[0, 1].legend()

# Latent space (2D encoding)
encoder_model = keras.Model(autoencoder.input,
                            autoencoder.get_layer('encoder').output)
encoded = encoder_model.predict(X_test_scaled)

axes[1, 0].scatter(encoded[y_true==0, 0], encoded[y_true==0, 1],
                   c='green', s=20, alpha=0.6, label='Normal')
axes[1, 0].scatter(encoded[y_true==1, 0], encoded[y_true==1, 1],
                   c='red', s=50, alpha=0.8, marker='x', label='Anomaly')
axes[1, 0].set_xlabel('Latent Dimension 1')
axes[1, 0].set_ylabel('Latent Dimension 2')
axes[1, 0].set_title('Latent Space Representation')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)

# ROC curve
from sklearn.metrics import roc_curve, auc
fpr, tpr, _ = roc_curve(y_true, reconstruction_errors)
roc_auc = auc(fpr, tpr)

axes[1, 1].plot(fpr, tpr, label=f'ROC (AUC = {roc_auc:.2f})')
axes[1, 1].plot([0, 1], [0, 1], 'k--', label='Random')
axes[1, 1].set_xlabel('False Positive Rate')
axes[1, 1].set_ylabel('True Positive Rate')
axes[1, 1].set_title('ROC Curve')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('autoencoder_anomaly.png', dpi=300)
print(f"\nResult: AUC={roc_auc:.3f}, learns nonlinear patterns in high-dimensional data")

Example 5: Time Series Anomaly Detection with LSTM

We learn temporal dependencies in time series and detect anomalies from prediction errors.

# ===================================
# Example 5: LSTM Time Series Anomaly Detection
# ===================================

# Generate time series data (seasonality + trend + anomalies)
np.random.seed(42)
n_timesteps = 1000
t = np.arange(n_timesteps)

# Normal pattern
normal_series = 100 + 0.01*t + 10*np.sin(2*np.pi*t/50) + np.random.normal(0, 1, n_timesteps)

# Inject anomalies
anomaly_ranges = [(200, 220), (500, 530), (800, 815)]
for start, end in anomaly_ranges:
    normal_series[start:end] += np.random.uniform(15, 25, end-start)

# Create sequence data
def create_sequences(data, seq_length=50):
    X, y = [], []
    for i in range(len(data) - seq_length):
        X.append(data[i:i+seq_length])
        y.append(data[i+seq_length])
    return np.array(X), np.array(y)

seq_length = 50
X_seq, y_seq = create_sequences(normal_series, seq_length)
X_seq = X_seq.reshape((X_seq.shape[0], X_seq.shape[1], 1))

# Train/test split (train on parts without anomalies)
split_idx = 150 - seq_length
X_train_lstm = X_seq[:split_idx]
y_train_lstm = y_seq[:split_idx]
X_test_lstm = X_seq
y_test_lstm = y_seq

# Build LSTM model
lstm_model = keras.Sequential([
    keras.layers.LSTM(32, activation='relu', input_shape=(seq_length, 1)),
    keras.layers.Dense(16, activation='relu'),
    keras.layers.Dense(1)
])

lstm_model.compile(optimizer='adam', loss='mse')

# Training
history_lstm = lstm_model.fit(
    X_train_lstm, y_train_lstm,
    epochs=30,
    batch_size=32,
    validation_split=0.2,
    verbose=0
)

# Prediction and error calculation
predictions_lstm = lstm_model.predict(X_test_lstm).flatten()
prediction_errors = np.abs(y_test_lstm - predictions_lstm)

# Anomaly threshold (95th percentile of training data)
train_errors = prediction_errors[:split_idx]
threshold_lstm = np.percentile(train_errors, 95)
anomalies_lstm = prediction_errors > threshold_lstm

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

# Original data
axes[0].plot(normal_series, 'b-', alpha=0.7, label='Actual values')
for start, end in anomaly_ranges:
    axes[0].axvspan(start, end, color='red', alpha=0.2)
axes[0].set_ylabel('Process Value')
axes[0].set_title('Time Series Data (Anomaly periods shown in red)')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Prediction vs actual
axes[1].plot(y_test_lstm, 'b-', alpha=0.5, label='Actual values')
axes[1].plot(predictions_lstm, 'g--', alpha=0.7, label='Predicted values')
axes[1].set_ylabel('Process Value')
axes[1].set_title('LSTM Prediction Results')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

# Prediction error and anomaly detection
axes[2].plot(prediction_errors, 'k-', alpha=0.6, label='Prediction error')
axes[2].axhline(threshold_lstm, color='red', linestyle='--', label=f'Threshold={threshold_lstm:.2f}')
axes[2].fill_between(range(len(anomalies_lstm)), 0, prediction_errors,
                     where=anomalies_lstm, color='red', alpha=0.3, label='Detected anomalies')
axes[2].set_xlabel('Time')
axes[2].set_ylabel('Prediction Error')
axes[2].set_title('LSTM Anomaly Detection Results')
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('lstm_anomaly.png', dpi=300)
print("\nResult: Learns time series patterns, detects anomalous periods with high accuracy")

3.5 Fault Classification and Ensemble Methods

Not only detecting anomalies but also classifying fault types enables rapid determination of appropriate countermeasures. We build a robust classification system using Random Forest and ensemble voting.

Example 6: Fault Classification with Random Forest

We classify multiple fault types using machine learning.

# ===================================
# Example 6: Fault Classification (Random Forest)
# ===================================
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import classification_report, confusion_matrix
import seaborn as sns

# Generate fault data (4 classes: normal, sensor fault, process fault, control fault)
np.random.seed(42)
n_per_class = 200

# Class 0: Normal
normal = np.random.multivariate_normal([350, 5, 100], [[4, 0.1, 5], [0.1, 0.04, 0.5], [5, 0.5, 25]], n_per_class)

# Class 1: Sensor fault (increased noise)
sensor_fault = np.random.multivariate_normal([350, 5, 100], [[16, 0.1, 5], [0.1, 0.16, 0.5], [5, 0.5, 100]], n_per_class)

# Class 2: Process fault (temperature rise)
process_fault = np.random.multivariate_normal([365, 5.5, 120], [[4, 0.2, 8], [0.2, 0.04, 0.6], [8, 0.6, 25]], n_per_class)

# Class 3: Control fault (changed variation pattern)
control_fault = np.random.multivariate_normal([350, 4.5, 90], [[9, -0.3, -10], [-0.3, 0.09, 2], [-10, 2, 64]], n_per_class)

# Combine data
X_fault = np.vstack([normal, sensor_fault, process_fault, control_fault])
y_fault = np.array([0]*n_per_class + [1]*n_per_class + [2]*n_per_class + [3]*n_per_class)

# Train/test split
X_train_fault, X_test_fault, y_train_fault, y_test_fault = train_test_split(
    X_fault, y_fault, test_size=0.3, random_state=42, stratify=y_fault
)

# Train Random Forest
rf_classifier = RandomForestClassifier(n_estimators=100, random_state=42)
rf_classifier.fit(X_train_fault, y_train_fault)

# Prediction
y_pred_rf = rf_classifier.predict(X_test_fault)
y_pred_proba = rf_classifier.predict_proba(X_test_fault)

# Evaluation
fault_names = ['Normal', 'Sensor Fault', 'Process Fault', 'Control Fault']
print("\nRandom Forest Fault Classification Results:")
print(classification_report(y_test_fault, y_pred_rf, target_names=fault_names))

# Confusion matrix
cm = confusion_matrix(y_test_fault, y_pred_rf)
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

sns.heatmap(cm, annot=True, fmt='d', cmap='Blues', xticklabels=fault_names,
            yticklabels=fault_names, ax=axes[0])
axes[0].set_xlabel('Predicted')
axes[0].set_ylabel('True')
axes[0].set_title('Confusion Matrix')

# Feature importance
importances = rf_classifier.feature_importances_
feature_names = ['Temperature', 'Pressure', 'Flow']
axes[1].barh(feature_names, importances, color='teal')
axes[1].set_xlabel('Importance')
axes[1].set_title('Feature Importance')

plt.tight_layout()
plt.savefig('fault_classification.png', dpi=300)
print("\nResult: Classifies multiple fault types with high accuracy, provides diagnostic rationale from feature importance")

Example 7: Ensemble Anomaly Detection (Voting)

We combine multiple anomaly detection methods to improve robustness.

# ===================================
# Example 7: Ensemble Anomaly Detection
# ===================================

# Train multiple detectors (using previously defined data)
# 1. Isolation Forest
iso_pred = (iso_forest.predict(X) == -1).astype(int)

# 2. One-Class SVM
svm_pred = (oc_svm.predict(X_scaled) == -1).astype(int)

# 3. Autoencoder (reconstruction error)
ae_pred = (reconstruction_errors > threshold).astype(int)

# Ensemble voting
ensemble_votes = np.column_stack([iso_pred, svm_pred, ae_pred])
ensemble_pred = (ensemble_votes.sum(axis=1) >= 2).astype(int)  # Majority vote

# Evaluation
from sklearn.metrics import accuracy_score, f1_score

methods = {
    'Isolation Forest': iso_pred,
    'One-Class SVM': svm_pred,
    'Autoencoder': ae_pred,
    'Ensemble (Voting)': ensemble_pred
}

print("\nEnsemble Anomaly Detection Results:")
print(f"{'Method':<25} {'Accuracy':<10} {'F1-score':<10}")
print("-" * 45)
for name, pred in methods.items():
    acc = accuracy_score(y_true, pred)
    f1 = f1_score(y_true, pred)
    print(f"{name:<25} {acc:.4f}     {f1:.4f}")

# Agreement analysis between individual methods
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Distribution of votes
vote_counts = ensemble_votes.sum(axis=1)
axes[0].hist([vote_counts[y_true==0], vote_counts[y_true==1]],
             bins=4, label=['Normal', 'Anomaly'], color=['green', 'red'], alpha=0.6)
axes[0].set_xlabel('Number of Detector Votes')
axes[0].set_ylabel('Frequency')
axes[0].set_title('Ensemble Vote Distribution')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Correlation between methods
method_matrix = np.column_stack([iso_pred, svm_pred, ae_pred, ensemble_pred])
correlation = np.corrcoef(method_matrix.T)
sns.heatmap(correlation, annot=True, fmt='.2f', cmap='coolwarm',
            xticklabels=['Iso Forest', 'SVM', 'AE', 'Ensemble'],
            yticklabels=['Iso Forest', 'SVM', 'AE', 'Ensemble'],
            ax=axes[1], vmin=-1, vmax=1)
axes[1].set_title('Correlation Between Detection Methods')

plt.tight_layout()
plt.savefig('ensemble_anomaly.png', dpi=300)
print("\nResult: Ensemble complements weaknesses of individual methods and improves accuracy")

3.6 Root Cause Analysis

When anomalies are detected, it is important to identify their causes. The Granger causality test estimates causal relationships between variables from time series data, narrowing down root causes.

Example 8: Root Cause Analysis with Granger Causality Test

We analyze causal relationships among multiple process variables to identify fault causes.

# ===================================
# Example 8: Granger Causality Analysis
# ===================================
from statsmodels.tsa.stattools import grangercausalitytests

# Generate process data with causal relationships
np.random.seed(42)
n = 500

# Cause variable (catalyst temperature)
catalyst_temp = np.zeros(n)
catalyst_temp[0] = 300
for i in range(1, n):
    catalyst_temp[i] = 0.95 * catalyst_temp[i-1] + 300 * 0.05 + np.random.normal(0, 1)

# Catalyst temperature affects reaction rate (with delay)
reaction_rate = np.zeros(n)
for i in range(3, n):
    reaction_rate[i] = 50 + 0.3 * catalyst_temp[i-2] + np.random.normal(0, 2)

# Reaction rate affects product yield
product_yield = np.zeros(n)
for i in range(2, n):
    product_yield[i] = 70 + 0.5 * reaction_rate[i-1] + np.random.normal(0, 1.5)

# Create dataframe
df_causal = pd.DataFrame({
    'catalyst_temp': catalyst_temp,
    'reaction_rate': reaction_rate,
    'product_yield': product_yield
})

# Granger causality test function
def granger_causality_matrix(data, variables, max_lag=5):
    """Test Granger causality between variables and display as matrix"""
    df_results = pd.DataFrame(np.zeros((len(variables), len(variables))),
                              columns=variables, index=variables)

    for c in variables:
        for r in variables:
            if c != r:
                test_result = grangercausalitytests(
                    data[[r, c]], max_lag, verbose=False
                )
                # Use minimum p-value across lags
                p_values = [test_result[lag][0]['ssr_ftest'][1] for lag in range(1, max_lag+1)]
                min_p = np.min(p_values)
                df_results.loc[r, c] = min_p

    return df_results

# Perform causality test
variables = ['catalyst_temp', 'reaction_rate', 'product_yield']
causality_matrix = granger_causality_matrix(df_causal[50:], variables, max_lag=5)

print("\nGranger Causality Test Results (p-values):")
print("Column → Row indicates causality (significant if p < 0.05)")
print(causality_matrix.round(4))

# Extract significant causal relationships
print("\nSignificant Causal Relationships (p < 0.05):")
for cause in variables:
    for effect in variables:
        p_value = causality_matrix.loc[effect, cause]
        if p_value < 0.05 and cause != effect:
            print(f"  {cause} → {effect} (p={p_value:.4f})")

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

# Time series plots
axes[0, 0].plot(df_causal['catalyst_temp'], label='Catalyst Temperature')
axes[0, 0].set_ylabel('Temperature (°C)')
axes[0, 0].set_title('Catalyst Temperature (Cause Variable)')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

axes[0, 1].plot(df_causal['reaction_rate'], label='Reaction Rate', color='orange')
axes[0, 1].set_ylabel('Rate (mol/s)')
axes[0, 1].set_title('Reaction Rate (Intermediate Variable)')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

axes[1, 0].plot(df_causal['product_yield'], label='Product Yield', color='green')
axes[1, 0].set_xlabel('Time')
axes[1, 0].set_ylabel('Yield (%)')
axes[1, 0].set_title('Product Yield (Effect Variable)')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)

# Causality matrix heatmap
significance_matrix = (causality_matrix < 0.05).astype(int)
sns.heatmap(significance_matrix, annot=causality_matrix.round(3), fmt='',
            cmap='RdYlGn_r', cbar_kws={'label': 'Significance'},
            xticklabels=variables, yticklabels=variables,
            ax=axes[1, 1], vmin=0, vmax=1)
axes[1, 1].set_title('Causal Relationship Map (green=significant)')
axes[1, 1].set_xlabel('Cause →')
axes[1, 1].set_ylabel('← Effect')

plt.tight_layout()
plt.savefig('granger_causality.png', dpi=300)

print("\nResult: Identified causal chain: Catalyst temperature → Reaction rate → Product yield")
print("Fault diagnosis: Root cause of product yield decline is likely catalyst temperature anomaly")

3.7 Strategy for Real Process Application

💡 Best Practices for Implementation

• Phased Introduction : Try in order: statistical methods → machine learning → deep learning
• Threshold Tuning : Optimize thresholds with operational data to reduce false alarms
• Leverage Ensembles : Improve robustness by combining multiple methods
• Emphasize Visualization : Visualize detection reasons so operators can make decisions
• Continuous Learning : Retrain models regularly in response to process changes

Method Selection Guide

3.8 Summary

In this chapter, we implemented 8 anomaly detection and fault diagnosis techniques ranging from statistical methods to deep learning. By understanding the characteristics of each method and appropriately selecting and combining them according to process characteristics and operational conditions, you can build a high-precision anomaly detection system.

Skills Acquired

• ✅ Basic anomaly detection using statistical methods (Z-score, Modified Z-score)
• ✅ Efficient anomaly detection in high-dimensional data with Isolation Forest
• ✅ Novelty detection and decision boundary visualization with One-Class SVM
• ✅ Reconstruction error-based anomaly detection with Autoencoder
• ✅ Time series pattern learning and prediction error-based detection with LSTM
• ✅ Multi-class classification of fault types with Random Forest
• ✅ Robust anomaly detection system through ensemble voting
• ✅ Root cause analysis with Granger causality test

📚 Next Steps

In Chapter 4, we will learn process optimization and soft sensor techniques. We will implement quality prediction using machine learning models, real-time optimization, and construction of virtual sensors.

3.9 Exercises

Exercise 1 (Easy): Comparison of Statistical Anomaly Detection

Modify the code in Example 1 to compare the performance of Z-score and Modified Z-score on data that does not follow a normal distribution (e.g., log-normal distribution). Explain which is more robust.

💡 Hint

Think about the basic principles covered in the chapter examples.

📝 Sample Solution

Implementation approach:

• Step 1: [Key implementation point]
• Step 2: [Analysis or comparison]
• Step 3: [Validation and interpretation]

Exercise 2 (Medium): Autoencoder Architecture Optimization

For the Autoencoder in Example 4, conduct the following experiments to determine the optimal architecture:

• Compare encoding dimensions of 1, 2, 4, and 8
• Compare number of hidden layers: 1, 2, and 3
• Compare activation functions (ReLU, tanh, ELU)

💡 Hint

Consider the trade-offs between different approaches and parameter settings.

📝 Sample Solution

Implementation approach:

• Step 1: [Key implementation point]
• Step 2: [Analysis or comparison]
• Step 3: [Validation and interpretation]

Exercise 3 (Hard): Hybrid Anomaly Detection System

Design and implement a practical anomaly detection system that meets the following requirements:

1. Primary screening with statistical methods (reduce computational cost)
2. Detailed inspection with machine learning for suspicious data
3. Estimate causes with Granger causality when anomalies are detected
4. Set thresholds to keep false alarm rate below 5%

💡 Hint

Break down the problem into smaller steps and validate each component.

📝 Sample Solution

Implementation approach:

• Step 1: [Key implementation point]
• Step 2: [Analysis or comparison]
• Step 3: [Validation and interpretation]

💡 Hint

For Exercise 3, a cascade detection system is effective. By performing high-speed screening with statistical methods and verifying only suspicious data with deep learning, you can balance computational cost and accuracy. Draw ROC curves to determine optimal thresholds.

← Back to Chapter 2 Proceed to Chapter 4 →

References

1. Montgomery, D. C. (2019). Design and Analysis of Experiments (9th ed.). Wiley.
2. Box, G. E. P., Hunter, J. S., & Hunter, W. G. (2005). Statistics for Experimenters: Design, Innovation, and Discovery (2nd ed.). Wiley.
3. Seborg, D. E., Edgar, T. F., Mellichamp, D. A., & Doyle III, F. J. (2016). Process Dynamics and Control (4th ed.). Wiley.
4. McKay, M. D., Beckman, R. J., & Conover, W. J. (2000). “A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code.” Technometrics , 42(1), 55-61.

Disclaimer

• This content is provided solely for educational, research, and informational purposes and does not constitute professional advice (legal, accounting, technical warranty, etc.).
• This content and accompanying code examples are provided “AS IS” without any warranty, express or implied, including but not limited to merchantability, fitness for a particular purpose, non-infringement, accuracy, completeness, operation, or safety.
• The author and Tohoku University assume no responsibility for the content, availability, or safety of external links, third-party data, tools, libraries, etc.
• To the maximum extent permitted by applicable law, the author and Tohoku University shall not be liable for any direct, indirect, incidental, special, consequential, or punitive damages arising from the use, execution, or interpretation of this content.
• The content may be changed, updated, or discontinued without notice.
• The copyright and license of this content are subject to the stated conditions (e.g., CC BY 4.0). Such licenses typically include no-warranty clauses.