17 min reading time

Soft sensors and inferential sensors

7.4. Soft sensors and inferential sensors

A soft sensor (also called an inferential sensor) infers a hard-to-measure quality variable from cheap, real-time process measurements that are already on the data historian. The use of PLS for exactly this purpose was introduced by Kresta, Marlin and MacGregor (1994). The next subsection outlines what these sensors are and how they are built; the rest of this section is a worked example: predicting the Kappa number on a continuous Kamyr pulp digester, where the Kappa number is reported less often, and with more delay, than the process tags around it.

7.4.1. What a soft sensor is

The intention of an inferential sensor is to infer a hard-to-measure property, usually a lab measurement or an expensive measurement, using a combination of process data and software-implemented algorithms. These sensors also go by the names of software sensors or just soft sensors.

Consider a distillation column where various automatic measurements are used to predict the vapour pressure. The actual vapour pressure is a lab measurement, usually taken 3 or 4 times per week, and takes several hours to complete. The soft sensor can predict the lab value from the real-time process measurements with sufficient accuracy. This is a common soft sensor on distillation columns. The lab values are used to build (train) the software sensor and to update in periodically.

Other interesting examples use camera images to predict hard-to-measure values. In the paper by Honglu Yu, John MacGregor, Gabe Haarsma and Wilfred Bourg (Ind. Eng. Chem. Res., 42, 3036–3044, 2003), the authors describe how machine vision is used to predict, in real-time, the seasoning of various snack-food products. This sensors uses the colour information of the snacks to infer the amount of seasoning dispensed onto them. The dispenser is controlled via a feedback loop to ensure the seasoning is at target.

Once validated, a soft sensor can also reduce costs of a process by allowing for rapid feedback control of the inferred property, so that less off-specification product is produced. They also often have the side-effect that reduced lab sampling is required; this saves on manpower costs.

Soft sensors using latent variables will almost always be PLS models. Once the model has been built, it can be applied in real-time. The \(T^2\) and SPE value for each new observation is checked for consistency with the model before a prediction is made. Contribution plots are used to diagnose unusual observations.

It is an indication that the predictive models need to be updated if the SPE and/or \(T^2\) values are consistently above the limits. This is a real advantage over using an MLR-based model, which has no such consistency checks.

7.4.2. Recap: process monitoring and the lab-measurement gap

Recall from the earlier chapter on process monitoring that control charts are built in two phases: phase 1 fits a chart to a stretch of stable historical data, and phase 2 uses that chart on new, real-time data to flag unusual variability so that an operator can intervene. We were careful to point out that monitoring is not feedback control: the adjustments that follow an alarm are infrequent, manual, and made only when a special cause has been identified.

This procedure has a practical requirement: the variable being charted must be available in real-time. A temperature, a flow rate, a pressure or an on-line composition analyser update every few seconds and present no difficulty. The case that does not fit is the lab measurement: a final quality property that is sampled at intervals of an hour or more, and – depending on whether the mill has an on-line analyser or a wet-chemistry titration in a busy lab – reported with a turnaround of anywhere from a minute to a couple of hours. Even when the chemistry itself is fast, the sampling interval is enough to leave a chart that updates only a few times per shift, which is too sparse to react to disturbances before the off-spec material has moved downstream.

A soft sensor solves this problem by inferring the lab value in real-time from the process tags that are available at that moment. The model is built on historical data where both the process tags and the lab values were collected; once it has been validated the prediction is used in place of the lab value on the monitoring chart, on the on-line trend, or in a feedback loop. The phase 1 and phase 2 requirements still apply: the model is fit on a representative stretch of data, then tested on data it has not seen before it is deployed. For a step-by-step procedure that covers the engineering choices along the way, see Lin, Recke, Knudsen and Jørgensen (2007).

7.4.3. Case study: predicting Kappa number on a Kamyr digester

A continuous Kamyr digester cooks wood chips under pressure in white-liquor to dissolve the lignin that binds the cellulose fibres. The amount of lignin that remains in the pulp is summarised by a single quality number, the Kappa number. A high Kappa number indicates a brown, paperboard-grade pulp; a low Kappa number indicates a pulp that is closer to a bleachable grade. The mill aims to hold the Kappa number on target with as little variability as possible.

