Triển khai Predictive Maintenance 4.0: IoT và Machine Learning trong bảo trì thiết bị

┌─────────────────────────────────────────────────────────┐
│                    Cloud/Edge Analytics                  │
│  ┌──────────────┐  ┌──────────────┐  ┌──────────────┐  │
│  │   Machine    │  │  Predictive  │  │  Dashboard   │  │
│  │   Learning   │  │    Models    │  │  & Alerts    │  │
│  └──────────────┘  └──────────────┘  └──────────────┘  │
└─────────────────────────────────────────────────────────┘
                            ▲
                            │ MQTT/HTTP/OPC UA
                            │
┌─────────────────────────────────────────────────────────┐
│                    IoT Gateway/Edge                      │
│  ┌──────────────┐  ┌──────────────┐  ┌──────────────┐  │
│  │ Data Buffer  │  │ Edge ML      │  │ Protocol     │  │
│  │ & Filtering  │  │ Inference    │  │ Translate    │  │
│  └──────────────┘  └──────────────┘  └──────────────┘  │
└─────────────────────────────────────────────────────────┘
                            ▲
                ┌───────────┼───────────┐
                │           │           │
        ┌───────▼──┐  ┌────▼────┐  ┌──▼──────┐
        │Vibration │  │  Temp   │  │ Current │
        │ Sensor   │  │ Sensor  │  │ Sensor  │
        └──────────┘  └─────────┘  └─────────┘
                │          │           │
        ┌───────▼──────────▼───────────▼────────┐
        │        Industrial Equipment            │
        │  (Motor / Pump / Gearbox / Bearing)    │
        └────────────────────────────────────────┘

// Arduino/PlatformIO code for IoT sensor node
#include 
#include 
#include 
#include 
#include 
#include 

// Sensor objects
Adafruit_ADXL345_Unified accel = Adafruit_ADXL345_Unified(12345);
OneWire oneWire(TEMP_PIN);
DallasTemperature tempSensors(&oneWire);

// MQTT client
WiFiClient espClient;
PubSubClient mqtt(espClient);

// Configuration
const char* mqtt_server = "mqtt.predictive-maintenance.com";
const char* device_id = "MOTOR_101";
const int SAMPLE_RATE = 10000; // 10kHz for vibration
const int PUBLISH_INTERVAL = 60000; // 1 minute

void setup() {
    Serial.begin(115200);
    
    // Init sensors
    if(!accel.begin()) {
        Serial.println("ADXL345 not found!");
    }
    accel.setRange(ADXL345_RANGE_16_G);
    
    tempSensors.begin();
    
    // Connect WiFi
    WiFi.begin(WIFI_SSID, WIFI_PASSWORD);
    while (WiFi.status() != WL_CONNECTED) {
        delay(500);
    }
    
    // Connect MQTT
    mqtt.setServer(mqtt_server, 1883);
}

void loop() {
    // Sample vibration at high frequency
    float vibration_rms = sampleVibration(1000); // 1000 samples
    
    // Read temperature
    tempSensors.requestTemperatures();
    float temperature = tempSensors.getTempCByIndex(0);
    
    // Read current (via analog input from current sensor)
    float current = analogRead(CURRENT_PIN) * CURRENT_SCALE;
    
    // Calculate features
    float vibration_peak = calculatePeak();
    float vibration_crest = vibration_peak / vibration_rms;
    
    // Publish to MQTT
    String payload = createJSON(vibration_rms, vibration_peak, 
                                vibration_crest, temperature, current);
    mqtt.publish("sensors/motor_101/data", payload.c_str());
    
    delay(PUBLISH_INTERVAL);
}

float sampleVibration(int samples) {
    float sum_squares = 0;
    
    for(int i = 0; i < samples; i++) {
        sensors_event_t event;
        accel.getEvent(&event);
        
        // Calculate magnitude
        float magnitude = sqrt(event.acceleration.x * event.acceleration.x +
                              event.acceleration.y * event.acceleration.y +
                              event.acceleration.z * event.acceleration.z);
        
        sum_squares += magnitude * magnitude;
        
        delayMicroseconds(100); // 10kHz sampling
    }
    
    return sqrt(sum_squares / samples); // RMS
}

