Batch PLS fault diagnosis: the simulated SBR reactor#

Case study for issue #156.

Styrene-butadiene rubber (SBR) is made by emulsion polymerization in a batch reactor. Six trajectories are recorded during each batch (reactor, cooling water and jacket temperatures, latex density, conversion and the energy released), and five quality attributes of the latex are measured at the end (composition, particle size, branching, cross-linking and polydispersity). The 53 batches were simulated from a first-principles model, which makes this a rare kind of case study: the fault is known (Nomikos and MacGregor, 1994). Batch 37 received 30% more organic impurity in the butadiene feed than the normal batches from its very start, and batch 34 50% more from midway through.

The complete script is sbr_batch_pls.py in this directory:

uv run python docs/user_guide/case_studies/batch/sbr_batch_pls.py --output-dir case-study-output/sbr

It prints the numbers quoted below and writes every figure as an HTML file to the output directory.

Data#

SBR batch reactor: 53 batches of 200 samples with nine trajectories, and five quality attributes per batch. process_improve.batch.load_sbr() downloads the workbook and returns the batch dictionary, the quality table, and the list of the six trajectories the original study modelled. The two feed flows are constant in the simulation and the feed temperature barely moves, so they are left out.

CONF_LEVEL = 0.95
FAULT_FROM_START = 37
FAULT_MID_BATCH = 34
FAULT_BATCHES = [FAULT_MID_BATCH, FAULT_FROM_START]
INSPECT_SAMPLE = 120  # a sample well after the fault in batch 34 has developed
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) -> tuple[dict, pd.DataFrame]:
    """Download the batches and keep the six trajectories the original study modelled.

    The two feed flows are constant in the simulation and the feed
    temperature barely moves, so they carry no batch-to-batch information.
    """
    sbr = load_sbr(url=url)
    trajectories = {batch_id: batch[sbr.trajectory_tags] for batch_id, batch in sbr.X.items()}
    first = next(iter(trajectories.values()))
    print(
        f"{len(trajectories)} batches x {first.shape[0]} samples x {first.shape[1]} tags; quality block {sbr.Y.shape}"
    )
    return trajectories, sbr.Y



