7.9. Diagnosing a known fault with batch PLS: the SBR reactor

This is the second of three case studies on batch data. The batches now have their final quality measured, so the model is a PLS from the trajectories to that quality, and it can be asked what the quality will be while a batch is still running.

The first set up the batchwise unfolding, the scores, the SPE and the contribution plots, which are used here without being explained again. It had to work out which batches were unusual without knowing; here the fault is known, so each diagnosis can be checked against it.

Styrene-butadiene rubber (SBR) is made by emulsion polymerization in a batch reactor. Six trajectories of the reactor are used here: the reactor, cooling-water and jacket temperatures, the latex density, the conversion, and the energy released.

Five quality attributes of the latex are measured at the end: composition, particle size, branching, cross-linking and polydispersity. Branching and cross-linking are one measurement in two units, in a fixed ratio, so the quality block carries four independent attributes rather than five.

The 53 batches were simulated from a first-principles model of the reactor (Nomikos and MacGregor, 1994; Nomikos, 1995), so the fault is known. Batch 37 received 30% more organic impurity in the butadiene feed from its very start, batch 34 received 50% more from midway through. Simulated data lets a model be checked against what is known to have happened, before it is trusted on plant data.

The model is a PLS from the unfolded trajectories to the five attributes. Three questions are asked of it:

• Does it single out the two faulty batches?

• Does it say what went wrong, and when?

• Could it have said so while they were still running?

7.9.1. The data

The SBR batch reactor dataset is a workbook of two sheets: the trajectories of 53 batches over 200 samples, and the quality attributes, one row per batch.

Three of its nine trajectories, the two feed flow rates and the feed temperature, carry only the noise the simulation adds and nothing that differs between batches. Scaled to unit variance they would weigh as much as the informative ones, so the model uses the six trajectories of the reactor itself, which load_sbr reads from the workbook.

Python, 23 lines

import numpy as np
import pandas as pd
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from process_improve.batch import (
    BatchMonitor, BatchPLS, contribution_at_time_plot, load_sbr,
    online_monitoring_plot, time_varying_loading_plot, unfolded_contribution_plot)
from process_improve.multivariate import PLS
from process_improve.univariate import median_absolute_deviation

# https://openmv.net/file/sbr-batch-reactor.xlsx
sbr = load_sbr()
trajectories = {batch_id: batch[sbr.trajectory_tags]
                for batch_id, batch in sbr.X.items()}
quality = sbr.Y
print(len(trajectories), "batches;", sbr.trajectory_tags,
      "; quality block", quality.shape)
# 53 batches; ['ReactorTemp', 'CoolingTemp', 'JacketTemp', 'LatexDensity',
# 'Conversion', 'EnergyReleased'] ; quality block (53, 5)

ratio = quality["CrossLinking"] / quality["Branching"]
print(f"CrossLinking / Branching: {ratio.min():.3f} to {ratio.max():.3f}")
# CrossLinking / Branching: 3.877 to 3.878

Python, 20 lines

# figure colours
GREY, ORANGE, AQUA, BLUE, PURPLE = "#c8c8c8", "#c55a11", "#1baf7a", "#1f3d7a", "#6f42c1"

def overlay(batches, tag, highlight):
    """One tag for every batch in grey, with the batches in `highlight` (id -> colour)
    drawn on top."""
    fig = go.Figure()
    for batch_id, batch in batches.items():
        if batch_id not in highlight:
            fig.add_trace(go.Scatter(y=batch[tag], mode="lines",
                                     line=dict(color=GREY, width=1), showlegend=False))
    for batch_id, colour in highlight.items():
        fig.add_trace(go.Scatter(y=batches[batch_id][tag], mode="lines",
                                 name=f"batch {batch_id}",
                                 line=dict(color=colour, width=3)))
    fig.update_layout(title=tag, xaxis_title="Sample [aligned time]", height=320)
    return fig

for tag in sbr.trajectory_tags:
    overlay(trajectories, tag, {34: ORANGE, 37: AQUA}).show()

The six trajectories of all 53 batches (grey) with batch 34 (orange) and batch 37 (aqua) drawn on top. Batch 37 converts more slowly than the others from the start. Batch 34, the one whose fault starts midway through, follows the others until about sample 100, after which its cooling-water and jacket temperatures rise above the band of the other batches and its energy released falls below it.

The overlay already shows both faults: batch 37 from the start, in the conversion, and batch 34 from about sample 100, in the two service temperatures (cooling water and jacket) and the energy released. Neither shows in the reactor temperature, which is held within a narrow range in every batch.

7.9.2. The batch PLS model

BatchPLS unfolds each batch batchwise into one row of 6 tags by 200 samples, 1200 columns, centres and scales every column, and fits a PLS from that row to the five quality attributes, also centred and scaled. Two components are used. As in the first case study, this first model is to see what is going on, not to settle the number of components; the predicted quality section asks what the second one buys, against the usual cross-validation.

Python, 51 lines

model = BatchPLS(n_components=2).fit(trajectories, quality)
# R2 of the quality block, per component
r2y = np.diff([0.0, *model.r2_cumulative_])
# R2 of the trajectories, per component
r2x = np.diff([0.0, *model.r2_per_variable_.mean(axis=0)])
print("R2X per component:", r2x.round(3), " R2Y per component:", r2y.round(3))
# R2X per component: [0.245 0.127]  R2Y per component: [0.653 0.069]
spe = model.spe_.iloc[:, -1]
t1, t2 = model.scores_.iloc[:, 0], model.scores_.iloc[:, 1]
marked = {34: ORANGE, 37: AQUA, 4: PURPLE}
others = [batch_id for batch_id in t1.index if batch_id not in marked]
# The marker area is the batch's squared residual. `sizemode="area"` gives the area,
# not the diameter, to the value; squaring spreads a factor of two in SPE over a factor
# of four in area, and the sum of squared residuals is what adds up over the cells in
# any case.
area = dict(sizemode="area", sizeref=2 * (spe**2).max() / 26**2, sizemin=3)
ex, ey = model.ellipse_coordinates(score_horiz=1, score_vert=2, conf_level=0.95)
fig = go.Figure()
fig.add_trace(go.Scatter(x=ex, y=ey, mode="lines", name="95% confidence ellipse",
                         line=dict(color=GREY, dash="dash")))
fig.add_trace(go.Scatter(x=t1.loc[others], y=t2.loc[others], mode="markers",
                         showlegend=False, text=others, hovertemplate="batch %{text}",
                         marker=dict(size=spe.loc[others] ** 2, color=BLUE, **area)))
for batch_id, colour in marked.items():
    fig.add_trace(go.Scatter(x=[t1.loc[batch_id]], y=[t2.loc[batch_id]], mode="markers",
                             name=f"batch {batch_id}",
                             marker=dict(size=[spe.loc[batch_id] ** 2],
                                         color=colour, **area,
                                         line=dict(color="#404040", width=1.5))))
fig.update_layout(xaxis_title=f"t1 [R2X {r2x[0]:.1%}]",
                  yaxis_title=f"t2 [R2X {r2x[1]:.1%}]", height=520).show()
print(f"SPE rank of batch 37: {int(spe.rank().loc[37])} of {len(spe)};",
      f"batch 34: {int(spe.rank().loc[34])}")            # 1 = the smallest residual
# SPE rank of batch 37: 1 of 53; batch 34: 12
print(f"t1 rank of batch 37: {int(t1.rank().loc[37])} of {len(t1)};",
      f"t2 rank of batch 34: {int(t2.rank().loc[34])}")   # 1 = the lowest score