The Kappa number is conventionally a wet-chemistry lab measurement defined by the ISO 302 titration method, which is itself about 30 minutes of bench work; once sampling, transport, sample preparation, lab queueing and reporting are added, the practical turnaround on a busy mill is typically one to two hours, and the sampling frequency is usually no faster than that. Modern installations also offer on-line NIR analysers that can return a Kappa-equivalent value in under a minute, calibrated against ISO 302, but these are still far from universal; on mills that do not have one the fastest information about the digester is the existing process historian, which is what a soft sensor lets us turn into a real-time Kappa estimate. Dayal, MacGregor, Taylor, Kildaw and Marcikic (1994) demonstrated PLS and neural networks side by side on exactly this problem and on the same data set we use here.

The data set used here is the Kamyr digester data from openmv.net, an hourly record from a kraft mill in Alberta. The file is 301 rows by 23 columns: the first column is an Observation timestamp, Y-Kappa is the lab measurement that we want to predict, and the remaining 21 columns are process tags. Two of the process columns, AAWhiteSt-4 and SulphidityL-4, are missing on about half of the rows and we drop them, leaving 19 process tags in \(\mathbf{X}\).

Several of the tags have already been pre-aligned in time: a column name that ends in a digit indicates how many hours that tag has been shifted to line up with the Kappa number it eventually produces. For example, ChipLevel4 is the chip level lagged 4 hours, BlackFlow-2 is the black-liquor flow lagged 2 hours, SteamHeatF-3 is the steam-heater flow lagged 3 hours. The unlagged tags (ChipRate, BF-CMratio, BlowFlow, UCZAA, WeakLiquorF, WeakWashF) measure quantities that affect the pulp downstream of the digester or that change slowly enough that the lag is negligible. The residence-time estimates were provided by the mill along with the data.

We load the file, drop the timestamp and the two mostly-missing columns, and median-impute the remaining gaps. The same digester dataframe is reused throughout this section, so the plotting and modelling blocks below can be pasted in order to reproduce every figure:

Python, 25 lines
import numpy as np
import pandas as pd
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from process_improve.multivariate import PLS, MCUVScaler

digester = pd.read_csv("https://openmv.net/file/kamyr-digester.csv")
digester.columns = [c.strip() for c in digester.columns]
digester = digester.drop(columns=["Observation", "AAWhiteSt-4", "SulphidityL-4"])
digester = digester.fillna(digester.median(numeric_only=True))

sample = np.arange(len(digester))
fig = make_subplots(rows=3, cols=1, shared_xaxes=True, vertical_spacing=0.05)
fig.add_trace(go.Scatter(x=sample, y=digester["Y-Kappa"], mode="lines",
        line=dict(color="black"), showlegend=False), row=1, col=1)
fig.add_trace(go.Scatter(x=sample, y=digester["ChipLevel4"], mode="lines",
        line=dict(color="#1f77b4"), showlegend=False), row=2, col=1)
fig.add_trace(go.Scatter(x=sample, y=digester["BlackFlow-2"], mode="lines",
        line=dict(color="#d62728"), showlegend=False), row=3, col=1)
fig.update_yaxes(title_text="Y-Kappa", row=1, col=1)
fig.update_yaxes(title_text="ChipLevel4", row=2, col=1)
fig.update_yaxes(title_text="BlackFlow-2", row=3, col=1)
fig.update_xaxes(title_text="Sample (1 hour spacing)", row=3, col=1)
fig.update_layout(height=620, margin=dict(l=70, r=20, t=20, b=50))
fig.show()
Raw Kappa number and two of the lagged process tags plotted against sample number.

The Kappa number in the top panel covers a range of about 15 units (from 12 to 28 Kappa). The two pre-shifted process tags below it move on the same time-scale and provide the real-time information that the soft sensor will use.

7.4.4. Building the soft sensor

We build the model with PLS, using the Kappa number as the \(y\)-variable and the nineteen process tags as the \(\mathbf{X}\) block. The process-improve package provides a PLS class and an MCUVScaler for centring and scaling. The data are centred to zero mean and scaled to unit standard deviation before fitting: PLS scores and loadings are only interpretable in that form. The PLS class scales both blocks itself by default (scale=True); since the blocks are already standardised by MCUVScaler, the model is fitted with scale=False so the data are not scaled a second time.

Python, 8 lines
X = digester.drop(columns=["Y-Kappa"])
y = digester[["Y-Kappa"]]

scaler_x = MCUVScaler().fit(X)
scaler_y = MCUVScaler().fit(y)

model = PLS(n_components=2, scale=False).fit(scaler_x.transform(X),
                                             scaler_y.transform(y))

