Batch PCA outlier diagnosis: the DuPont polymerization reactor#

Case study for issue #155.

Industrial nylon is made in a batch autoclave. The critical quality of a batch is measured in the laboratory about twelve hours after the batch has ended, so there is no feedback that could correct a running batch, and a long hold-up before its disposition is known. What the plant does have is the trajectory of ten process measurements (temperatures, pressures and flows) over every batch. The question of this case study is what those trajectories can tell about a batch, and what they cannot.

The data were supplied by DuPont and are the worked example of Nomikos and MacGregor (1995). Batches 40, 41, 42, 50, 51, 53, 54 and 55 had a final quality well outside the acceptable limit; batches 38, 45, 46, 49 and 52 were above or very close to it.

The complete script is dupont_batch_pca.py in this directory:

uv run python docs/user_guide/case_studies/batch/dupont_batch_pca.py --output-dir case-study-output/dupont

It prints the numbers quoted below and writes every figure as an HTML file to the output directory. The sections below quote it piece by piece.

Data#

Industrial batch polymerization: 55 batches, each already aligned to 100 equal time intervals, with ten tags per interval. The values are scaled for confidentiality. process_improve.batch.load_dupont() downloads the file and returns the standard batch dictionary, one 100-by-10 frame per batch.

CONF_LEVEL = 0.95
SPE_OUTLIER = 49
SCORE_OUTLIERS = [50, 51, 52, 53, 54, 55]
DIFFERENT_BUT_ACCEPTABLE = [37, 39, 43, 44, 45, 46, 47, 48]
POOR_QUALITY_NOT_VISIBLE = [38, 40, 41, 42]
HIGHLIGHT = '{"color": "red", "width": 4}'  # Plotly line style, JSON-encoded, for the highlighted batches
LABELS = {"show_labels": True}

def load_data(url: str | None = None) -> dict:
    """Download the 55 aligned batches; each is a 100-sample by 10-tag frame."""
    batches = load_dupont(url=url)
    first = next(iter(batches.values()))
    print(f"{len(batches)} batches x {first.shape[0]} samples x {first.shape[1]} tags")
    return batches

Looking at one tag across all batches is the first check. The trajectories overlay well (the alignment has already been done), and a few batches are visibly unusual, but a plot per tag cannot rank 55 batches on ten variables at once.

def plot_raw(batches: dict, tag: str, highlight: list[int]) -> go.Figure:
    """Overlay one tag for every batch, with the batches of interest drawn in red."""
    return plot_all_batches_per_tag(
        batches, tag, batches_to_highlight={HIGHLIGHT: highlight}, extra_info=f"highlighting {highlight}"
    )

Model A: batch PCA on all 55 batches#

process_improve.batch.BatchPCA unfolds the batches batchwise: each batch becomes one row of 10 tags times 100 samples, 1000 columns. Centring the columns removes the average trajectory, and scaling them to unit variance gives every (tag, time) cell the same weight, so the two components describe how the batches deviate from the average batch.

def fit_model_a(batches: dict) -> BatchPCA:
    """Two-component batch PCA on all 55 batches.

    Every batch becomes one row of 10 tags x 100 samples = 1000 columns; the
    columns are centred (removing the average trajectory) and scaled to unit
    variance, so the model describes batch-to-batch deviations.
    """
    model = BatchPCA(n_components=2).fit(batches)
    per_component = model.r2_per_component_.to_numpy()
    print(
        f"Model A: R2 per component = {per_component[0]:.3f}, {per_component[1]:.3f}; cumulative = {per_component.sum():.3f}"
    )
    print(
        f"Model A: largest |t1| = {abs_top(model.scores_.iloc[:, 0], 4)}; largest |t2| = {abs_top(model.scores_.iloc[:, 1], 3)}"
    )
    spe = model.spe_.iloc[:, -1]
    print(
        f"Model A: largest SPE = batch {spe.idxmax()} ({spe.max():.1f} vs {CONF_LEVEL:.0%} limit {model.spe_limit(conf_level=CONF_LEVEL):.1f})"
    )
    return model