# t1 rank of batch 37: 1 of 53; t2 rank of batch 34: 53
# the batch nearest the average quality
standardised = (quality - quality.mean()) / quality.std(ddof=1)
distance = (standardised ** 2).sum(axis=1).pow(0.5) / np.sqrt(quality.shape[1])
print("nearest the average quality:", distance.nsmallest(2).round(2).to_dict())
# nearest the average quality: {4: 0.16, 27: 0.26}
for batch_id in (34, 37):
    t2_end = model.hotellings_t2_.loc[batch_id].iloc[-1]
    spe_end = model.spe_.loc[batch_id].iloc[-1]
    t2_lim = model.hotellings_t2_limit(conf_level=0.95)
    spe_lim = model.spe_limit(conf_level=0.95)
    print(f"batch {batch_id}: T2 = {t2_end:.1f} (limit {t2_lim:.1f}),",
          f"SPE = {spe_end:.1f} (limit {spe_lim:.1f})")
    # batch 34: T2 = 28.2 (limit 6.6), SPE = 23.1 (limit 34.6)
    # batch 37: T2 = 19.2 (limit 6.6), SPE = 18.7 (limit 34.6)

The first component explains 65.3% of the variance in the quality block and the second 6.9%; their shares of the variance in the trajectories are on the axes of the score plot.

Scores of the batch PLS model. Batch 37 (aqua) has the lowest \(t_1\) of all batches and batch 34 (orange) the highest \(t_2\); both lie far outside the 95% confidence ellipse. The area of each marker is that batch’s squared residual, on the scale of the three grey circles in the legend, so the plot answers both questions asked of a batch: the two faulty ones are extreme along the components and carry small residuals. Batch 4 (purple) is the batch nearest the average quality, used below to show the prediction while the batch runs.

What the plot cannot show is that both faulty batches are inside this fit. Each is extreme along a component that its own deviation helped to define, which is why the fault lands in the scores and leaves so little in the residual, and why the two largest values of \(T^2\) (Hotelling’s statistic) belong to two of the smallest residuals. The 95% limits inherit the same circularity, since they are computed from the 53 batches that include the two faults being judged.

The SPE answers the other question about a batch, how far it sits away from the components. Drawing it against \(T^2\) puts both questions in one figure, each axis carrying its own 95% limit.

Python, 29 lines

def influence_plot(model, highlight, labels, conf_level=0.95):
    """Hotelling's T2 against SPE, one dot per batch, with both limits drawn."""
    t2, spe = model.hotellings_t2_.iloc[:, -1], model.spe_.iloc[:, -1]
    others = [batch_id for batch_id in t2.index if batch_id not in highlight]
    fig = go.Figure()
    fig.add_trace(go.Scatter(x=t2.loc[others], y=spe.loc[others], mode="markers",
                             showlegend=False, marker=dict(size=8, color=BLUE),
                             text=others, hovertemplate="batch %{text}"))
    for batch_id, colour in highlight.items():
        fig.add_trace(go.Scatter(x=[t2.loc[batch_id]], y=[spe.loc[batch_id]],
                                 mode="markers", name=f"batch {batch_id}",
                                 marker=dict(size=12, color=colour)))
    fig.add_vline(x=model.hotellings_t2_limit(conf_level=conf_level),
                  line_dash="dash", line_color=GREY)
    fig.add_hline(y=model.spe_limit(conf_level=conf_level),
                  line_dash="dash", line_color=GREY)
    for batch_id in labels:
        fig.add_annotation(x=t2.loc[batch_id], y=spe.loc[batch_id], text=str(batch_id),
                           showarrow=False, xshift=13, yshift=9)
    fig.update_layout(xaxis_title="Hotelling's T\u00b2", yaxis_title="SPE", height=420)
    return fig

spe, t2 = model.spe_.iloc[:, -1], model.hotellings_t2_.iloc[:, -1]
above_spe = sorted(spe.index[spe > model.spe_limit(conf_level=0.95)])
print("above the SPE limit:", above_spe, "  above the T2 limit:",
      sorted(t2.index[t2 > model.hotellings_t2_limit(conf_level=0.95)]))
# above the SPE limit: [8, 15, 16]   above the T2 limit: [34, 37]
influence_plot(model, highlight={34: ORANGE, 37: AQUA},
               labels=[34, 37, *above_spe]).show()

Hotelling’s \(T^2\) against the SPE for every batch, with both 95% limits. Batches 34 (orange) and 37 (aqua) are far to the right and below the SPE limit. The batches the SPE flags, 8, 15 and 16, are different batches, in the upper left.

A 95% limit is crossed by about three of 53 batches when nothing is wrong, so the three the SPE flags are what the limit allows, not evidence of a fault.

The figure cannot show why the two known faults are not among them. The SPE sums the residuals of all 1200 cells of a batch after the components have taken their share. A deviation along a component is absorbed by the scores, and a deviation confined to a few cells is small against a sum over 1200.

The scores say a batch moved in a direction the model knows; the SPE, that it moved in one the model does not.

7.9.3. Where the model explains the trajectories