String createJSON(float vib_rms, float vib_peak, float crest, 
                  float temp, float current) {
    String json = "{";
    json += "\"device\":\"" + String(device_id) + "\",";
    json += "\"timestamp\":" + String(millis()) + ",";
    json += "\"vibration_rms\":" + String(vib_rms, 4) + ",";
    json += "\"vibration_peak\":" + String(vib_peak, 4) + ",";
    json += "\"crest_factor\":" + String(crest, 4) + ",";
    json += "\"temperature\":" + String(temp, 2) + ",";
    json += "\"current\":" + String(current, 2);
    json += "}";
    return json;
}

# Python code for feature extraction
import numpy as np
from scipy import signal
from scipy.stats import kurtosis, skew

def extract_features(vibration_signal, fs=10000):
    """
    Extract time-domain and frequency-domain features
    
    Args:
        vibration_signal: 1D array of acceleration data
        fs: Sampling frequency
    
    Returns:
        Dictionary of features
    """
    features = {}
    
    # Time-domain features
    features['rms'] = np.sqrt(np.mean(vibration_signal**2))
    features['peak'] = np.max(np.abs(vibration_signal))
    features['crest_factor'] = features['peak'] / features['rms']
    features['kurtosis'] = kurtosis(vibration_signal)
    features['skewness'] = skew(vibration_signal)
    features['shape_factor'] = features['rms'] / np.mean(np.abs(vibration_signal))
    
    # Frequency-domain features (FFT)
    fft = np.fft.fft(vibration_signal)
    freqs = np.fft.fftfreq(len(vibration_signal), 1/fs)
    psd = np.abs(fft)**2
    
    # Find dominant frequencies (bearing fault frequencies)
    dominant_freq = freqs[np.argmax(psd[:len(psd)//2])]
    features['dominant_frequency'] = dominant_freq
    
    # Bearing fault frequency bands (example for 1800 RPM motor)
    bpfo = 1800/60 * 3.5  # Ball Pass Frequency Outer race
    bpfi = 1800/60 * 5.5  # Ball Pass Frequency Inner race
    bsf = 1800/60 * 2.3   # Ball Spin Frequency
    
    features['bpfo_magnitude'] = extract_band_power(psd, freqs, bpfo, 10)
    features['bpfi_magnitude'] = extract_band_power(psd, freqs, bpfi, 10)
    features['bsf_magnitude'] = extract_band_power(psd, freqs, bsf, 5)
    
    # Envelope analysis for bearing faults
    analytic_signal = signal.hilbert(vibration_signal)
    envelope = np.abs(analytic_signal)
    features['envelope_rms'] = np.sqrt(np.mean(envelope**2))
    
    return features

def extract_band_power(psd, freqs, center_freq, bandwidth):
    """Extract power in specific frequency band"""
    idx = np.where((freqs >= center_freq - bandwidth) & 
                   (freqs <= center_freq + bandwidth))
    return np.sum(psd[idx])

# Train anomaly detection model
from sklearn.ensemble import IsolationForest
from sklearn.preprocessing import StandardScaler
import pandas as pd

# Load training data (normal operation only)
df_train = pd.read_csv('normal_operation_data.csv')

# Feature columns
feature_cols = ['rms', 'peak', 'crest_factor', 'kurtosis', 
                'temperature', 'current', 'bpfo_magnitude']

X_train = df_train[feature_cols]

# Normalize features
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)

# Train Isolation Forest
model = IsolationForest(
    contamination=0.01,  # 1% expected anomalies
    max_samples=256,
    random_state=42
)
model.fit(X_train_scaled)

# Inference on new data
def predict_anomaly(new_features):
    """
    Returns:
        1: Normal
        -1: Anomaly (potential fault)
    """
    X_new = scaler.transform([new_features])
    prediction = model.predict(X_new)
    anomaly_score = model.score_samples(X_new)
    
    return prediction[0], anomaly_score[0]

# Example usage
features = extract_features(live_vibration_signal)
feature_vector = [features[col] for col in feature_cols]
is_normal, score = predict_anomaly(feature_vector)

if is_normal == -1:
    print(f"⚠️ ANOMALY DETECTED! Score: {score:.4f}")
    print("Recommended action: Schedule inspection")

# LSTM model for RUL prediction
import tensorflow as tf
from tensorflow.keras import layers, models

def build_rul_model(sequence_length, n_features):
    """
    Build LSTM model for Remaining Useful Life prediction
    
    Args:
        sequence_length: Number of time steps in sequence
        n_features: Number of features per time step
    """
    model = models.Sequential([
        layers.LSTM(128, return_sequences=True, 
                   input_shape=(sequence_length, n_features)),
        layers.Dropout(0.2),
        layers.LSTM(64, return_sequences=True),
        layers.Dropout(0.2),
        layers.LSTM(32),
        layers.Dense(16, activation='relu'),
        layers.Dense(1)  # Output: RUL in hours
    ])
    
    model.compile(optimizer='adam', loss='mse', metrics=['mae'])
    return model

# Training data preparation
def create_sequences(data, rul_values, sequence_length=50):
    """
    Create sequences for LSTM training
    
    Args:
        data: Feature matrix (samples × features)
        rul_values: RUL labels for each sample
        sequence_length: Window size
    """
    X, y = [], []
    
    for i in range(len(data) - sequence_length):
        X.append(data[i:i+sequence_length])
        y.append(rul_values[i+sequence_length])
    
    return np.array(X), np.array(y)

# Train model
model = build_rul_model(sequence_length=50, n_features=7)
history = model.fit(
    X_train, y_train,
    validation_split=0.2,
    epochs=100,
    batch_size=64,
    callbacks=[
        tf.keras.callbacks.EarlyStopping(patience=10),
        tf.keras.callbacks.ModelCheckpoint('best_rul_model.h5', 
                                          save_best_only=True)
    ]
)

# Prediction
def predict_rul(model, recent_data):
    """Predict RUL from recent sensor data"""
    # recent_data: Last 50 time steps of features
    rul_hours = model.predict(recent_data.reshape(1, 50, 7))
    return rul_hours[0][0]

// InfluxDB query example
SELECT 
    mean("vibration_rms") as "Vibration",
    mean("temperature") as "Temperature",
    mean("current") as "Current"
FROM "sensor_data"
WHERE 
    "device" = 'MOTOR_101'
    AND time >= now() - 24h
GROUP BY time(1m)

// Alerting rules in Grafana
WHEN avg() OF query(A, 5m, now) IS ABOVE 0.5  // Anomaly threshold
THEN
    ALERT "Motor 101 Abnormal Vibration"
    MESSAGE "Vibration RMS exceeded threshold. Potential bearing fault."
    SEND TO maintenance_team_webhook

# Python backend API
from fastapi import FastAPI
from pydantic import BaseModel
import requests

app = FastAPI()

class AnomalyAlert(BaseModel):
    device_id: str
    anomaly_type: str
    severity: str
    predicted_rul: float

@app.post("/api/anomaly/detected")
async def handle_anomaly(alert: AnomalyAlert):
    """Create maintenance work order"""
    
    # Determine urgency
    if alert.predicted_rul < 24:  # Less than 24 hours
        priority = "URGENT"
    elif alert.predicted_rul < 168:  # Less than 1 week
        priority = "HIGH"
    else:
        priority = "NORMAL"
    
    # Create work order in CMMS
    work_order = {
        "title": f"{alert.anomaly_type} detected on {alert.device_id}",
        "description": f"Predictive maintenance alert. Estimated RUL: {alert.predicted_rul:.1f} hours",
        "asset_id": alert.device_id,
        "priority": priority,
        "assigned_to": "maintenance_team",
        "due_date": calculate_due_date(alert.predicted_rul)
    }
    
    response = requests.post(
        "https://cmms.company.com/api/work_orders",
        json=work_order,
        headers={"Authorization": f"Bearer {CMMS_API_KEY}"}
    )
    
    # Send notification
    send_notification(alert, priority)
    
    return {"status": "Work order created", "id": response.json()['id']}