def abs_top(values: pd.Series, n: int) -> list:
    """Batch identifiers with the largest absolute values, largest first."""
    return values.abs().nlargest(n).index.tolist()



Model A: R2 per component = 0.383, 0.176; cumulative = 0.559
Model A: largest |t1| = [54, 52, 50, 51]; largest |t2| = [53, 55, 50]
Model A: largest SPE = batch 49 (39.3 vs 95% limit 29.1)

Two components explain 55.9% of the variance in the unfolded matrix. The score plot shows batches 50 to 55 far from the rest; they pull the model towards themselves, which is the first sign that the model needs to be rebuilt without them. The SPE plot flags a different batch, 49, which sits in the score plot among the normal batches: its problem is not a large deviation along the main directions of variation but a break in the correlation structure.

Batch 49: which variables, and when#

The raw data are ambiguous about batch 49. Flow-1 looks suspicious in the overlay, but it is a noisy tag in every batch. The SPE contributions settle the question. The signed contributions are the residuals of every (tag, time) cell after the two-component reconstruction; their squares add up to the SPE, so the squares are each cell’s share of it. Summing the shares per tag ranks the variables, and summing them per time sample locates the event.

def diagnose_spe_outlier(
    model: BatchPCA, batches: dict, batch_id: int = SPE_OUTLIER
) -> tuple[pd.DataFrame, pd.Series, pd.Series]:
    """Which variables, and when, push one batch off the model plane.

    The signed SPE contributions are the residuals of every (tag, time) cell;
    their squares add up to the batch's SPE, so the squares are each cell's
    share of it. Summing the shares per tag ranks the variables; summing them
    per time sample locates the event.
    """
    spe_share = model.spe_contributions(model.unfold_and_scale(batches)) ** 2
    by_tag = spe_share.loc[batch_id].groupby(level="tag", sort=False).sum()
    by_time = spe_share.loc[batch_id].groupby(level="sequence").sum()
    print(
        f"Batch {batch_id}: SPE share per tag = "
        + ", ".join(f"{tag} {share:.0%}" for tag, share in (by_tag / by_tag.sum()).items())
    )
    peak = by_time.nlargest(7).index
    print(f"Batch {batch_id}: the seven largest per-sample shares sit at samples {sorted(int(k) for k in peak)}")
    return spe_share, by_tag, by_time


def plot_share_over_time(by_time: pd.Series, batch_id: int) -> go.Figure:
    """Bar chart of a batch's SPE share per time sample."""
    fig = go.Figure(go.Bar(x=list(by_time.index), y=by_time.to_numpy()))
    fig.update_layout(
        title=f"SPE share per time sample, batch {batch_id}",
        xaxis_title="Time [sequence order]",
        yaxis_title="Share of SPE",
    )
    return fig



Batch 49: SPE share per tag = Flow-1 3%, Flow-2 18%, Press-1 5%, Press-2 15%, Press-3 12%,
          TempC-1 19%, TempH-1 14%, TempR-1 7%, TempR-2 3%, TempR-3 4%
Batch 49: the seven largest per-sample shares sit at samples [57, 58, 59, 60, 61, 62, 63]

Flow-1 carries 3% of the residual. The residual belongs to the cooling and heating medium temperatures and to the pressures, and it is concentrated in a short window around samples 57 to 63: a small disturbance in the heating, cooling and pressure systems during that stretch of the batch. Nomikos and MacGregor report that the quality of batch 49 was barely acceptable, which is consistent with a short event rather than a batch that was wrong throughout. process_improve.batch.unfolded_contribution_plot() draws the full vector of 1000 bars grouped by tag, and the same data summed per tag.

The score outliers#