The cumulative \(R^2_Y\) measures how much of the Kappa variability the model accounts for as each latent variable is added:

Python, 2 lines
print(model.r2_cumulative_.round(3).to_list())
# [0.35, 0.503]

The first component picks up 35% of the Kappa variability, and the second adds another 15%, for a cumulative 50%. This is not a high-\(R^2_Y\) model in absolute terms, but for a soft sensor on a continuous process it is enough to be useful, as we will see when we evaluate it on held-out data below.

The regression coefficients show which tags drive the model:

Python, 7 lines
coefs = model.beta_coefficients_.iloc[:, 0]
colors = ["#1f77b4" if c >= 0 else "#d62728" for c in coefs.values]
fig = go.Figure(go.Bar(x=coefs.index, y=coefs.values, marker_color=colors))
fig.add_hline(y=0, line_color="black", line_width=0.6)
fig.update_layout(yaxis_title="Coefficient on scaled X", xaxis_tickangle=-40,
        height=420, margin=dict(l=70, r=20, t=20, b=120))
fig.show()
Bar chart of PLS regression coefficients onto Y-Kappa for the nineteen process tags.

PLS regression coefficients on the centred and scaled \(\mathbf{X}\). The bar heights are directly comparable. BF-CMratio, ChipLevel4, SteamHeatF-3, ChipRate, SteamFlow-4 and WhiteFlow-4 carry most of the relationship to Y-Kappa.

The signs match what is known about the chemistry of the cook. A higher SteamHeatF-3, SteamFlow-4 or WhiteFlow-4 puts more heat or active alkali into the digester, increases the rate of delignification, and pulls the Kappa number down – their coefficients are negative. A higher ChipRate or ChipLevel4 moves more wood through the digester for the same heat and alkali, leaves more residual lignin, and pushes the Kappa number up. The coefficients are a correlation, not a causation, but they agree with what a process engineer would predict from first principles.

The cumulative \(R^2_Y\) reports the fit on the data we trained on. To know whether the model will be useful as a live soft sensor we have to evaluate it on data it has never seen. We split the 301 rows in time order: the first 70% (211 rows) for training, the last 30% (90 rows) as a held-out test set, and we report the root-mean-square error of prediction (RMSEP) on the test set:

Python, 18 lines
def evaluate_split(df, x_cols, y_col, frac=0.70):
        n_train = int(round(frac * len(df)))
        train, test = df.iloc[:n_train], df.iloc[n_train:]
        sx = MCUVScaler().fit(train[x_cols])
        sy = MCUVScaler().fit(train[[y_col]])
        m = PLS(n_components=2, scale=False).fit(sx.transform(train[x_cols]),
                sy.transform(train[[y_col]]))
        y_hat_scaled = pd.DataFrame(
                np.asarray(m.predict(sx.transform(test[x_cols]))),
                index=test.index, columns=[y_col])
        y_hat = sy.inverse_transform(y_hat_scaled).values.ravel()
        y_obs = test[y_col].values
        rmsep = float(np.sqrt(np.mean((y_obs - y_hat) ** 2)))
        return rmsep, y_obs, y_hat

x_cols = list(X.columns)
rmsep_base, y_obs_base, y_hat_base = evaluate_split(digester, x_cols, "Y-Kappa")
print(f"RMSEP (process tags only): {rmsep_base:.2f} Kappa units")

The same predicted-vs-observed scatter is reused below for the lag-augmented model, so we define it once as a helper:

Python, 15 lines
def plot_obs_pred(y_obs, y_hat, title):
        lo = float(min(y_obs.min(), y_hat.min()))
        hi = float(max(y_obs.max(), y_hat.max()))
        pad = 0.05 * (hi - lo)
        fig = go.Figure()
        fig.add_trace(go.Scatter(x=[lo - pad, hi + pad], y=[lo - pad, hi + pad],
                mode="lines", line=dict(color="black", dash="dash", width=1), name="ideal"))
        fig.add_trace(go.Scatter(x=y_obs, y=y_hat, mode="markers",
                marker=dict(size=8, opacity=0.8), name="predictions"))
        fig.update_layout(title=title, xaxis_title="Observed Kappa",
                yaxis_title="Predicted Kappa", height=520, width=580)
        fig.update_yaxes(scaleanchor="x", scaleratio=1)
        fig.show()