Every unfolded column has its own \(R^2\), so the fit reads per tag and per time sample. So do the weights \(\mathbf{w}_1\) and \(\mathbf{w}_2\), each a vector of 1200 entries, which time_varying_loading_plot draws as six curves over the batch. PLS weights play the role PCA loadings do, and how the PLS model is calculated explains the difference. The contribution plots later on the page are built from the direct weights \(\mathbf{W}(\mathbf{P}'\mathbf{W})^{-1}\), the combination the PLS chapter prefers for interpretation, rather than the \(\mathbf{w}\) drawn here.

Python, 14 lines

# rows = tags, columns = time
r2_grid = model.r2_per_variable_.iloc[:, -1].unstack(level="sequence")
fig = go.Figure()
for tag, row in r2_grid.iterrows():
    fig.add_trace(go.Scatter(x=row.index, y=row.values, mode="lines", name=tag))
fig.update_layout(title="R2 of each (tag, time) cell after two components",
                  xaxis_title="Sample [aligned time]", yaxis_title="R2", height=360)
fig.show()
print("R2 per tag, averaged over time:", r2_grid.mean(axis=1).round(2).to_dict())
# R2 per tag, averaged over time: {'Conversion': 0.75, 'CoolingTemp': 0.23,
# 'EnergyReleased': 0.26, 'JacketTemp': 0.24, 'LatexDensity': 0.67,
# 'ReactorTemp': 0.08}
time_varying_loading_plot(model, component=1).show()
time_varying_loading_plot(model, component=2).show()

\(R^2\) of every (tag, time) cell after two components, one panel per tag. Latex density and conversion are well explained from about sample 50 onwards; the two temperatures and the energy released only in the second half of the batch; the reactor temperature hardly at all.

Latex density and conversion are the trajectories the model uses most, with an \(R^2\) averaged over the batch of 0.67 and 0.75, against 0.23 to 0.26 for the two service temperatures and the energy released, and 0.08 for the reactor temperature. Every curve is low at the start. All batches begin alike, so after centring and scaling the first samples hold little but noise.

Time-varying weights \(\mathbf{w}_1\) (blue) and \(\mathbf{w}_2\) (orange) of the two components, one panel per tag. The first is dominated by latex density and conversion over the whole batch; the second by the cooling-water and jacket temperatures, the energy released, the latex density and the conversion in the second half of the batch, with opposite signs.

The weights say what a score means:

• a low \(t_1\) is a below-average latex density and conversion throughout the batch;

• a high \(t_2\) is cooling-water and jacket temperatures above their average trajectories over the second half, with the energy released, the latex density and the conversion below theirs, and the reactor temperature carrying weights that change sign from sample to sample, so the component says nothing consistent about it.

A component says these tags move together, not which of them drives the others. Those are predictions about batches 37 and 34 respectively, and the contribution plots test them.

7.9.4. Batch 37: the fault from the start

Python, 22 lines

# the 53 x 1200 matrix the model was fitted on
scaled = model.unfold_and_scale(trajectories)
t1 = model.score_contributions(scaled, component=1)
unfolded_contribution_plot(t1, batch_id=37).show()
unfolded_contribution_plot(t1, batch_id=37, by_tag=True).show()
print("batch 37, t1 contributions per tag:",
      t1.loc[37].groupby(level="tag", sort=False).sum().round(1).to_dict())
# batch 37, t1 contributions per tag: {'Conversion': -33.4, 'CoolingTemp': -4.2,
# 'EnergyReleased': -5.2, 'JacketTemp': -4.3, 'LatexDensity': -25.5,
# 'ReactorTemp': -1.3}

def share_per_fifth(row):
    """Split a contribution vector into five equal time blocks and give each block's
    share of the total."""
    by_time = row.groupby(level="sequence").sum()
    fifths = by_time.groupby(np.arange(len(by_time)) * 5 // len(by_time)).sum()
    return (fifths / fifths.sum()).round(2).tolist()

print("batch 37: share of the t1 contribution per fifth of the batch:",
      share_per_fifth(t1.loc[37]))
# batch 37: share of the t1 contribution per fifth of the batch:
# [0.15, 0.18, 0.19, 0.26, 0.22]

Top: the contribution of each (tag, time) cell of batch 37 to \(t_1\). Middle: the same contributions summed per tag. Bottom: summed per sample. Conversion and latex density carry the contribution, and every part of the batch contributes.

As in the first case study, the contribution vector has one entry per (tag, time) cell, here \(K = 6\) tags by \(J = 200\) time samples, and summing it per tag or per sample says which trajectory and which part of the batch the score came from.

Batch 37 sits at the low end of \(t_1\) because its conversion and latex density were below average, and the contribution is spread over the whole batch. Each fifth carries between 15% and 26% of the total, which is what a fault present from the start looks like. The raw trajectory overlay agrees. The impurity slowed the reaction from the moment the batch began.

7.9.5. Batch 34: the same fault, from the middle of the batch

Python, 12 lines

t2 = model.score_contributions(scaled, component=2)
unfolded_contribution_plot(t2, batch_id=34).show()
contribution_at_time_plot(t2, k=120, batch_id=34).show()
print("batch 34, t2 contributions per tag:",
      t2.loc[34].groupby(level="tag", sort=False).sum().round(1).to_dict())
# batch 34, t2 contributions per tag: {'Conversion': 8.1, 'CoolingTemp': 12.2,
# 'EnergyReleased': 14.3, 'JacketTemp': 12.3, 'LatexDensity': 7.2,
# 'ReactorTemp': 4.0}
print("batch 34: share of the t2 contribution per fifth of the batch:",
      share_per_fifth(t2.loc[34]))
# batch 34: share of the t2 contribution per fifth of the batch:
# [0.06, 0.09, 0.17, 0.39, 0.29]

Top: the contribution of each (tag, time) cell of batch 34 to \(t_2\). Middle: the same contributions summed per tag. Bottom: summed per sample. The energy released and the two temperatures lead, and the contribution per sample is small until sample 100 and large after it.

Batch 34 is high on \(t_2\) through the energy released and the two service temperatures, with the conversion and latex density behind them. Timing distinguishes it from batch 37. The first two fifths of the batch carry 15% of the \(t_2\) contribution and the last two fifths 68%. contribution_at_time_plot shows the same three tags carrying the deviation at sample 120.

The raw data can be asked the same question directly, tag by tag: at which sample does each trajectory of a faulty batch leave the band of the other batches?

The code below answers that for batches 34 and 37. Each tag becomes a \(z\) value at every sample: a distance from a centre, divided by a spread, both taken from the 51 batches with no known fault. Two versions are shown:

• classical: centred on the mean of those 51 at that sample, divided by their standard deviation;

• robustified: centred on their median, divided by their median absolute deviation, scaled to agree with the standard deviation on normally distributed values, then EWMA-smoothed.

A tag has left the band when its \(z\) value stays beyond plus or minus two for 20 samples in a row, a tenth of the batch.

Twenty samples is longer than the reactor temperature ever stays outside the band in a normal batch, but not longer than the other tags do: the latex density and the conversion of a normal batch can stay outside it for 50 samples or more. A sustained departure therefore says where to look, not that there is a fault.

Python, 96 lines

others = np.stack([batch.to_numpy() for batch_id, batch in trajectories.items()
                   if batch_id not in (34, 37)])
z = {batch_id: (trajectories[batch_id].to_numpy() - others.mean(axis=0))
               / others.std(axis=0, ddof=1)
     for batch_id in (37, 34)}
# 1.4826 x MAD
spread = median_absolute_deviation(others, axis=0, scale="normal")
# EWMA-smoothed
z_robust = {batch_id: pd.DataFrame((trajectories[batch_id].to_numpy()
                                    - np.median(others, axis=0)) / spread)
            .ewm(alpha=0.3, adjust=False).mean().to_numpy()
            for batch_id in (37, 34)}

def sustained_departure(z_batch, n_sd=2.0, run=20):
    """First sample, counting from one, from which each tag stays more than `n_sd`
    scale units away for `run` samples."""
    # 1 where the tag is outside the band, 0 inside
    outside = (np.abs(z_batch) > n_sd).astype(int)
    onset = {}
    for j, tag in enumerate(sbr.trajectory_tags):
        # Convolving the 0/1 column with `run` ones is a moving sum over every
        # window of `run` consecutive samples: entry i is how many of samples i to
        # i+run-1 lie outside the band. "valid" keeps only the windows that fit
        # entirely inside the batch, so entry i starts at sample i. A window summing
        # to `run` is an unbroken stretch, and argmax finds the first True, which is
        # the sample the stretch starts at. The +1 counts samples from one, as the
        # rest of the page does. Without a True, argmax returns 0, so test `any()`
        # first and report None instead.
        runs = np.convolve(outside[:, j], np.ones(run, dtype=int), mode="valid") == run
        onset[tag] = int(runs.argmax()) + 1 if runs.any() else None
    return onset

for batch_id in (37, 34):
    print(f"batch {batch_id}, first sustained departure:",
          sustained_departure(z[batch_id]))
    # batch 37, first sustained departure: {'ReactorTemp': None, 'CoolingTemp': None,
    # 'JacketTemp': None, 'LatexDensity': 14, 'Conversion': 10,
    # 'EnergyReleased': None}
    # batch 34, first sustained departure: {'ReactorTemp': None, 'CoolingTemp': 104,
    # 'JacketTemp': 105, 'LatexDensity': 130, 'Conversion': 124,
    # 'EnergyReleased': 106}
    print("  robust, smoothed:", sustained_departure(z_robust[batch_id]))
    #   robust, smoothed: {'ReactorTemp': None, 'CoolingTemp': 1, 'JacketTemp': 1,
    #   'LatexDensity': 16, 'Conversion': 1, 'EnergyReleased': None}
    #   robust, smoothed: {'ReactorTemp': None, 'CoolingTemp': 106,
    #   'JacketTemp': 108, 'LatexDensity': 135, 'Conversion': 127,
    #   'EnergyReleased': 108}

# Why 20: it has to outlast what the noise alone produces, and the answer must not
# depend on it.
def longest_run(flags):
    """Length of the longest unbroken stretch of True."""
    best = current = 0
    for flag in flags:
        current = current + 1 if flag else 0
        best = max(best, current)
    return best

quiet = [np.abs((batch.to_numpy() - others.mean(axis=0))
                / others.std(axis=0, ddof=1)) > 2
         for batch_id, batch in trajectories.items() if batch_id not in (34, 37)]
print("longest run outside the band with no fault:",
      {tag: max(longest_run(batch[:, j]) for batch in quiet)
       for j, tag in enumerate(sbr.trajectory_tags)})
# longest run outside the band with no fault: {'ReactorTemp': 9, 'CoolingTemp': 36,
# 'JacketTemp': 37, 'LatexDensity': 51, 'Conversion': 72, 'EnergyReleased': 18}
print("run lengths giving the same onsets:",
      [run for run in (15, 20, 25, 30, 40)
       if all(sustained_departure(z[b], run=run) == sustained_departure(z[b])
              for b in (37, 34))])
# run lengths giving the same onsets: [20, 25, 30, 40]

fig = make_subplots(rows=2, cols=6, subplot_titles=sbr.trajectory_tags,
                    shared_yaxes=True)
for row, (batch_id, colour) in enumerate([(37, AQUA), (34, ORANGE)], start=1):
    for col in range(6):
        # signed distances, not their absolute value: the sign says whether the tag
        # ran above or below the others
        fig.add_trace(go.Scatter(y=z_robust[batch_id][:, col], mode="lines",
                                 name=f"batch {batch_id}, robust",
                                 line=dict(color=colour, width=1.5),
                                 showlegend=(col == 0)),
                      row=row, col=col + 1)
        fig.add_trace(go.Scatter(y=z[batch_id][:, col], mode="lines",
                                 name=f"batch {batch_id}, mean and sd",
                                 line=dict(color=colour, width=1, dash="dash"),
                                 showlegend=(col == 0)),
                      row=row, col=col + 1)
        fig.add_hrect(y0=-2, y1=2, fillcolor=GREY, opacity=0.15, line_width=0,
                      row=row, col=col + 1)
        for level in (-2, 2):
            fig.add_hline(y=level, line_dash="dot", line_color=GREY,
                          row=row, col=col + 1)
title = "Distance from the other batches: robust (solid) and mean-and-sd (dashed)"
fig.update_layout(title=title, height=420)
fig.show()

Signed distance of batch 37 (top, aqua) and batch 34 (bottom, orange) from the other batches for each of the six tags. Solid: the robust distance, the median and 1.4826 MAD of the others, EWMA-smoothed. Dashed: the distance from the mean of the others in units of their standard deviation. The band between plus and minus two is shaded.

Batch 37 leaves the band in the conversion and latex density within the first 20 samples and stays out. Batch 34 leaves it midway, first in the two service temperatures and the energy released, then some 20 samples later in the conversion and latex density. The robust distance moves batch 34’s onsets by two to five samples, and for batch 37 it also puts the cooling-water and jacket temperatures outside the band from the first sample, where the standard-deviation distance never holds them there long enough to count. This demonstrates the added value of robust methods.

7.9.6. One fault, two places in the score plot

The same kind of fault appears in two places of the score plot because it started at two different times. A batch model describes deviations in (tag, time) cells, so the time of an event is part of its signature: a slow reaction from the start deviates along \(\mathbf{w}_1\), one from the middle of the batch along \(\mathbf{w}_2\). The size of the fault sets how far along that direction the batch lies, and the two differ there too: 30% more impurity in batch 37 against 50% in batch 34.

That makes batch models useful for diagnosis, and is also a caution. A library of known faults built in score space needs the onset time as a coordinate, since a fault seen at one onset appears as a new fault at another. The on-line section returns to both batches with a model that has seen neither fault, where the two places become two statistics.

7.9.7. Predicted quality

A PLS model predicts the quality attributes as well, and the predicted quality of the two faulty batches, from a model fitted without each of them, shows what the trajectories knew before the laboratory did.

Python, 72 lines

# RMSEE measures the fitted values. RMSEP is the same error when each batch in turn
# is left out of the fit and predicted by a model that never saw it: the loop takes
# about two minutes. Both component counts are fitted in the same pass, to see what
# the second one buys.
rmsee = np.sqrt(((quality - model.predictions_) ** 2).mean())
held_out, held_out_one = {}, {}
for batch_id in trajectories:
    rest = {b: t for b, t in trajectories.items() if b != batch_id}
    left_out, target = {batch_id: trajectories[batch_id]}, quality.loc[list(rest)]
    for components, store in ((2, held_out), (1, held_out_one)):
        model_without = BatchPLS(n_components=components).fit(rest, target)
        store[batch_id] = model_without.predict(left_out).y_hat.loc[batch_id]
held_out = pd.DataFrame(held_out).T
rmsep = np.sqrt(((held_out - quality) ** 2).mean())
rmsep_one = np.sqrt(((pd.DataFrame(held_out_one).T - quality) ** 2).mean())
sd = quality.std(ddof=1)
for variable in ("Composition", "ParticleSize"):
    fitted = rmsee[variable]
    print(f"{variable}: RMSEE {fitted:.3g} ({fitted / sd[variable]:.2f} sd),",
          f"RMSEP {rmsep[variable]:.3g} ({rmsep[variable] / sd[variable]:.2f} sd)")
# Composition: RMSEE 0.00106 (0.71 sd), RMSEP 0.00122 (0.81 sd)
# ParticleSize: RMSEE 1.87 (0.60 sd), RMSEP 2.42 (0.78 sd)
print("RMSEP / sd, one component: ", (rmsep_one / sd).round(2).to_dict())
# RMSEP / sd, one component:  {'Composition': 0.8, 'ParticleSize': 0.87,
# 'Branching': 0.28, 'CrossLinking': 0.28, 'Polydispersity': 0.67}
print("RMSEP / sd, two components:", (rmsep / sd).round(2).to_dict())
# RMSEP / sd, two components: {'Composition': 0.81, 'ParticleSize': 0.78,
# 'Branching': 0.28, 'CrossLinking': 0.28, 'Polydispersity': 0.73}
# The RMSEP says by how far the model misses; the cross-validated R2 says how much
# of the attribute it predicts. Five folds over whole batches, which takes about 40
# seconds.
unfolded = pd.DataFrame({b: t.to_numpy().ravel(order="F")
                         for b, t in trajectories.items()}).T
q2 = PLS.select_n_components(unfolded, quality.loc[unfolded.index], max_components=2,
                             cv=5, random_state=0).r2y_validated.loc[2]
print("Q2 per attribute:", q2[quality.columns].round(3).to_dict())
# Q2 per attribute: {'Composition': 0.322, 'ParticleSize': 0.339, 'Branching': 0.914,
# 'CrossLinking': 0.914, 'Polydispersity': 0.476}
for variable in ("Composition", "ParticleSize"):
    error = rmsep[variable]
    fig = model.predictions_vs_observed_plot(quality, variable=variable)
    fig.add_annotation(xref="paper", yref="paper", x=0.02, y=0.98,
                       showarrow=False, align="left",
                       text=f"RMSEP {error:.3g} ({error / sd[variable]:.2f} sd)"
                            f"<br>Q² {q2[variable]:.3f}")
    lo, hi = quality[variable].min(), quality[variable].max()
    half = 2 * rmsep[variable]
    # Two prediction errors either side of y = x, as a shape rather than a trace so that
    # `layer="below"` keeps it behind the batches it is there to be read against.
    fig.add_shape(type="path", layer="below", line_width=0,
                  fillcolor="rgba(31, 61, 122, 0.10)",
                  path=f"M {lo},{lo - half} L {hi},{hi - half}"
                       f" L {hi},{hi + half} L {lo},{lo + half} Z")
    fig.show()
faulty = [34, 37]
# Where the measured value ranks among the 53 batches, and where the prediction from
# a model that never saw the batch ranks among the 53 predictions. Rank 1 is the
# lowest.
observed_rank, held_out_rank = quality.rank().astype(int), held_out.rank().astype(int)
moves = pd.DataFrame({attribute: [f"{observed_rank.loc[b, attribute]} to"
                                  f" {held_out_rank.loc[b, attribute]}"
                                  for b in faulty]
                      for attribute in quality.columns},
                     index=[f"batch {b}" for b in faulty])
print(moves.to_string())
#          Composition ParticleSize Branching CrossLinking Polydispersity
# batch 34     5 to 26       1 to 1    4 to 2       4 to 2        17 to 2
# batch 37      4 to 1       2 to 2    1 to 1       1 to 1         1 to 1
shift = (model.scores_.loc[faulty] * model.y_loadings_.loc["Composition"]).round(2)
print("shift of the scaled composition per component:", shift.to_dict("index"))
# shift of the scaled composition per component:
# {34: {1: -1.39, 2: 0.99}, 37: {1: -2.89, 2: -0.14}}

Observed against fitted composition (left) and particle size (right) for the 53 batches, with batches 34 (orange) and 37 (aqua) marked. Each panel lists two errors: RMSEE, the scatter of these fitted values about the \(y = x\) line, and RMSEP, the same scatter when every batch is left out of the fit in turn, each in the attribute’s own units and in standard deviations of it. Beside them is \(Q^2\), the cross-validated \(R^2\) of that attribute: the errors say by how far the model misses, \(Q^2\) how much of the attribute it predicts. The shaded band is two RMSEP either side of the line: the scatter to expect from a batch the model has not seen, and so wider than the scatter of the fitted values drawn here.

The two components explain 72% of the quality block when they are fitted to it. What they predict in a batch they have not seen is the \(Q^2\) beside each error, and the two numbers carry the same information.

Write the prediction error as a fraction of the attribute’s own standard deviation \(s\), so \(f = \text{RMSEP}/s\). Predicting the average of every batch would give \(f = 1\), and a model is worth having only below that. Squaring gives the share of the variance the predictions still miss, \(f^2\), and what remains is the share the model accounts for:

\[Q^2 \approx 1 - f^2\]

Both attributes here miss by about 0.8 standard deviations, so \(f^2\) is about 0.64: that is the share of the variance still missed, leaving \(1 - 0.64 = 0.36\), about a third of each attribute predicted.

The relation is approximate, and writing both quantities over the same sums of squares shows why. Both are built from the prediction error sum of squares, PRESS, against the total sum of squares of the attribute, TSS:

\[Q^2 = 1 - \frac{\text{PRESS}}{\text{TSS}} \qquad\text{and}\qquad f^2 = \left(\frac{\text{RMSEP}}{s}\right)^2 = \frac{\text{PRESS}}{\text{TSS}} \cdot \frac{N-1}{N}\]

The two therefore differ by the factor \((N-1)/N\), which is 52/53 here: the RMSEP divides by \(N\) and the standard deviation \(s\) by \(N-1\). A second, larger difference is deliberate: the RMSEP leaves out one batch at a time while the \(Q^2\) uses five folds.

Composition and particle size are the two attributes the model predicts least well: branching and cross-linking reach 0.91 and polydispersity 0.48.

The second component earns almost nothing in prediction: leaving each batch out and refitting with one component gives the same average error over the five attributes, buying particle size (0.87 down to 0.78) and losing polydispersity (0.67 up to 0.73). It is kept because it carries the fault of batch 34, which a one-component model would push into the residual, and because a score plot needs a second axis.

Both faulty batches were in the model, so their fitted values prove nothing on their own. Leaving each one out and refitting is the test: batch 37 is measured lowest of the 53 on branching, cross-linking and polydispersity, and a model that never saw it predicts it lowest on all three. Either batch would have been flagged from the trajectories alone, before the laboratory results arrived.

Batch 34 is the interesting one, because one attribute does not follow.

Attribute	Measured rank, of 53	Held-out predicted rank
Composition	5	26
Particle size	1	1
Polydispersity	17	2

Where batch 34 sits among the 53, rank 1 the lowest value, measured and as predicted by a model fitted on the other 52. Branching and cross-linking are left out: they are one measurement in two units, so they repeat the particle-size row exactly.

Composition is the exception, and the reason is in the scores rather than in the fault. Batch 34’s two scores move its composition in opposite directions, \(t_1\) by 1.4 standard deviations down and \(t_2\) by 1.0 up, so the prediction lands mid-pack while the measurement sits near the bottom. A component’s share of the quality block says how much it explains across all batches, not how far it moves any one of them.

7.9.8. Predicting quality before the batch ends

The fitted values of the previous section used the whole batch. It is far more interesting if it is possible to predict what its final quality will be while the batch is still running:

• the unfolded row is complete up to the current sample and empty after it, so the scores are estimated from the observed cells alone, the rest treated as missing data (Wold and co-workers, 2009, who set out the projection to the model plane and point to the literature for the other estimators);

• the model’s regression from scores to quality turns those scores into a prediction;

• the prediction error after \(k\) samples, RMSEP, comes from refitting without each batch in turn and tracing the held-out batch (online_rmse, about a minute).

The estimator used here is trimmed score regression (Arteaga and Ferrer, 2002), which García-Muñoz, Kourti and MacGregor (2004) found stable from the first samples of a batch. Applying the loadings of the measured cells to those cells alone gives a trimmed score: the score the row would have if every cell still to come sat at its average.

It is biased, so the same trimming is done to every training batch, where the score from the complete row is known as well, and those true scores are regressed on the trimmed ones. Applying that regression to the new row is the estimate.

Fitting the regression over history is what keeps it steady early. While few cells have been measured they say little about the final scores, the coefficients are small, and the estimate stays near the average batch; as more cells arrive the same regression carries it out to where the batch is.

Python, 35 lines

squared = 0
# leave one batch out: about a minute
for held_out in trajectories:
    rest = {b: t for b, t in trajectories.items() if b != held_out}
    model_wo = BatchPLS(n_components=2).fit(rest, quality.loc[list(rest)])
    squared = squared + model_wo.online_rmse({held_out: trajectories[held_out]},
                                             quality.loc[[held_out]]) ** 2
# after 1, 2, ..., 200 samples
rmsep_k = np.sqrt(squared / len(trajectories))
relative = rmsep_k / quality.std()
at = [10, 25, 50, 150, 200]
# RMSEP / sd averaged over the five attributes
print(relative.loc[at].mean(axis=1).round(2).tolist())
# [0.97, 0.89, 0.8, 0.62, 0.58]
# The sample from which each curve stays below one standard deviation, which is not the
# first sample it dips below: a curve that crosses and comes back has not settled.
stays = {a: int(relative.index[np.flatnonzero((relative[a] >= 1).to_numpy())[-1] + 1])
         if (relative[a] >= 1).any() else 1 for a in relative.columns}
print("stays below the standard deviation from sample:", stays)
# stays below the standard deviation from sample: {'Composition': 1,
# 'ParticleSize': 125, 'Branching': 2, 'CrossLinking': 2, 'Polydispersity': 2}

# one more figure colour
MAGENTA = "#b03a78"
# one per attribute; branching and cross-linking coincide, so they share aqua
COLOURS = (BLUE, ORANGE, AQUA, AQUA, MAGENTA)
fig = go.Figure()
for attribute, colour in zip(quality.columns, COLOURS):
    fig.add_trace(go.Scatter(x=relative.index[9:], y=relative[attribute].iloc[9:],
                             name=attribute, line=dict(color=colour)))
fig.add_hline(y=1.0, line_color=GREY, annotation_text="as good as the average batch")
fig.update_layout(xaxis_title="Samples observed",
                  yaxis_title="RMSEP / standard deviation of the attribute",
                  height=420)
fig.show()

Root-mean-square error of the leave-one-batch-out prediction after \(k\) samples, divided by the standard deviation of each attribute. Branching and cross-linking coincide. A curve above 1 is worse than predicting the average batch.

No curve begins far above 1, because the estimator does not extrapolate from the few cells it has: with little observed it falls back on the average batch, whose error is one standard deviation by construction. What the curves show is therefore when each attribute becomes predictable, not how badly it starts.

The particle size is the one attribute that does not settle early. Its error hovers at about one standard deviation through the first half of the batch and stays below it only from sample 125, where the fit of every (tag, time) cell also placed the information.

Nomikos and MacGregor (1995) give the reason: particle size is set largely by the seed particles charged before the batch starts, which the trajectories do not record. It was the attribute their model explained least.

Kourti, Nomikos and MacGregor (1995) set out how such a block of before-the-batch measurements enters the model beside the trajectories. The third case study does that with two of them.

Two batches show what the curves summarise, both for the particle size: batch 4, nearest the average quality, and batch 34, whose particle size is the lowest of the 53.

Each prediction is shown against the number of samples observed, the prediction from the complete batch dashed and the measured value solid. The band is one prediction error at that sample, read off the RMSEP curve; it is not a prediction interval. Unlike that curve, these two traces are of batches the model was fitted to.

Python, 45 lines

print("average particle size over the 53 batches:",
      round(quality["ParticleSize"].mean(), 1))
# average particle size over the 53 batches: 1257.1
band = rmsep_k["ParticleSize"]
for batch_id, colour in ((4, BLUE), (34, ORANGE)):
    trace = model.predict_online_trace(trajectories[batch_id])
    evolving = trace.y_hat["ParticleSize"]
    fig = go.Figure()
    fig.add_trace(go.Scatter(x=trace.time[4:], y=(evolving + band).iloc[4:],
                             line=dict(width=0), showlegend=False))
    fig.add_trace(go.Scatter(x=trace.time[4:], y=(evolving - band).iloc[4:],
                             fill="tonexty", fillcolor="rgba(200, 200, 200, 0.35)",
                             line=dict(width=0), name="one prediction error"))
    fig.add_trace(go.Scatter(x=trace.time[4:], y=evolving.iloc[4:],
                             name="prediction so far", line=dict(color=colour)))
    fig.add_hline(y=model.predictions_.loc[batch_id, "ParticleSize"],
                  line_dash="dash", line_color=GREY,
                  annotation_text="final prediction")
    fig.add_hline(y=quality.loc[batch_id, "ParticleSize"], line_color="black",
                  annotation_text="measured")
    # The printed figure puts the two panels on one axis; pinning the range here
    # keeps this plot comparable with it, and the two batches with each other.
    fig.update_layout(title=f"Batch {batch_id}: ParticleSize",
                      xaxis_title="Samples observed", yaxis_range=[1242, 1261],
                      height=380)
    fig.show()
    steps = evolving.loc[[10, 50, 100, 150]].round(1).tolist()
    print(f"batch {batch_id}: measured {quality.loc[batch_id, 'ParticleSize']:.1f},",
          f"final prediction {model.predictions_.loc[batch_id, 'ParticleSize']:.1f},",
          f"after 10/50/100/150 samples {steps}")
    # batch 4: measured 1256.9, final prediction 1257.1, after 10/50/100/150 samples
    # [1256.8, 1256.5, 1256.4, 1257.4]
    # batch 34: measured 1243.7, final prediction 1245.3, after 10/50/100/150 samples
    # [1257.2, 1255.3, 1254.9, 1245.3]

print(f"band: widest {band.max():.2f} at sample {band.idxmax()},",
      f"narrowest {band.iloc[-1]:.2f} at the end")
# band: widest 3.45 at sample 103, narrowest 2.42 at the end
average = quality["ParticleSize"].mean()
moved = (model.predict_online_trace(trajectories[34]).y_hat["ParticleSize"]
         - average).abs()
outside = moved > band
print("batch 34 first sits further from the average than the band at sample",
      int(outside.idxmax()))
# batch 34 first sits further from the average than the band at sample 107

The predicted particle size as each batch is observed: batch 4 (left), nearest the average quality, and batch 34 (right), the lowest particle size of the 53. Dotted: the average of the 53 batches. Dashed: the prediction from the complete batch. Solid: the measured value. Band: one prediction error at that sample, the same band in both panels, centred on each batch’s own prediction. Both predictions start on the dotted rule, because that is where the estimator sits at the start of the batch when no data is available. Gradually, as data comes in, the batch prediction changes. On the left the three rules lie within 0.2 units of each other and draw as one.

The band is not a property of either batch. It is the leave-one-batch-out prediction error after that many samples, so it says how far a prediction made at that point in the batch has missed over all 53, and the same band is drawn in both panels; only its centre differs.

The band is widest in the middle of the batch, at sample 103, and narrowest at the end, 2.42 units.

While a prediction sits inside the band it says only that the batch looks average. Batch 34’s leaves the band at sample 107, seven samples after the impurity enters. The SPE chart flagged the same batch at 105.

Batch 4 is already at the average, so its prediction is close to right from the tenth sample and moves little after it. Batch 34 starts in the same place and leaves it between samples 100 and 150, arriving 12 units lower, which is where the RMSEP curve for the particle size settles below one standard deviation.

The two panels are the same estimator on the same model. What separates them is when the batch’s own cells began to say something the average did not. A prediction that has not moved off the average is not yet a statement about the batch.

7.9.9. Would the model have caught it on-line?

The model on this page holds batches 34 and 37, which is right for diagnosing them after the fact and wrong for monitoring. A monitoring model must describe normal operation, as the first case study does when it removes its outliers. The reference model here is fitted to the 51 batches without 34 and 37, and the two are run through it sample by sample, as if they were new batches on a running plant.

Two statistics are tracked. Hotelling’s \(T^2\) of the score estimate says how far the batch so far sits along the model’s components; the SPE of the newest sample says how far that sample sits off them.

The \(T^2\) limit is the same at every sample, depending only on the number of components and of reference batches, and it assumes the score estimates of the reference batches at each sample are normally distributed. What changes is the spread of the score estimates. A score estimated from the first few samples uses only a few cells of the unfolded row, and across the reference batches such estimates scatter far more widely than the final scores.

\(T^2\) at each sample is therefore computed against the covariance of the reference estimates at that same sample, which the monitor stores. Nomikos and MacGregor, in their Technometrics paper of 1995, set their score-chart limits from that spread and note that \(T^2\) needs the covariance at each sample too; García-Muñoz, Kourti and MacGregor (2004) compute it. The SPE limit is fitted sample by sample the same way: at each sample it is a chi-square approximation matched to the mean and the variance of the 51 reference values there.

The same normalisation pins the reference batches’ average. Measuring a set of batches against their own covariance forces their mean \(T^2\) to \(A(N-1)/N\), 1.96 here for two components and 51 batches, whatever the data: it is a property of the arithmetic, not of the process. The grey line in the monitoring figure below is flat for that reason, and it is what lets a single limit serve every sample. The SPE is not normalised this way, so its reference mean does vary.

Python, 18 lines

normal = {b: t for b, t in trajectories.items() if b not in (34, 37)}
reference = BatchPLS(n_components=2).fit(normal, quality.loc[list(normal)])
monitor = BatchMonitor(reference, conf_level=0.99,
                       spe_statistic="instantaneous").fit(normal)
# sd of the estimates, per sample
spread = np.sqrt(np.diagonal(monitor.score_covariance_over_time_, axis1=1, axis2=2))
# relative to the final scores
spread = spread / reference.scores_.std(ddof=1).to_numpy()
fig = go.Figure()
for a, colour in enumerate((BLUE, ORANGE)):
    fig.add_trace(go.Scatter(x=np.arange(1, len(spread) + 1), y=spread[:, a],
                             name=f"t{a + 1}", line=dict(color=colour)))
fig.add_hline(y=1.0, line_dash="dash", line_color=GREY,
              annotation_text="spread of the final scores")
fig.update_layout(xaxis_title="Samples observed",
                  yaxis_title="Spread relative to the final scores",
                  yaxis_type="log", height=380)
fig.show()

The spread of the on-line score estimates across the 51 reference batches, divided by the spread of their final scores, against the number of samples observed (logarithmic axis). The estimator pulls the scores toward the average batch while little has been observed, so early estimates lie closer together than the final scores do, and the spread grows to match them as the batch fills in. \(T^2\) at each sample is scaled by this spread, which is what lets the limit stay the same throughout the batch.

The limits are set at 99%, and an alarm here means three consecutive samples above the limit, the same kind of rule the departure analysis used.

Python, 77 lines

def first_sustained(alarm, run=3):
    """First of the `run` consecutive samples above the limit, counting from one,
    or None."""
    runs = np.convolve(alarm.astype(int), np.ones(run, dtype=int), mode="valid") == run
    return int(runs.argmax()) + 1 if runs.any() else None

print(f"T2 limit at 99%: {monitor.t2_limit_over_time_[0]:.1f}")
# T2 limit at 99%: 10.5
# A(N-1)/N at every sample, whatever the data do
mean_t2 = np.asarray(monitor.t2_mean_over_time_)
print(f"reference mean T2: {mean_t2.min():.2f} to {mean_t2.max():.2f}")
# reference mean T2: 1.96 to 1.96
for batch_id in (37, 34, 4):
    result = monitor.monitor(trajectories[batch_id])
    t2_alarm_at = first_sustained(result.t2_alarm)
    print(f"batch {batch_id}: T2 alarm after {t2_alarm_at} samples;",
          f"SPE alarm after {first_sustained(result.spe_alarm)} samples")
    # batch 37: T2 alarm after 23 samples; SPE alarm after 15 samples
    # batch 34: T2 alarm after 190 samples; SPE alarm after 105 samples
    # batch 4: T2 alarm after None samples; SPE alarm after None samples

def longest_run(alarm):
    """Longest stretch of consecutive samples above the limit."""
    best = current = 0
    for above in alarm:
        current = current + 1 if above else 0
        best = max(best, current)
    return best

traced = [monitor.monitor(t) for t in normal.values()]
t2_rate = np.mean([r.t2_alarm.mean() for r in traced])
spe_rate = np.mean([r.spe_alarm.mean() for r in traced])
print(f"reference batches: {t2_rate:.2%} of T2 values and",
      f"{spe_rate:.2%} of SPE values above their limits")
# reference batches: 0.17% of T2 values and 1.13% of SPE values above their limits
for run in (3, 5, 10):
    print(f"  with an SPE alarm of {run} consecutive samples somewhere:",
          sum(first_sustained(r.spe_alarm, run) is not None for r in traced), "of 51;",
          "T2:", sum(first_sustained(r.t2_alarm, run) is not None for r in traced),
          "of 51")
    #   with an SPE alarm of 3 consecutive samples somewhere: 13 of 51; T2: 1 of 51
    #   with an SPE alarm of 5 consecutive samples somewhere: 8 of 51; T2: 1 of 51
    #   with an SPE alarm of 10 consecutive samples somewhere: 3 of 51; T2: 1 of 51
print("  longest SPE alarm run in a reference batch:",
      max(longest_run(r.spe_alarm) for r in traced), "samples")
#   longest SPE alarm run in a reference batch: 15 samples
print("batch 34: above the SPE limit for",
      int(monitor.monitor(trajectories[34]).spe_alarm[104:].sum()),
      "of its last 96 samples; batch 37: longest T2 alarm run",
      longest_run(monitor.monitor(trajectories[37]).t2_alarm), "samples")
# batch 34: above the SPE limit for 87 of its last 96 samples;
# batch 37: longest T2 alarm run 178 samples
online_monitoring_plot(monitor, trajectories[37], "t2").show()
fig = online_monitoring_plot(monitor, trajectories[34], "spe")
fig.add_vline(x=100, line_dash="dash", line_color=ORANGE,
              annotation_text="impurity enters")
fig.show()
alarm_k = first_sustained(monitor.monitor(trajectories[34]).spe_alarm)
for k in (alarm_k, alarm_k + 4):
    residuals = reference.predict_online(trajectories[34], upto_k=k).residuals
    shares = residuals.xs(k - 1, level="sequence") ** 2
    shares = (shares / shares.sum() * 100).round(0).astype(int)
    print(f"batch 34 after {k} samples, share of the squared residual per tag [%]:",
          shares.to_dict())
    # batch 34 after 105 samples, share of the squared residual per tag [%]:
    # {'Conversion': 1, 'CoolingTemp': 31, 'EnergyReleased': 11, 'JacketTemp': 17,
    # 'LatexDensity': 0, 'ReactorTemp': 41}
    # batch 34 after 109 samples, share of the squared residual per tag [%]:
    # {'Conversion': 2, 'CoolingTemp': 36, 'EnergyReleased': 19, 'JacketTemp': 30,
    # 'LatexDensity': 1, 'ReactorTemp': 12}
alarm_residuals = reference.predict_online(trajectories[34], upto_k=alarm_k).residuals
at_alarm = alarm_residuals.xs(alarm_k - 1, level="sequence") ** 2
at_alarm = (at_alarm / at_alarm.sum() * 100).round(0).astype(int)
fig = go.Figure(go.Bar(x=at_alarm.index, y=at_alarm.values, marker_color=BLUE))
fig.update_layout(title=f"Batch 34 after {alarm_k} samples",
                  yaxis_title="Share of the squared residual [%]", height=320)
fig.show()

On-line monitoring of the two faulty batches against the reference model of 51 normal batches. Left: Hotelling’s \(T^2\) of batch 37 (aqua) with the 99% limit (dashed) and the reference-batch mean (grey, flat by construction); the first sustained alarm is after 23 samples. Middle: the SPE of the newest sample of batch 34 (orange) with its per-sample limit; the impurity enters at sample 100 (dotted vertical) and the first sustained alarm is after 105 samples. Right: the share of the squared residual per tag at that alarm sample.

Batch 37 raises both alarms and batch 34 only the SPE, each with most of the batch still to run.

The SPE limit does not hold on this reference set: thirteen of the fifty-one normal batches raise a sustained alarm on it, against one on \(T^2\). Both statistics are autocorrelated, so crossings arrive in runs rather than singly.

Three ways to make the limit harder to trip by chance. Pooling it over five neighbouring samples does not help, leaving all thirteen. The cumulative SPE leaves two and a 99.9% limit none, and both pay for it in the true alarm: batch 34 is caught after 112 and 106 samples rather than 105. That zero is not a false-alarm rate, because a 99.9% quantile fitted from 51 values lies beyond the largest of them.

Python, 25 lines

pooled = BatchMonitor(reference, conf_level=0.99, spe_statistic="instantaneous",
                      spe_window=2).fit(normal)
pooled_traced = [pooled.monitor(t) for t in normal.values()]
pooled_rate = np.mean([r.spe_alarm.mean() for r in pooled_traced])
pooled_count = sum(first_sustained(r.spe_alarm) is not None for r in pooled_traced)
pooled_at = first_sustained(pooled.monitor(trajectories[34]).spe_alarm)
print(f"limit pooled over five samples: {pooled_rate:.2%} of the reference",
      f"SPE values above it; {pooled_count} of 51 batches",
      f"with a three-sample run; batch 34 flagged after {pooled_at} samples")
# limit pooled over five samples: 1.16% of the reference SPE values above it;
# 13 of 51 batches with a three-sample run; batch 34 flagged after 105 samples
# the SPE over every sample observed so far
cumulative = BatchMonitor(reference, conf_level=0.99).fit(normal)
# reference batches with an alarm; batch 34
print([sum(first_sustained(cumulative.monitor(t).spe_alarm) is not None
           for t in normal.values()),
       first_sustained(cumulative.monitor(trajectories[34]).spe_alarm)])
# [2, 112]
tight = BatchMonitor(reference, conf_level=0.999,
                     spe_statistic="instantaneous").fit(normal)   # the 99.9% limit
# reference batches with an alarm; batch 34
print([sum(first_sustained(tight.monitor(t).spe_alarm) is not None
           for t in normal.values()),
       first_sustained(tight.monitor(trajectories[34]).spe_alarm)])
# [0, 106]

Which statistic catches a fault says what kind of fault it is; how much of the batch has been seen says when. A fault along a direction the reference batches already vary in is carried by the score, a large \(T^2\), which is batch 37 from sample 23. A fault the model has no component for leaves the newest samples not fitting, which is the SPE, and that is batch 34 from sample 105.

Early in a batch every fault looks like the second kind, because the model cannot yet place the batch along its components. That is why batch 37’s SPE crosses at 15, before its \(T^2\) at 23.

The contribution shares name which tags stopped fitting, not what happened in the reactor: the departure analysis puts the onsets in the two service temperatures, while the reactor temperature is simply the tag the model explains least.

The same score estimate that predicts the quality also predicts the rest of the trajectories, as the model’s reconstruction \(\hat{\boldsymbol{\tau}} \mathbf{P}^{T}\) read off for the samples not yet seen (Wold and co-workers, 2009, Eq. 4). García-Muñoz, Kourti and MacGregor (2004) show it to be an adaptive time-series forecast, built on how the variables co-vary over the whole batch rather than on the trend so far.

It is shown in the z form of the departure analysis, each tag as a distance from the 51 normal batches at that sample in their standard deviations, so a departure reads directly against the zero line.

Python, 60 lines

def z_form(frame):
    """Each tag as a distance from the 51 normal batches at that sample, in their
    standard deviations."""
    return pd.DataFrame((frame.to_numpy() - others.mean(axis=0))
                        / others.std(axis=0, ddof=1),
                        columns=frame.columns, index=frame.index)

def forecast_panel(batch_id, tag, from_samples, colour):
    """The normal batches as a band, what the batch did, and the forecast of the rest
    from two points."""
    # A band, not 51 lines: on a noisy tag the lines fill the panel and the forecasts
    # have to be read through them. It covers the middle 90% of the normal batches at
    # each sample.
    spread = np.stack([z_form(t)[tag].to_numpy() for t in normal.values()])
    lo, hi = np.percentile(spread, 5, axis=0), np.percentile(spread, 95, axis=0)
    fig = go.Figure([go.Scatter(x=[*range(len(lo)), *reversed(range(len(hi)))],
                                y=[*lo, *reversed(hi)], fill="toself",
                                fillcolor="rgba(200, 200, 200, 0.45)", line_width=0,
                                name="normal batches, middle 90%")])
    fig.update_layout(title=tag, xaxis_title="Sample [aligned time]", height=320)
    z_actual = z_form(trajectories[batch_id])[tag]
    fig.add_trace(go.Scatter(y=z_actual, mode="lines",
                             name=f"batch {batch_id}, what happened",
                             line=dict(color=colour, width=1), opacity=0.4))
    # the later one apart
    for k, dash, line_colour, width in zip(from_samples, ("dash", "dot"),
                                           (colour, BLUE), (2, 3)):
        forecast = z_form(reference.predict_online(trajectories[batch_id],
                                                   upto_k=k).forecast)[tag]
        fig.add_trace(go.Scatter(x=forecast.index[k:], y=forecast.iloc[k:],
                                 mode="lines", name=f"forecast from sample {k}",
                                 line=dict(color=line_colour, width=width, dash=dash)))
        # the jump from the data used
        fig.add_trace(go.Scatter(x=[k - 1, k - 1],
                                 y=[z_actual.iloc[k - 1], forecast.iloc[k]],
                                 mode="lines",
                                 line=dict(color=line_colour, width=1.5),
                                 showlegend=False))
    fig.add_trace(go.Scatter(y=z_actual.iloc[:from_samples[0]], mode="lines",
                             name=f"batch {batch_id}, observed",
                             line=dict(color=colour, width=3)))
    fig.add_hline(y=0, line_color=GREY)
    fig.update_layout(title=f"Batch {batch_id}: {tag}",
                      yaxis_title="Distance from the normal batches [sd]")
    return fig

forecast_panel(37, "Conversion", (30, 60), AQUA).show()
panel = forecast_panel(34, "CoolingTemp", (60, 115), ORANGE)
panel.add_vline(x=100, line_dash="dash", line_color=ORANGE).show()
for batch_id, tag, k in ((37, "Conversion", 30), (37, "Conversion", 60),
                         (34, "CoolingTemp", 60), (34, "CoolingTemp", 115)):
    forecast = z_form(reference.predict_online(trajectories[batch_id],
                                               upto_k=k).forecast)[tag].iloc[k:]
    print(f"batch {batch_id}, {tag}, from sample {k} onwards:"
          f" forecast mean {forecast.mean():.2f} sd,",
          f"actual {z_form(trajectories[batch_id])[tag].iloc[k:].mean():.2f} sd")
# batch 37, Conversion, from sample 30 onwards: forecast mean -1.81 sd, actual -4.81 sd
# batch 37, Conversion, from sample 60 onwards: forecast mean -3.62 sd, actual -4.82 sd
# batch 34, CoolingTemp, from sample 60 onwards: forecast mean 0.01 sd, actual 2.36 sd
# batch 34, CoolingTemp, from sample 115 onwards: forecast mean 0.36 sd, actual 3.30 sd

Forecast of the rest of the batch from the score estimate, in z form: each tag as a distance from the normal batches at that sample, in their standard deviations, so the zero line is the average batch and the grey band holds the middle 90% of them. Left: the conversion of batch 37 (aqua), observed for 30 samples, and the forecasts from sample 30 onwards (dashed aqua) and from sample 60 onwards (dotted dark blue); what batch 37 did is the faint line. Right: the cooling-water temperature of batch 34 (orange) with the forecasts from sample 60 onwards (dashed orange) and from sample 115 onwards (dotted dark blue); the impurity enters at sample 100. A vertical line joins each forecast to the observed value at the sample it was made from: the forecast uses the batch’s own data up to that sample, then continues along the model’s components. Batch 37’s two forecasts lie below zero and above what happened; the later one, with more of the batch behind it, lies closer to it.

A fault along the model’s components can be forecast; one off them can only be flagged. Batch 37’s two forecasts run in the right direction, the later one closer, while batch 34’s stay on the average batch.

Every number here comes from the reference set, easy to choose because the simulation says which batches are faulty. On plant data those batches are the first thing to get right, and the first case study shows how much of the work that is. The third case study turns to a process where the trajectories are one of four blocks about a batch, and asks which block matters most.

7.9.10. References and readings

The full list of readings on batch data is on the batch process monitoring page; this page lists what it draws on.

• Paul Nomikos, Statistical process control of batch processes, Ph.D thesis, McMaster University, 1995. The source of the simulation and of the two faulty batches.

• Paul Nomikos and John F. MacGregor, “Monitoring batch processes using multiway principal component analysis”, AIChE Journal, 40, 1361-1375, 1994. Describes the simulation and the two faulty batches, and monitors them with a batch PCA model of all nine trajectories.

• Paul Nomikos and John F. MacGregor, “Multi-way partial least squares in monitoring batch processes”, Chemometrics and Intelligent Laboratory Systems, 30, 97-108, 1995.

• Theodora Kourti, Paul Nomikos and John F. MacGregor, “Analysis, monitoring and fault diagnosis of batch processes using multiblock and multiway PLS”, Journal of Process Control, 5, 277-284, 1995.

• Paul Nomikos and John F. MacGregor, “Multivariate SPC charts for monitoring batch processes”, Technometrics, 37, 41-59, 1995. The on-line monitoring scheme, with the statistics compared with their limits at every sample.

• Paul Nomikos, “Detection and diagnosis of abnormal batch operations based on multi-way principal component analysis”, ISA Transactions, 35, 259-266, 1996. The 99.9% limit for the residual chart.

• Henk-Jan Ramaker, Eric N. M. Van Sprang, Johan A. Westerhuis, Stephen P. Gurden, Age K. Smilde and Frank H. Van Der Meulen, “Performance assessment and improvement of control charts for statistical batch process monitoring”, Statistica Neerlandica, 60, 339-360, 2006. Assesses a batch chart by its overall type I error and by the time it takes to signal.

• Francisco Arteaga and Alberto Ferrer, “Dealing with missing data in MSPC: several methods, different interpretations, some examples”, Journal of Chemometrics, 16, 408-418, 2002. Trimmed score regression.

• Salvador García-Muñoz, Theodora Kourti and John F. MacGregor, “Model predictive monitoring for batch processes”, Industrial and Engineering Chemistry Research, 43, 5929-5941, 2004. Compares the estimators of the scores of a batch so far, shows that the forecast of the rest of the batch is an adaptive time-series forecast, and computes Hotelling’s \(T^2\) against the time-varying covariance of the score estimates.

• Svante Wold, Nouna Kettaneh-Wold, John F. MacGregor and Kevin G. Dunn, “Batch process modeling and MSPC”, Comprehensive Chemometrics, 2, chapter 2.10, 163-197, 2009. Mid-batch prediction of quality and of the remaining trajectories.