Batches 50 to 55 are far out along the components, so the tool here is the score contribution: how much every (tag, time) cell contributes to \(t_1\) or \(t_2\). The loading \(p_1\), drawn as a function of time with process_improve.batch.time_varying_loading_plot(), shows which parts of the batch each component describes.

def diagnose_score_outliers(model: BatchPCA, batches: dict) -> dict[str, pd.DataFrame]:
    """Score contributions for the batches that stand out on t1 and t2."""
    scaled = model.unfold_and_scale(batches)
    contributions = {
        "t1": model.score_contributions(scaled, component=1),
        "t2": model.score_contributions(scaled, component=2),
    }
    for batch_id, component in ((54, "t1"), (55, "t2")):
        by_tag = contributions[component].loc[batch_id].groupby(level="tag", sort=False).sum()
        print(
            f"Batch {batch_id}: {component} contributions per tag = "
            + ", ".join(f"{tag} {value:+.1f}" for tag, value in by_tag.items())
        )
    return contributions



Batch 54: t1 contributions per tag = Flow-1 +5.8, Flow-2 +4.6, Press-1 +5.9, Press-2 +7.4, Press-3 +8.1,
          TempC-1 +7.4, TempH-1 +5.2, TempR-1 +7.5, TempR-2 +8.7, TempR-3 +6.9
Batch 55: t2 contributions per tag = Flow-1 +1.4, Flow-2 +1.7, Press-1 +3.2, Press-2 +6.8, Press-3 +8.5,
          TempC-1 +6.3, TempH-1 +3.8, TempR-1 +0.8, TempR-2 +2.6, TempR-3 +0.6

Batch 54 has a high \(t_1\) because every tag contributes in the same direction: the whole batch ran away from the average trajectory. Batch 55 stands out on \(t_2\) through the pressures and the cooling-medium temperature.

Model B: exclude batches 49 to 55 and rebuild#

A reference model must describe normal operation, so the batches identified so far are removed and the model refitted on the remaining 48 batches, now with three components.

def fit_model_b(batches: dict) -> BatchPCA:
    """Rebuild without batches 49 to 55: a second cluster appears on t2 and t3."""
    kept = {batch_id: batch for batch_id, batch in batches.items() if batch_id < SPE_OUTLIER}
    model = BatchPCA(n_components=3).fit(kept)
    per_component = model.r2_per_component_.to_numpy()
    print(f"Model B ({len(kept)} batches): R2 per component = " + ", ".join(f"{value:.3f}" for value in per_component))
    print(
        f"Model B: largest |t2| = {abs_top(model.scores_.iloc[:, 1], 6)}; largest |t3| = {abs_top(model.scores_.iloc[:, 2], 6)}"
    )
    return model


def diagnose_different_batch(model: BatchPCA, batches: dict, batch_id: int = 39) -> pd.DataFrame:
    """Return the t3 contributions of a batch from the second cluster."""
    kept = {key: value for key, value in batches.items() if key in model.batch_ids_}
    contributions = model.score_contributions(model.unfold_and_scale(kept), component=3)
    by_tag = contributions.loc[batch_id].groupby(level="tag", sort=False).sum()
    print(
        f"Batch {batch_id}: t3 contributions per tag = "
        + ", ".join(f"{tag} {value:+.1f}" for tag, value in by_tag.items())
    )
    return contributions



Model B (48 batches): R2 per component = 0.333, 0.133, 0.085
Model B: largest |t2| = [37, 48, 44, 46, 16, 29]; largest |t3| = [45, 39, 46, 47, 43, 14]
Batch 39: t3 contributions per tag = Flow-1 +0.5, Flow-2 +3.5, Press-1 +0.8, Press-2 +2.3, Press-3 +6.5,
          TempC-1 +4.9, TempH-1 +1.7, TempR-1 -0.1, TempR-2 +0.9, TempR-3 +0.2