plot_obs_pred(y_obs_base, y_hat_base, "Soft sensor predictions: process tags only")
Predicted vs observed Kappa number on the held-out test set with the process-tag model.

Predicted vs observed Kappa number on the held-out test set, using only the 19 process tags. The RMSEP is 1.96 Kappa units.

The held-out RMSEP is 1.96 Kappa units, on a Kappa number that varies between 12 and 28 in this dataset (standard deviation about 3 units). A soft sensor that is within 2 Kappa of the lab value is already useful, and a common refinement makes it a fair bit better: include the previous Kappa value as another predictor. Each time the lab returns a value we keep it as a one-step memory and feed it back into \(\mathbf{X}\) until the next lab value arrives. We add a column Kappa_lag1 to the data, drop the first row (which has no previous Kappa), and re-fit:

Python, 10 lines
df_lag = digester.copy()
df_lag["Kappa_lag1"] = df_lag["Y-Kappa"].shift(1)
df_lag = df_lag.dropna(subset=["Kappa_lag1"]).reset_index(drop=True)

rmsep_lag, y_obs_lag, y_hat_lag = evaluate_split(df_lag, x_cols + ["Kappa_lag1"],
                                                 "Y-Kappa")
print(f"RMSEP (with one-step Kappa lag): {rmsep_lag:.2f} Kappa units")

plot_obs_pred(y_obs_lag, y_hat_lag,
              "Soft sensor predictions: process tags + 1-step Kappa lag")
Predicted vs observed Kappa with one-step Kappa lag in X.

Adding the previous Kappa value as a twentieth predictor reduces the RMSEP from 1.96 to 1.73 Kappa units, an improvement of about 12% on the same underlying data.

The scatter plots tell us how close the predictions are on average, but they hide where the soft sensor disagrees with the lab. Plotting the two test-set predictions and the lab value against time shows both stories on one figure:

Python, 13 lines
sample = np.arange(len(y_obs_base))
fig = go.Figure()
fig.add_trace(go.Scatter(x=sample, y=y_obs_base, mode="lines",
        line=dict(color="black", width=2.4), name="Lab (actual)"))
fig.add_trace(go.Scatter(x=sample, y=y_hat_base, mode="lines+markers",
        line=dict(color="#1f77b4", dash="dash"), marker=dict(symbol="circle", size=6),
        name="Soft sensor: process tags only"))
fig.add_trace(go.Scatter(x=sample[-len(y_hat_lag):], y=y_hat_lag, mode="lines+markers",
        line=dict(color="#d62728", dash="dot"), marker=dict(symbol="square", size=6),
        name="Soft sensor: process tags + Kappa lag"))
fig.update_layout(xaxis_title="Test sample index (1 hour spacing)",
        yaxis_title="Kappa number", height=440)
fig.show()
Time-series overlay of held-out Kappa predictions and the actual lab values.

Held-out test predictions over time. The actual lab Kappa is the solid black line; the process-tags-only prediction is the blue dashed line with round markers; the lag-augmented prediction is the red dotted line with square markers. Both soft sensors track the slow upward drift in the second half, but the lag-augmented model stays closer to the lab values on sample-to-sample swings, which is what drives the RMSEP down from 1.96 to 1.73.

Note

The PLS coefficients are a measure of correlation, not causation. They are useful as a starting point for a discussion with the operators about what is driving the variability in the process, but they should not be interpreted as a cause-and-effect statement on their own.

This soft sensor can now be deployed in the same way as any other monitoring artefact: the model is fit on a representative stretch of historical operation, tested on a held-out tail, and once the RMSEP is small enough compared to the Kappa specification limits, the prediction is wired into the same chart that would have been used for the lab value. The chart now updates every hour rather than once per shift. Tzovla and Mehta (2002) describe how this build-test-deploy loop can be automated for soft sensors inside commercial control systems.

Two refinements are worth noting. First, the relationship between \(\mathbf{X}\) and the Kappa number changes over time, with chip species, mill upsets and seasonal raw-material variation. A soft sensor that is never re-fit will gradually lose accuracy. The standard practice is to re-fit the model on a rolling window of the most recent lab values; Kadlec, Grbić and Gabrys (2011) review the variations on this adaptation theme. Second, the Kappa number is autocorrelated in time: we exploited this with a single one-step lag of \(y\) above, but adding lags of two or three hours, and lags on the most influential \(x\) variables, typically reduces the prediction error further.

7.4.5. Further reading

Contents