4. Exploratory data analysis for industrial datasets

Industrial datasets rarely arrive as clean, stationary, single-regime tables. A reactor historian export may contain 40 tags sampled at 1 Hz across three production campaigns, with a grade change in the middle, a sensor that drifted for two weeks, and a compressor that only starts when ambient temperature exceeds 30 °C. Exploratory data analysis (EDA) is the disciplined first pass that reveals these structures before any model is fitted. In this chapter we focus on four techniques that carry most of the diagnostic weight in process and manufacturing data: multivariate visualization, correlation analysis, lag plots, and operational regime identification.

Multivariate Visualization: Seeing More Than Two Tags at Once

Univariate histograms and time trends are necessary but insufficient. Most industrial faults and quality deviations are multivariate: each tag stays inside its own specification limit while the joint relationship between tags breaks down. The practical toolkit includes scatterplot matrices, parallel coordinates, and small multiples of pairwise scatter plots colored by a categorical variable such as product grade or shift.

A scatterplot matrix of, say, reactor temperature, pressure, feed rate, and coolant flow immediately exposes clusters, outliers, and nonlinear manifolds. Parallel coordinates are stronger when the number of variables exceeds six: each observation becomes a polyline crossing vertical axes, and regime changes appear as bundles of lines that separate between two axes. A useful convention is to normalize each axis to its 5th–95th percentile range so that one high-variance tag does not dominate the visual scale.

For high-dimensional data, project to two dimensions with PCA or UMAP and color by a known label (grade, campaign, fault flag). If the projection shows clean separation, a simple classifier will likely work; if it shows a continuous manifold, regime detection should be treated as a segmentation problem rather than classification.

Correlation Analysis: Beyond Pearson

Pearson correlation assumes linearity and is highly sensitive to outliers. In industrial data, both assumptions fail routinely. Use Spearman rank correlation for monotonic but nonlinear relationships, and always inspect the scatter plot behind any correlation coefficient above 0.6 or below −0.6. A correlation matrix should be accompanied by a p-value or bootstrap confidence interval; with 10,000 samples, a coefficient of 0.03 can be "significant" yet meaningless.

Two industrial-specific cautions matter. First, autocorrelation inflates apparent correlation between two slow-moving tags; two independent random walks will show |r| > 0.7 over a few hundred samples. Second, common-cause drivers create spurious pairwise correlations: feed rate drives both reactor temperature and downstream pressure, so temperature and pressure correlate even with no direct link. Partial correlation or a simple causal graph helps disentangle these.

Lag Plots: Quantifying Dynamic Response

Process variables respond to each other with delay. A lag plot of x(t) versus y(t+τ) for a range of τ reveals the transport delay and the direction of influence. The standard workflow is:

1. Detrend and, if needed, difference both series to remove slow drift.
2. Compute cross-correlation for τ from −N to +N (N typically 10–30 % of the dominant time constant).
3. Identify the τ* that maximizes |cross-correlation|; that is the estimated lag.
4. Inspect the scatter at τ* to confirm the relationship is monotonic and not driven by a few outliers.

For a temperature loop with a 90-second residence time sampled at 1 Hz, expect a peak near τ = 90 samples. If the peak appears at τ = 0 with a broad plateau, the sampling rate is too slow or the variables are dominated by a common slow driver. Lag plots also detect feedback: if x leads y at positive τ and y leads x at negative τ with similar magnitude, you are likely looking at a closed loop, not a causal chain.

Identifying Operational Regimes

An operational regime is a persistent mode of operation with distinct statistical behavior—startup, steady state at grade A, steady state at grade B, turndown, cleaning cycle. Regimes must be identified before building soft sensors or anomaly detectors, because a model trained on mixed regimes will average incompatible dynamics.

A practical detection procedure combines segmentation and clustering:

1. Feature windows: slide a window (e.g., 10 minutes) and compute mean, standard deviation, and slope for each tag.
2. Dimensionality reduction: apply PCA to the windowed features, retaining components explaining ≥ 90 % variance.
3. Segmentation: run change-point detection (PELT or binary segmentation) on the first few principal components to find candidate boundaries.
4. Clustering: cluster the segments with k-means or HDBSCAN; assign each segment a regime label.
5. Validation: overlay labels on the original time trends and check against operator logs, batch records, or grade-change events.

Typical failure modes: too-short windows produce noisy labels; too-long windows smear transitions; unnormalized features let one high-variance tag dominate. Always validate against an independent source—shift reports, DCS event logs, or lab sample timestamps.

Putting It Together: A Reproducible EDA Checklist

The following checklist is designed to be run on any historian export before modeling. It assumes a pandas DataFrame df with a DatetimeIndex and numeric tag columns.

import pandas as pd, numpy as np
from scipy.stats import spearmanr
from scipy.signal import correlate

# 1. Multivariate overview
pd.plotting.scatter_matrix(df.sample(2000), alpha=0.3, figsize=(12,12))

# 2. Correlation (Spearman, robust to nonlinearity)
corr = df.corr(method='spearman')
# Flag pairs with |r| > 0.6 for manual scatter inspection
high = corr.where(np.triu(np.ones(corr.shape), k=1).astype(bool)).stack()
print(high[high.abs() > 0.6])

# 3. Lag analysis for a candidate pair
x = df['reactor_temp'].diff().dropna()
y = df['product_quality'].diff().dropna()
cc = correlate(x - x.mean(), y - y.mean(), mode='full')
lags = np.arange(-len(x)+1, len(x))
tau_star = lags[np.argmax(np.abs(cc))]
print(f"Estimated lag: {tau_star} samples")

# 4. Regime segmentation
from ruptures import Pelt
algo = Pelt(model="rbf").fit(df[['PC1','PC2']].values)
breakpoints = algo.predict(pen=10)
print(f"Change points at indices: {breakpoints}")

Run this on a 24-hour window first, then on the full campaign. If the lag estimate shifts by more than 20 % between windows, the process is non-stationary and regime-aware modeling is mandatory. If the correlation matrix changes sign between windows, you are almost certainly mixing regimes.

Interpreting the Results

After running the checklist, you should be able to answer four questions: (1) Are there distinct clusters or manifolds in the multivariate space? (2) Which tag pairs are strongly coupled, and is that coupling linear or monotonic? (3) What is the dominant transport delay between key inputs and quality outputs? (4) How many operational regimes exist, and where are the boundaries? These answers determine whether the next step is a global model, a regime-switching model, or a dynamic model with explicit lag terms. EDA is not a preliminary formality—it is the step that decides which modeling assumptions are defensible.