With the extreme batches gone, a second group separates on the \(t_2\) and \(t_3\) plane: batches 37, 39 and 43 to 48. A contribution is the weighted difference between two points, and either point can be the average of a group. The group’s mean contribution vector against the model centre, averaged over the eight members, names Press-3 and TempC-1 on both components, with the contribution concentrated in the first 25 samples, so the difference lies in how these batches were started; the \(t_3\) contributions of batch 39 printed above point at the same two tags. Not every member is consistent with the average: TempH-1 takes both signs across the eight. The raw overlay shows that these batches were run on a slightly different pressure profile. Their quality was acceptable. They are not bad batches, they were operated differently, and a model of normal operation should either contain enough of them to describe that mode or leave them out; the course notes leave them out.

Model C: the final model, used to verify the unusual batches#

Model C is fitted on the 40 batches that remain once batch 49, batches 50 to 55 and the eight batches of the second group are removed. The 15 left-out batches are then projected onto it (predict_online at the last sample gives a complete batch’s scores, \(T^2\) and SPE against a model it was not part of) and compared with the 95% limits of the 40 training batches.

def fit_model_c(batches: dict) -> BatchPCA:
    """Rebuild once more without the second cluster: the reference model."""
    excluded = set(range(SPE_OUTLIER, 56)) | set(DIFFERENT_BUT_ACCEPTABLE)
    kept = {batch_id: batch for batch_id, batch in batches.items() if batch_id not in excluded}
    model = BatchPCA(n_components=3).fit(kept)
    per_component = model.r2_per_component_.to_numpy()
    print(f"Model C ({len(kept)} batches): R2 per component = " + ", ".join(f"{value:.3f}" for value in per_component))
    return model


def verify_left_out(model: BatchPCA, batches: dict) -> pd.DataFrame:
    """Project the 15 batches left out of model C onto it and compare them with its 95% limits.

    The on-line projection at the last sample of a complete batch gives its scores, T2 and
    SPE against the centre and scale of the model it was not part of.
    """
    left_out = [SPE_OUTLIER, *SCORE_OUTLIERS, *DIFFERENT_BUT_ACCEPTABLE]
    projected = {b: model.predict_online(batches[b], upto_k=model.n_timesteps_) for b in left_out}
    table = pd.DataFrame(
        {b: (float(r.hotellings_t2), float(r.spe)) for b, r in projected.items()}, index=["T2", "SPE"]
    ).T
    t2_limit, spe_limit = model.hotellings_t2_limit(conf_level=CONF_LEVEL), model.spe_limit(conf_level=CONF_LEVEL)
    print(
        f"Model C: left-out batches above the SPE limit ({spe_limit:.1f}): {sorted(table.index[table['SPE'] > spe_limit])}"
    )
    print(
        f"Model C: left-out batches above the T2 limit ({t2_limit:.2f}): {sorted(table.index[table['T2'] > t2_limit])}"
    )
    return table


def observability_table(model: BatchPCA) -> pd.DataFrame:
    """T2 and SPE of the poor-quality batches that the trajectories do not reveal, against the limits."""
    table = pd.DataFrame(
        {
            "T2": model.hotellings_t2_.loc[POOR_QUALITY_NOT_VISIBLE].iloc[:, -1],
            "T2 limit": model.hotellings_t2_limit(conf_level=CONF_LEVEL),
            "SPE": model.spe_.loc[POOR_QUALITY_NOT_VISIBLE].iloc[:, -1],
            "SPE limit": model.spe_limit(conf_level=CONF_LEVEL),
        }
    )
    inside = bool((table["T2"] < table["T2 limit"]).all() and (table["SPE"] < table["SPE limit"]).all())
    print(f"Model C: batches {POOR_QUALITY_NOT_VISIBLE} inside both limits: {inside}")
    return table