def plot_raw(trajectories: 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(
        trajectories, tag, batches_to_highlight={HIGHLIGHT: highlight}, extra_info=f"highlighting {highlight}"
    )

The batch PLS model#

process_improve.batch.BatchPLS unfolds each batch into one row of 6 tags times 200 samples, scales every column to unit variance, and fits a PLS model from that row to the five quality attributes.

def fit_model(trajectories: dict, quality: pd.DataFrame) -> BatchPLS:
    """Two-component batch PLS from the unfolded trajectories to the five quality attributes.

    Each batch is one row of 6 tags x 200 samples = 1200 columns. Every
    column is scaled to unit variance, so the mean of the per-column R2 is the
    R2 of the whole trajectory block.
    """
    model = BatchPLS(n_components=2).fit(trajectories, quality)
    r2y = per_component(model.r2_cumulative_.to_numpy())
    r2x = per_component(model.r2_per_variable_.mean(axis=0).to_numpy())
    print(f"R2X per component = {r2x[0]:.3f}, {r2x[1]:.3f}; R2Y per component = {r2y[0]:.3f}, {r2y[1]:.3f}")
    t1, t2 = model.scores_.iloc[:, 0], model.scores_.iloc[:, 1]
    print(f"lowest t1 = {t1.nsmallest(3).index.tolist()}; highest t2 = {t2.nlargest(2).index.tolist()}")
    t2_limit = model.hotellings_t2_limit(conf_level=CONF_LEVEL)
    spe_limit = model.spe_limit(conf_level=CONF_LEVEL)
    for batch_id in FAULT_BATCHES:
        hotelling = model.hotellings_t2_.loc[batch_id].iloc[-1]
        spe = model.spe_.loc[batch_id].iloc[-1]
        print(f"batch {batch_id}: T2 = {hotelling:.1f} (limit {t2_limit:.1f}), SPE = {spe:.1f} (limit {spe_limit:.1f})")
    return model


def per_component(cumulative: np.ndarray) -> np.ndarray:
    """Turn a cumulative R2 vector into the increment each component adds."""
    return np.diff(np.concatenate([[0.0], cumulative]))



R2X per component = 0.245, 0.127; R2Y per component = 0.653, 0.069
lowest t1 = [37, 34, 38]; highest t2 = [34, 9]
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 24.5% of the trajectories and 65.3% of the quality block; the second adds 12.7% and 6.9%. The score plot flags both faulty batches, which is encouraging: batch 37 has the lowest \(t_1\) of all batches, and batch 34 the highest \(t_2\). Both are far outside the Hotelling’s \(T^2\) limit. The SPE, on the other hand, flags neither. The SPE of a whole batch averages the residuals over 200 samples, so a deviation that the model can describe (a shift along the components) does not show up there. SPE and scores answer different questions.

Where the model explains the trajectories#

Every unfolded column has its own \(R^2\), so the fit can be read per tag and per time sample.

def r2_breakdown(model: BatchPLS) -> pd.DataFrame:
    """R2 of every (tag, time) cell after two components, as a tags x time grid."""
    grid = model.r2_per_variable_.iloc[:, -1].unstack(level="sequence")  # noqa: PD010 - inverse of the unfold
    print(
        "R2 per tag, averaged over time: " + ", ".join(f"{tag} {value:.2f}" for tag, value in grid.mean(axis=1).items())
    )
    return grid


def plot_r2_over_time(grid: pd.DataFrame) -> go.Figure:
    """One line per tag: how much of each trajectory the model explains at every sample."""
    fig = go.Figure()
    for tag, row in grid.iterrows():
        fig.add_trace(go.Scatter(x=list(grid.columns), y=row.to_numpy(), mode="lines", name=str(tag)))
    fig.update_layout(
        title="R2 of each trajectory over the batch", xaxis_title="Time [sequence order]", yaxis_title="R2"
    )
    return fig



R2 per tag, averaged over time: Conversion 0.75, CoolingTemp 0.23, EnergyReleased 0.26,
                                JacketTemp 0.24, LatexDensity 0.67, ReactorTemp 0.08

Latex density and conversion are the trajectories the model uses most, and \(R^2\) is low at the start of every trajectory because all batches begin alike. The time-varying weights \(w_1\) and \(w_2\), drawn with process_improve.batch.time_varying_loading_plot(), show the same picture per component.

Batch 37: the fault from the start#

def diagnose_batch_37(model: BatchPLS, trajectories: dict) -> pd.DataFrame:
    """Score contributions to t1: why batch 37 sits at the low end of t1."""
    contributions = model.score_contributions(model.unfold_and_scale(trajectories), component=1)
    row = contributions.loc[FAULT_FROM_START]
    by_tag = row.groupby(level="tag", sort=False).sum()
    print(f"batch {FAULT_FROM_START}: t1 contributions per tag = " + describe(by_tag))
    print(f"batch {FAULT_FROM_START}: share of the t1 contribution per fifth of the batch = {share_per_fifth(row)}")
    print(
        f"batch {FAULT_FROM_START}: first sustained departure from the other batches: {sustained_departure(trajectories, FAULT_FROM_START)}"
    )
    return contributions


def describe(by_tag: pd.Series) -> str:
    """Format per-tag contributions as one line."""
    return ", ".join(f"{tag} {value:+.1f}" for tag, value in by_tag.items())


def share_per_fifth(row: pd.Series) -> str:
    """Split a batch's contribution vector into five equal time blocks and give each block's share."""
    by_time = row.groupby(level="sequence").sum()
    fifths = by_time.groupby(np.arange(len(by_time)) * 5 // len(by_time)).sum()
    return ", ".join(f"{share:.0%}" for share in fifths / fifths.sum())


def sustained_departure(trajectories: dict, batch_id: int, n_sd: float = 2.0, run: int = 20) -> dict[str, int | None]:
    """Return the first sample from which each tag of one batch stays outside the band of the normal batches.

    The band is the mean of the batches without the fault plus or minus ``n_sd``
    standard deviations, sample by sample. A tag has departed once it stays
    outside the band for ``run`` consecutive samples, a tenth of the batch by
    default, so that a single noisy excursion does not count. ``None`` means
    the tag never departs.
    """
    others = np.stack([batch.to_numpy() for key, batch in trajectories.items() if key not in FAULT_BATCHES])
    z = (trajectories[batch_id].to_numpy() - others.mean(axis=0)) / others.std(axis=0, ddof=1)
    outside = (np.abs(z) > n_sd).astype(int)
    window = np.ones(run, dtype=int)
    onset: dict[str, int | None] = {}
    for j, tag in enumerate(trajectories[batch_id].columns):
        runs = np.convolve(outside[:, j], window, mode="valid") == run
        onset[tag] = int(np.argmax(runs)) if runs.any() else None
    return onset



batch 37: t1 contributions per tag = Conversion -33.4, CoolingTemp -4.2, EnergyReleased -5.2,
                                    JacketTemp -4.3, LatexDensity -25.5, ReactorTemp -1.3
batch 37: share of the t1 contribution per fifth of the batch = 15%, 18%, 19%, 26%, 22%
batch 37: first sustained departure from the other batches: {'ReactorTemp': None,
          'CoolingTemp': None, 'JacketTemp': None, 'LatexDensity': 13, 'Conversion': 9, 'EnergyReleased': None}

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: every fifth of it carries between 15% and 26% of the total. The raw data confirm it. sustained_departure expresses each trajectory of a faulty batch as a distance from the mean of the normal batches, in units of their standard deviation at that sample, and reports the first sample from which a tag stays more than two standard deviations away for twenty samples in a row (a single crossing is not informative, because a noisy tag crosses that line now and then in every batch). Conversion and latex density of batch 37 depart at samples 9 and 13 and run under the other batches to the end; none of its other four trajectories stays outside the band for twenty samples. The impurity slowed the reaction from the start, which is the injected fault.

Batch 34: the same fault, from the middle of the batch#

def diagnose_batch_34(model: BatchPLS, trajectories: dict) -> pd.DataFrame:
    """Score contributions to t2: the same fault, but starting midway through batch 34."""
    contributions = model.score_contributions(model.unfold_and_scale(trajectories), component=2)
    row = contributions.loc[FAULT_MID_BATCH]
    by_tag = row.groupby(level="tag", sort=False).sum()
    print(f"batch {FAULT_MID_BATCH}: t2 contributions per tag = " + describe(by_tag))
    print(f"batch {FAULT_MID_BATCH}: share of the t2 contribution per fifth of the batch = {share_per_fifth(row)}")
    print(
        f"batch {FAULT_MID_BATCH}: first sustained departure from the other batches: {sustained_departure(trajectories, FAULT_MID_BATCH)}"
    )
    return contributions



batch 34: t2 contributions per tag = Conversion +8.1, CoolingTemp +12.2, EnergyReleased +14.3,
                                    JacketTemp +12.3, LatexDensity +7.2, ReactorTemp +4.0
batch 34: share of the t2 contribution per fifth of the batch = 6%, 9%, 17%, 39%, 29%
batch 34: first sustained departure from the other batches: {'ReactorTemp': None,
          'CoolingTemp': 103, 'JacketTemp': 104, 'LatexDensity': 129, 'Conversion': 123, 'EnergyReleased': 105}

Batch 34 is high on \(t_2\), and the contributions come from the energy released, the jacket temperature and the cooling-water temperature. The timing is different from batch 37 in both views. The first two fifths of the batch carry only 15% of the \(t_2\) contribution and the last two fifths carry 68%, and in the raw data the cooling-water temperature, the jacket temperature and the energy released leave the band of the other batches at samples 103 to 105, the middle of the batch, while conversion and latex density only do so at samples 123 and 129. The same impurity, injected midway, shows up first in the heat balance of the reactor and only afterwards in the extent of reaction. process_improve.batch.contribution_at_time_plot() at sample 120 shows the same three tags carrying the deviation.

The same fault appears in two different 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. This is what makes batch models useful for diagnosis, and it is also why a library of “known faults” in score space needs the onset time as a coordinate.

Predicted quality#

def compare_predictions(model: BatchPLS, quality: pd.DataFrame) -> pd.DataFrame:
    """Observed and fitted quality of the two faulty batches, with the rank of each observed value."""
    table = pd.concat(
        {
            "observed": quality.loc[FAULT_BATCHES],
            "predicted": model.predictions_.loc[FAULT_BATCHES],
            "rank of observed": quality.rank().loc[FAULT_BATCHES].astype(int),
        },
        names=["value", "batch_id"],
    )
    print(f"quality of the faulty batches (rank 1 = lowest of {len(quality)} batches)")
    print(table.to_string(float_format=lambda value: f"{value:.4g}"))
    return table



quality of the faulty batches (rank 1 = lowest of 53 batches)
                           Composition  ParticleSize  Branching  CrossLinking  Polydispersity
value            batch_id
observed         34             0.4525          1244  1.234e-05     4.784e-05           3.599
                 37             0.4525          1247  1.173e-05     4.549e-05           3.462
predicted        34              0.454          1245  1.228e-05     4.761e-05           3.577
                 37               0.45          1250  1.183e-05     4.585e-05           3.491
rank of observed 34                  5             1          4             4              17
                 37                  4             2          1             1               1

Both batches produced poor latex: batch 37 has the lowest branching, cross-linking and polydispersity of all 53 batches and batch 34 the smallest particle size. The fitted values from the PLS model place them at the same end of every attribute, so a quality prediction from the trajectories would have flagged both batches before the laboratory did. Batch 37 is predicted low on every attribute; batch 34 is predicted only mildly low on polydispersity, because the \(t_2\) direction that carries its fault explains 6.9% of the quality block.

Predicting quality before the batch ends#

The model above was fitted on complete batches, and a plant would like to know the quality while the batch is still running. The unfolded row of a running batch is complete up to the newest sample and missing after it, so this is a missing-data problem: after \(k\) samples, estimate the scores from the cells observed so far, and read the quality prediction off the scores through the model’s Y loadings. The scores are estimated with trimmed score regression (Arteaga and Ferrer, 2002), a regression of the scores on the observed cells built from the training batches, the estimator that Garcia-Munoz, Kourti and MacGregor (2004) found gives stable score estimates from the first samples of a batch. process_improve.batch.BatchPLS.predict_online() does this for one point in time, and process_improve.batch.BatchPLS.predict_online_trace() for every sample of a complete batch, as the prediction would have evolved in real time.

AVERAGE_BATCH = 4  # the batch whose trajectories lie closest to the average
PREDICTION_SAMPLES = [10, 25, 50, 100, 150, 200]
ERROR_SAMPLES = [10, 50, 100, 150, 200]


def evolving_prediction(
    model: BatchPLS, trajectories: dict, quality: pd.DataFrame, batch_id: int = AVERAGE_BATCH
) -> pd.DataFrame:
    """Predict the final quality after every sample of one batch, as it would have been seen in real time.

    The unfolded row of a running batch is complete up to the newest sample
    and missing after it. ``predict_online_trace`` estimates the scores from
    the observed cells alone (trimmed score regression by default) and maps
    them to the quality attributes, once per sample.
    """
    y_hat = model.predict_online_trace(trajectories[batch_id]).y_hat
    actual = quality.loc[batch_id, "ParticleSize"]
    fitted = model.predictions_.loc[batch_id, "ParticleSize"]
    print(f"batch {batch_id}: ParticleSize measured {actual:.1f}, fitted from the complete batch {fitted:.1f}")
    print(
        f"batch {batch_id}: ParticleSize predicted after "
        + ", ".join(f"{k} samples {y_hat.loc[k, 'ParticleSize']:.1f}" for k in PREDICTION_SAMPLES)
    )
    return y_hat


def evolving_error(model: BatchPLS, trajectories: dict, quality: pd.DataFrame) -> pd.DataFrame:
    """RMSEE of the evolving prediction at every sample, relative to the standard deviation of each attribute.

    ``online_rmse`` traces every training batch and pools the squared errors
    per sample, so this is the estimation error (RMSEE) as a function of how
    much of the batch has been observed. A ratio of about 1 is the error of
    predicting the average batch.
    """
    ratio = model.online_rmse(trajectories, quality) / quality.std(ddof=1)
    for target in ("ParticleSize", "Branching"):
        print(
            f"RMSEE / sd of {target} after "
            + ", ".join(f"{k} samples {ratio.loc[k, target]:.2f}" for k in ERROR_SAMPLES)
        )
    return ratio


def plot_online_rmse(ratio: pd.DataFrame) -> go.Figure:
    """One line per quality attribute: the RMSEE relative to that attribute's standard deviation."""
    fig = go.Figure()
    for target in ratio.columns:
        fig.add_trace(go.Scatter(x=ratio.index, y=ratio[target].to_numpy(), mode="lines", name=str(target)))
    fig.add_hline(y=1.0, line_dash="dash", annotation_text="predicting the average batch")
    fig.update_layout(
        title="RMSEE of the evolving quality prediction (training batches)",
        xaxis_title="Samples observed",
        yaxis_title="RMSEE / standard deviation",
    )
    return fig


def plot_online_prediction(y_hat: pd.DataFrame, quality: pd.DataFrame, batch_id: int) -> go.Figure:
    """Draw the evolving prediction of one batch's particle size against its measured value and the average batch."""
    fig = go.Figure()
    fig.add_trace(go.Scatter(x=y_hat.index, y=y_hat["ParticleSize"].to_numpy(), mode="lines", name="predicted"))
    fig.add_hline(y=quality.loc[batch_id, "ParticleSize"], line_dash="dash", annotation_text="measured")
    fig.add_hline(y=quality["ParticleSize"].mean(), line_dash="dot", annotation_text="average of all batches")
    fig.update_layout(
        title=f"ParticleSize of batch {batch_id}, predicted as the batch runs",
        xaxis_title="Samples observed",
        yaxis_title="ParticleSize",
    )
    return fig



batch 4: ParticleSize measured 1256.9, fitted from the complete batch 1257.1
batch 4: ParticleSize predicted after 10 samples 1251.0, 25 samples 1248.1, 50 samples 1255.3,
         100 samples 1254.9, 150 samples 1257.4, 200 samples 1257.1
RMSEE / sd of ParticleSize after 10 samples 2.84, 50 samples 1.40, 100 samples 1.29,
                               150 samples 0.62, 200 samples 0.60
RMSEE / sd of Branching after 10 samples 2.10, 50 samples 0.88, 100 samples 0.81,
                            150 samples 0.32, 200 samples 0.24

Batch 4 is the batch whose trajectories lie closest to the average. Its particle size is predicted 6 to 9 units away from the measured value in the first 25 samples and within about 2 units from sample 50 onwards, and after 200 samples the prediction equals the fitted value from the complete batch, as it must, because the row is then complete. One batch says little about the error, though. process_improve.batch.BatchPLS.online_rmse() traces every training batch and pools the squared errors sample by sample; on the training batches this is the root-mean-square error of estimation, RMSEE, as a function of how much of the batch has been observed, and its last value is the RMSEE of the model fitted on complete batches. Dividing by the standard deviation of each attribute puts the five attributes on one axis, where a ratio of about 1 is the error of predicting the average batch every time.

For particle size the ratio is 2.84 after 10 samples, still 1.29 at the halfway point and 0.62 after 150 samples: in the first half of the batch the prediction is no better than the average batch, and it improves in the second half. Branching, whose final RMSEE is a quarter of its standard deviation, is below 1 by sample 50. The early predictions are worse than the average because few cells have been observed, and they are the cells where every batch begins alike (the \(R^2\) breakdown showed this), so the score estimate carries little information about the batch. The RMSEE is measured on the batches the model was fitted to. Refitting the model with one batch left out and tracing that batch, which the script leaves out because it refits 53 models, gives a prediction error, RMSEP, of 3.10, 1.66, 1.65, 0.84 and 0.78 standard deviations at the same five samples: the same shape, at a higher level.

Would the model have caught it on-line?#

The score plot flagged both faulty batches once they were complete. The question here is whether a chart could have flagged them while they were running, and after how many samples. Two things change from the model above. A reference model must describe normal operation, so batches 34 and 37 are left out and a two-component model is refitted on the other 51 batches. And a running batch needs a limit at every sample rather than one limit for the whole batch: the score estimates early in a batch are shrunk and noisy compared with those near its end, so the scatter of the reference batches differs from sample to sample. process_improve.batch.BatchMonitor passes every reference batch through predict_online_trace and summarises the spread at each sample, following Nomikos and MacGregor (1995), who set the limits of their score charts from that spread and note that the \(T^2\) chart needs the score covariance at each sample as well. The \(T^2\) at sample \(k\) is standardised by the covariance of the reference batches’ score estimates at that sample, as Garcia-Munoz, Kourti and MacGregor (2004) compute it, which gives one limit for the whole batch, and the SPE limit is a chi-squared limit fitted to the reference batches’ SPE at that sample (spe_window pools the values of neighbouring samples into that fit, which steadies the limits when few reference batches are available; on these 51 batches a window of two samples either side leaves every alarm sample unchanged). The SPE charted is the instantaneous one, the residual of the newest sample only, which reacts in the sample a fault begins; the cumulative SPE over every cell observed so far is diluted by the earlier, normal samples, and here it reacts to batch 34 after 112 samples instead of 105. process_improve.batch.online_monitoring_plot() draws either chart.

MONITOR_CONF_LEVEL = 0.99
ALARM_RUN = 3  # consecutive samples above the limit before an alarm counts
FORECAST_FROM = {FAULT_FROM_START: ("Conversion", [30, 60]), FAULT_MID_BATCH: ("CoolingTemp", [60, 115])}


def normal_batches(trajectories: dict) -> dict:
    """Return the batches without the injected fault."""
    return {batch_id: batch for batch_id, batch in trajectories.items() if batch_id not in FAULT_BATCHES}


def online_monitoring(trajectories: dict, quality: pd.DataFrame) -> tuple[BatchMonitor, dict]:
    """Per-sample T2 and SPE charts of the two faulty batches against a model of the normal batches.

    A reference model must describe normal operation, so batches 34 and 37
    are left out of it. ``BatchMonitor`` projects every reference batch after
    1, 2, ..., 200 samples with the same missing-data estimator and builds a
    limit at each sample from the spread of those projections. The SPE
    charted is the instantaneous one, the residual of the newest sample only,
    which reacts in the sample a fault begins. An alarm counts once the
    statistic stays above its limit for three consecutive samples.
    """
    normal = normal_batches(trajectories)
    reference = BatchPLS(n_components=2).fit(normal, quality.loc[list(normal)])
    monitor = BatchMonitor(reference, conf_level=MONITOR_CONF_LEVEL, spe_statistic="instantaneous").fit(normal)
    print(f"reference model on {len(normal)} batches; T2 limit {monitor.t2_limit_over_time_[0]:.2f} at every sample")
    alarms: dict[int, dict[str, int | None]] = {}
    for batch_id in FAULT_BATCHES:
        result = monitor.monitor(trajectories[batch_id])
        alarms[batch_id] = {
            "T2": first_sustained_alarm(result.t2_alarm),
            "SPE": first_sustained_alarm(result.spe_alarm),
        }
        print(
            f"batch {batch_id}: first {ALARM_RUN} consecutive samples above the limit: "
            f"T2 after {alarms[batch_id]['T2']} samples, SPE after {alarms[batch_id]['SPE']} samples"
        )
    t2_rate, spe_rate = alarm_rates(monitor, normal)
    print(f"normal batches: {t2_rate:.2%} of the T2 values and {spe_rate:.2%} of the SPE values above their limits")
    spe_alarm = alarms[FAULT_MID_BATCH]["SPE"]
    if spe_alarm is None:
        raise RuntimeError(
            f"batch {FAULT_MID_BATCH} raised no sustained SPE alarm; the residual shares below need one."
        )
    for k in (spe_alarm, spe_alarm + 4):
        shares = residual_shares(reference, trajectories[FAULT_MID_BATCH], k)
        print(f"batch {FAULT_MID_BATCH}: share of the SPE after {k} samples per tag = " + describe_shares(shares))
    for batch_id, (tag, sample_points) in FORECAST_FROM.items():
        for k in sample_points:
            forecast = reference.predict_online(trajectories[batch_id], k).forecast
            print(
                f"batch {batch_id} {tag}, mean of the samples after {k}: "
                f"forecast {forecast[tag].iloc[k:].mean():.4f}, actual {trajectories[batch_id][tag].iloc[k:].mean():.4f}, "
                f"average of the normal batches {np.mean([batch[tag].iloc[k:].mean() for batch in normal.values()]):.4f}"
            )
    return monitor, alarms


def first_sustained_alarm(alarm: np.ndarray, run: int = ALARM_RUN) -> int | None:
    """Return the 1-based sample after which a statistic first stays above its limit for ``run`` samples in a row."""
    runs = np.convolve(np.asarray(alarm, dtype=int), np.ones(run, dtype=int), mode="valid") == run
    return int(np.argmax(runs)) + 1 if runs.any() else None


def alarm_rates(monitor: BatchMonitor, batches: dict) -> tuple[float, float]:
    """Fraction of the (batch, sample) points of these batches above the T2 limit and above the SPE limit."""
    results = [monitor.monitor(batch) for batch in batches.values()]
    t2_rate = float(np.mean([result.t2_alarm for result in results]))
    spe_rate = float(np.mean([result.spe_alarm for result in results]))
    return t2_rate, spe_rate


def residual_shares(model: BatchPLS, batch: pd.DataFrame, k: int) -> pd.Series:
    """Share of the squared instantaneous SPE after ``k`` samples carried by each tag."""
    squared = model.predict_online(batch, k).residuals.xs(k - 1, level="sequence") ** 2
    return squared / squared.sum()


def describe_shares(shares: pd.Series) -> str:
    """Format per-tag shares as one line, largest first."""
    return ", ".join(f"{tag} {share:.0%}" for tag, share in shares.sort_values(ascending=False).items())


def plot_forecasts(reference: BatchPLS, trajectories: dict) -> go.Figure:
    """Draw the rest of each faulty batch as forecast from its scores at two points, next to what happened."""
    normal = normal_batches(trajectories)
    titles = [f"batch {batch_id}: {tag}" for batch_id, (tag, _) in FORECAST_FROM.items()]
    fig = make_subplots(rows=1, cols=2, subplot_titles=titles)
    for col, (batch_id, (tag, sample_points)) in enumerate(FORECAST_FROM.items(), start=1):
        n = len(trajectories[batch_id])
        time = np.arange(1, n + 1)
        average = np.mean([batch[tag].to_numpy() for batch in normal.values()], axis=0)
        fig.add_trace(go.Scatter(x=time, y=average, mode="lines", name="average of the normal batches"), row=1, col=col)
        fig.add_trace(
            go.Scatter(x=time, y=trajectories[batch_id][tag].to_numpy(), mode="lines", name=f"batch {batch_id}"),
            row=1,
            col=col,
        )
        for k in sample_points:
            forecast = reference.predict_online(trajectories[batch_id], k).forecast[tag]
            fig.add_trace(
                go.Scatter(
                    x=time[k:],
                    y=forecast.iloc[k:].to_numpy(),
                    mode="lines",
                    line={"dash": "dash"},
                    name=f"forecast from {k} samples",
                ),
                row=1,
                col=col,
            )
        fig.update_xaxes(title_text="Time [sequence order]", row=1, col=col)
        fig.update_yaxes(title_text=tag, row=1, col=col)
    fig.update_layout(title="The rest of the batch, forecast from the reference model")
    return fig



reference model on 51 batches; T2 limit 10.54 at every sample
batch 34: first 3 consecutive samples above the limit: T2 after 190 samples,
          SPE after 105 samples
batch 37: first 3 consecutive samples above the limit: T2 after 23 samples,
          SPE after 145 samples
normal batches: 0.17% of the T2 values and 1.24% of the SPE values above their limits
batch 34: share of the SPE after 105 samples per tag = ReactorTemp 40%, CoolingTemp 30%,
          JacketTemp 16%, EnergyReleased 11%, Conversion 2%, LatexDensity 1%
batch 34: share of the SPE after 109 samples per tag = CoolingTemp 36%, JacketTemp 29%,
          EnergyReleased 19%, ReactorTemp 12%, Conversion 2%, LatexDensity 1%
batch 37 Conversion, mean of the samples after 30: forecast 0.6472, actual 0.6420,
                     average of the normal batches 0.6609
batch 37 Conversion, mean of the samples after 60: forecast 0.6480, actual 0.6477,
                     average of the normal batches 0.6656
batch 34 CoolingTemp, mean of the samples after 60: forecast 46.6357, actual 46.7744,
                      average of the normal batches 46.6152
batch 34 CoolingTemp, mean of the samples after 115: forecast 46.6366, actual 46.8173,
                      average of the normal batches 46.6015

An alarm counts once the statistic stays above its 99% limit for three consecutive samples. A single crossing is not informative: the normal batches put 0.17% of their \(T^2\) values and 1.24% of their SPE values above the limits, so a batch of 200 samples crosses the SPE limit now and then (batch 34 has one isolated crossing after 37 samples, and batch 4, the near-average batch, one after 26).

Batch 37 shows in the \(T^2\) chart after 23 samples: its \(T^2\) is 9.6 after 21 samples, 10.2 after 22, 10.7 after 23 and 11.9 after 25, against a limit of 10.54, and it stays above the limit for the rest of the batch. The impurity slowed the reaction from the first sample, and low conversion and latex density is the direction the reference model’s first component describes, so the score estimate moves along \(t_1\) as soon as enough samples have been observed to estimate it. The instantaneous SPE of batch 37 does not react until sample 145; the deviation lies in the model plane and leaves little residual.

Batch 34 shows in the SPE chart after 105 samples, within a couple of samples of the point where the cooling-water temperature, the jacket temperature and the energy released leave the band of the other batches in the raw data (samples 103 to 105). At the alarm sample the reactor temperature carries 40% of the squared residual and the cooling-water temperature 30%; four samples later the cooling-water temperature, the jacket temperature and the energy released carry 84% between them, the same three tags the whole-batch contribution plot named. The \(T^2\) of batch 34 does not alarm until sample 190. None of the reference batches has a heat-balance deviation that starts midway, so the reference model has no component for it; the deviation is off the model plane and shows in the residual, not in the scores. The two charts answer different questions: the \(T^2\) reacts to a deviation along the components the reference batches exhibited, the SPE to a direction the reference model never saw, and which chart catches a fault first depends on the fault.

The same distinction decides whether the model can forecast the rest of a running batch. predict_online also returns a forecast: the batch’s own values up to the newest sample, and beyond it the trajectories implied by the score estimate, Eq. 4 of Wold, Kettaneh-Wold, MacGregor and Dunn (2009). For batch 37 the forecast from 30 samples puts the mean conversion of the remaining samples at 0.6472, against 0.6420 observed and 0.6609 for the average normal batch; from 60 samples the forecast is 0.6480 against 0.6477 observed. The scores had picked up the slow reaction, and the forecast follows it. For batch 34 the forecast of the cooling-water temperature stays at the average of the normal batches, 46.64 from 60 samples and again from 115 samples, while the batch ran at 46.77 and 46.82. The fault of batch 34 does not move the scores of the reference model, and a forecast made from the scores cannot see it.

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/sbr"))
    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")

    trajectories, quality = load_data(args.data_url)
    for tag in ("Conversion", "LatexDensity", "CoolingTemp"):
        save(plot_raw(trajectories, tag, FAULT_BATCHES), f"raw-{tag}")

    model = fit_model(trajectories, quality)
    save(model.score_plot(settings=LABELS), "scores")
    save(model.spe_plot(settings=LABELS), "spe")
    save(plot_r2_over_time(r2_breakdown(model)), "r2-over-time")
    save(time_varying_loading_plot(model, component=1), "weights-w1")
    save(time_varying_loading_plot(model, component=2), "weights-w2")

    t1_contributions = diagnose_batch_37(model, trajectories)
    save(unfolded_contribution_plot(t1_contributions, FAULT_FROM_START), "contributions-37-t1")
    save(unfolded_contribution_plot(t1_contributions, FAULT_FROM_START, by_tag=True), "contributions-37-t1-by-tag")
    for tag in ("LatexDensity", "Conversion"):
        save(plot_raw(trajectories, tag, [FAULT_FROM_START]), f"raw-{tag}-batch-37")

    t2_contributions = diagnose_batch_34(model, trajectories)
    save(unfolded_contribution_plot(t2_contributions, FAULT_MID_BATCH), "contributions-34-t2")
    save(
        contribution_at_time_plot(t2_contributions, k=INSPECT_SAMPLE, batch_id=FAULT_MID_BATCH),
        "contributions-34-t2-at-sample",
    )
    for tag in ("CoolingTemp", "JacketTemp", "EnergyReleased"):
        save(plot_raw(trajectories, tag, [FAULT_MID_BATCH]), f"raw-{tag}-batch-34")

    compare_predictions(model, quality)
    for variable in ("Composition", "ParticleSize"):
        save(model.predictions_vs_observed_plot(quality, variable=variable), f"observed-vs-predicted-{variable}")

    y_hat = evolving_prediction(model, trajectories, quality)
    save(plot_online_prediction(y_hat, quality, AVERAGE_BATCH), f"online-prediction-batch-{AVERAGE_BATCH}")
    save(plot_online_rmse(evolving_error(model, trajectories, quality)), "online-rmse")

    monitor, _alarms = online_monitoring(trajectories, quality)
    save(online_monitoring_plot(monitor, trajectories[FAULT_FROM_START], "t2"), f"online-t2-batch-{FAULT_FROM_START}")
    save(online_monitoring_plot(monitor, trajectories[FAULT_MID_BATCH], "spe"), f"online-spe-batch-{FAULT_MID_BATCH}")
    save(plot_forecasts(monitor.model, trajectories), f"forecast-batch-{FAULT_FROM_START}-and-{FAULT_MID_BATCH}")
    print(f"figures written to {args.output_dir}")
    return 0

References#