Model C (40 batches): R2 per component = 0.375, 0.114, 0.064
Model C: left-out batches above the SPE limit (24.2): [37, 39, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55]
Model C: left-out batches above the T2 limit (9.27): [37, 50, 51, 52, 53, 54, 55]
Model C: batches [38, 40, 41, 42] inside both limits: True

Every one of the 15 lies above the SPE limit, and the six score outliers and batch 37 above the \(T^2\) limit as well: the model built without them flags them. Batches 38, 40, 41 and 42, known to have had poor quality and kept in the training set, sit inside both limits: nothing in the ten trajectories distinguishes them from the good batches. This is the lesson the case study is built around. A model can only detect what the measurements contain; if the cause of poor quality leaves no trace in the recorded variables, no amount of modelling will find it, and the fix is to measure something else. Nomikos and MacGregor (1995) also left these four batches out of their reference set, on the principle that it should hold only batches with acceptable operation and acceptable product, and built their monitoring model on the remaining 36 batches with three components; they are kept in model C here so that the check can be made. The SBR case-study page runs the same check sample by sample with BatchMonitor.

Running the script#

    """Run the whole case study and write its figures."""
    parser = argparse.ArgumentParser(description=__doc__.split("\n\n")[0])
    parser.add_argument("--output-dir", type=pathlib.Path, default=pathlib.Path("case-study-output/dupont"))
    parser.add_argument("--data-url", default=None, help="override the openmv.net URL (for example a file:// copy)")
    args = parser.parse_args(argv)
    args.output_dir.mkdir(parents=True, exist_ok=True)

    def save(fig: go.Figure, name: str) -> None:
        fig.write_html(args.output_dir / f"{name}.html", include_plotlyjs="cdn")

    batches = load_data(args.data_url)
    save(plot_raw(batches, "TempC-1", SCORE_OUTLIERS), "raw-TempC-1-score-outliers")
    save(plot_raw(batches, "Press-1", SCORE_OUTLIERS), "raw-Press-1-score-outliers")

    model_a = fit_model_a(batches)
    save(model_a.score_plot(settings=LABELS), "model-a-scores")
    save(model_a.spe_plot(settings=LABELS), "model-a-spe")

    save(plot_raw(batches, "Flow-1", [SPE_OUTLIER]), "raw-Flow-1-batch-49")
    save(plot_raw(batches, "TempC-1", [SPE_OUTLIER]), "raw-TempC-1-batch-49")
    spe_share, _by_tag, by_time = diagnose_spe_outlier(model_a, batches)
    save(unfolded_contribution_plot(spe_share, SPE_OUTLIER), "spe-contributions-49")
    save(unfolded_contribution_plot(spe_share, SPE_OUTLIER, by_tag=True), "spe-contributions-49-by-tag")
    save(plot_share_over_time(by_time, SPE_OUTLIER), "spe-contributions-49-by-time")

    save(time_varying_loading_plot(model_a, component=1), "model-a-loadings-p1")
    contributions = diagnose_score_outliers(model_a, batches)
    save(unfolded_contribution_plot(contributions["t1"], 54), "score-contributions-54-t1")
    save(plot_raw(batches, "Press-2", [54]), "raw-Press-2-batch-54")
    save(unfolded_contribution_plot(contributions["t2"], 55), "score-contributions-55-t2")

    model_b = fit_model_b(batches)
    save(model_b.score_plot(pc_horiz=2, pc_vert=3, settings=LABELS), "model-b-scores-t2-t3")
    save(
        unfolded_contribution_plot(diagnose_different_batch(model_b, batches), 39, by_tag=True),
        "score-contributions-39-t3",
    )
    save(plot_raw(batches, "Press-3", [39]), "raw-Press-3-batch-39")

    model_c = fit_model_c(batches)
    save(model_c.score_plot(settings=LABELS), "model-c-scores")
    verify_left_out(model_c, batches)
    observability_table(model_c)
    print(f"figures written to {args.output_dir}")
    return 0

References#