Problems: Dimensionality Reduction

Problems: Dimensionality Reduction#

Get this problem set

Download everything (Topic5.2-Dimensionality_Reduction_Problems.zip) — the notebook and tep_monitoring.csv, in a folder that is ready to run as-is.

The Download badge at the top of the page will also give you the notebook on its own. On Vocareum everything is already set up for you.

Before you start

This problem set accompanies Dimensionality Reduction. It is worth 100 points:

  • Part A — Skill Checks (30 pts) — short answers, auto-graded, resubmit as often as you like until they pass.

  • Part B — Visualization (35 pts) — plots plus written interpretation, peer graded.

  • Part C — Open Ended (35 pts) — one synthesis problem, peer graded.

The parts build on each other: Part A works out the syntax you need for Part B, and Part B produces the evidence you argue from in Part C. Do them in order.

Setup#

In this problem set we will work with the Tennessee Eastman Process (TEP), a simulated chemical plant that Eastman Chemical published in 1993 as a benchmark for process monitoring. The control system records 52 variables every 3 minutes: 22 flows, pressures, levels and temperatures (XMEAS_01–XMEAS_22), 19 analyzer compositions (XMEAS_23–XMEAS_41) and 11 valve positions (XMV_01–XMV_11). The full list is in the benchmark’s documentation. Two variables come up by name below: XMEAS_02 D feed and XMEAS_03 E feed, both in kg/h.

One of the most common uses of PCA in chemical plants is process monitoring. The idea is to fit PCA to data from normal operation and then check each new sample to see whether it still looks like normal operation. The file data/tep_monitoring.csv is set up to do exactly that:

run

split

rows

what it is

normal_train

train

500

normal operation, 25 h. Fit everything to this.

fault03_train, fault09_train, fault15_train

train

480 each

faults 3, 9, 15 active throughout (Part C only)

fault01_test, fault04_test, fault05_test

test

480 each

8 h normal, then the fault switched on for 16 h

fault03_test, fault09_test, fault15_test

test

480 each

the same, for faults 3, 9, 15 (Part C only)

The faulty column is 1 on rows where the fault is active, fault is the fault number (0 for normal), and time_h counts hours within each run. In every test run the fault switches on at 8 h.

fault

disturbance

1

step change in the A/C feed ratio

4

step change in the reactor cooling water inlet temperature

5

step change in the condenser cooling water inlet temperature

3

step change in the D feed temperature

9

random variation in the D feed temperature

15

sticking condenser cooling water valve

The 52 variables have very different units, so we will work with standardized data throughout, using the mean and standard deviation of normal_train only. The standardize helper below does this. Use it for every set of rows (including fault rows), so that everything is measured relative to normal operation.

%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
try:
    plt.style.use('../settings/plot_style.mplstyle')   # available inside the book
except OSError:
    pass                                              # downloaded notebook: use defaults

df = pd.read_csv('data/tep_monitoring.csv')
sensors = [c for c in df.columns if c.startswith('XM')]
normal_train = df[df['run'] == 'normal_train']
test = df[df['split'] == 'test']

mu, sd = normal_train[sensors].mean(), normal_train[sensors].std()

def standardize(rows):
    """Scale rows of df by the normal-operation mean and standard deviation."""
    return ((rows[sensors] - mu) / sd).values

Z_train = standardize(normal_train)
print(f"{len(df)} rows, {len(sensors)} sensors;  Z_train is {Z_train.shape}")

4820 rows, 52 sensors;  Z_train is (500, 52)

Run this once. It defines the check helper used after each Part A answer. It will not run inside the book — use the downloaded notebook.

# Run this once. It defines `check`, used after each Part A answer below.
import hashlib

_EXPECTED = {
    'q1': (6, ['2a6a9d08764d', '3683be74603d', 'f673f0fffee3', '8291fe5714d2', '5fd8456a51b5', '4793940c6051', '5dc8b9350272']),
    'q2': (5, ['1d18f3b85e46', '45747019f833', 'fae2b5aa0328', '8b813f3c3651', '81b7f35b45fc', 'b2499f4b883d', '6fbac7d622e6', '2de5bdcea5e5', 'e231df4b8242']),
    'q3': (3, ['7f0603b0ae1c', '06ff91cd5970', '189acfbd91c9', '46a5ca376bf7', '114e48c95b81', '90e380245f82', 'bf13f51d682e', '15d78a49a9b1', '9e165e3c5161', 'fef88969682a', '3d719c7957a8', '4aa931833115', 'f4b73f0ebdd9', '006fb43ed30d', 'd5aa20459b8c', '954a08408650', '18fdff6d133b', '73dbd6429689', 'f43895e95b5b', 'd6aa31a6b9f0', '348ef55b6386']),
}

def check(qid, value):
    """Compare an answer against a stored digest. Practice only -- nothing is recorded."""
    places, digests = _EXPECTED[qid]
    if value is Ellipsis:
        print(f"{qid}: not answered yet")
        return
    got = hashlib.sha256(f"{float(value):.{places}f}".encode()).hexdigest()[:12]
    print(f"{qid}: {'correct' if got in digests else 'not the expected value'}")

Part A — Skill Checks (30 pts)#

Three questions, 10 points each.

A PCA monitor relies on one key assumption: normal operation only occupies a few directions in the 52-dimensional space, and a fault pushes samples out of those directions. In this part we will check that assumption by asking how many directions normal operation really uses, how many components it takes to describe it, and whether a sample that leaves those components actually gets caught.

A1. Directions the plant never moves in (10 pts)#

Exercise 1

Compute the covariance matrix of Z_train with np.cov and find its eigenvalues, as in the chapter. Count the eigenvalues that are smaller than \(10^{-3}\).

Assign the count to n_tiny.

# YOUR CODE HERE
n_tiny = ...

check("q1", n_tiny)

A2. How many components does normal operation need? (10 pts)#

Exercise 2

Using the A1 eigenvalues, find the smallest number of principal components whose cumulative retained variance is at least 0.90.

Assign that integer to k90.

# YOUR CODE HERE
k90 = ...

check("q2", k90)

A3. Does the monitor catch fault 5? (10 pts)#

A sample that looks like normal operation can be reconstructed well from the first \(k\) principal components, while a sample that does not will have a large reconstruction error. The squared prediction error (SPE) of a standardized sample \(\mathbf{z}\) is this error, summed over the 52 sensors:

\[ \text{SPE} = \lVert \mathbf{z} - \hat{\mathbf{z}} \rVert^2, \qquad \hat{\mathbf{z}} = P P^{\mathsf T} \mathbf{z} \]

where the columns of \(P\) (52 × \(k\)) are the first \(k\) eigenvectors. This is just the reconstruction from the chapter, written as a formula. The monitor raises an alarm when the SPE is above a limit, which we will set as the 99th percentile of SPE over normal_train.

Exercise 3

Using \(k = 5\) components:

  1. compute the SPE of every normal_train sample and set the limit with np.percentile(..., 99);

  2. compute the SPE of the rows of fault05_test where faulty == 1, standardized with standardize;

  3. find the fraction of those rows with SPE above the limit.

Assign this fraction (the detection rate) to det_f5.

# YOUR CODE HERE
det_f5 = ...

check("q3", det_f5)

Part B — Visualization (35 pts)#

In Part A we counted how many components standardized normal operation needs, and how often a 5-component monitor catches fault 5. However, these numbers don’t tell us what PCA actually found or when the monitor fires. In this part we will look at what the variance ranking picks out with and without standardization, and how the three faults look to the monitor over time and in the score plot.

Exercise 4

1. Compare scree curves. On one set of axes, plot the cumulative retained variance against the number of components for normal_train, once standardized and once with the raw, unstandardized sensors. Report the two sensors with the largest absolute loadings on the raw PC1, and explain in one sentence why two raw components keeping about 92% of the variance is not a useful compression.

2. Plot the monitoring chart. For fault01_test, fault04_test and fault05_test, plot the SPE (\(k = 5\)) against time_h on a logarithmic y-axis, with the A3 limit and a line at the fault onset (8 h). For fault 5, report the detection rate in the first 8 hours after onset and in the last 8 hours, and explain the A3 value in one sentence.

3. Plot the scores. Plot the first two principal component scores of normal_train and of the faulty rows of the three test runs, colored by fault. In one to two sentences, say which fault the SPE chart catches but the score plot does not show, and why that is possible.

Label all axes.

# YOUR CODE HERE

Part C — Open Ended (35 pts)#

So far the monitor works well: it catches faults 1 and 4 almost all the time, and fault 5 for as long as the fault is visible. However, the file has three more faults (3, 9 and 15), which are a step change and a random variation in the D feed temperature and a sticking condenser valve. These are real disturbances that matter to the plant.

In this part we will see whether these faults can be detected, and if not, whether the problem is with the method or with the data. Perhaps a different number of components would help, or perhaps PCA is not the right tool since it is linear and does not use the labels. On the other hand, maybe there just isn’t enough information in the 52 sensors. One way to find out is to give the problem every advantage we can and see whether anything changes.

Exercise 5

1. Sweep the number of components. For \(k\) from 2 to 45, compute the SPE detection rates for faults 3, 9, 15 and 5, along with the monitor’s false-alarm rate (the fraction of all test rows with faulty == 0 that are above the limit). Plot all five rates against \(k\), and report the largest value of (detection rate − false-alarm rate) over all \(k\) for each of faults 3, 9 and 15.

2. Use the labels. Fit LinearDiscriminantAnalysis to all of the split == 'train' rows, using fault as the class (0, 3, 9 or 15). Score it on the test rows of the three hidden-fault runs, labeling each row with its fault where faulty == 1 and 0 otherwise. Report the test accuracy, the training accuracy, and the accuracy of always predicting “normal”.

3. Method or data? In two sentences, say whether these faults are a limitation of the method or a property of the data, citing your numbers from tasks 1 and 2, and say what you would recommend to the plant’s engineers.

There is more than one defensible answer. You are graded on whether your numbers support your conclusion.

# YOUR CODE HERE

Additional Reading#

  1. L. H. Chiang, E. L. Russell and R. D. Braatz, Fault Detection and Diagnosis in Industrial Systems (Springer, 2001), doi:10.1007/978-1-4471-0347-9 — PCA monitoring with the T² and Q (SPE) statistics, on these data files.

  2. J. J. Downs and E. F. Vogel, “A plant-wide industrial process control problem”, Comput. Chem. Eng. 17, 245–255 (1993), doi:10.1016/0098-1354(93)80018-I.

  3. S. J. Qin, “Statistical process monitoring: basics and beyond”, J. Chemometrics 17, 480–502 (2003), doi:10.1002/cem.800.

Contents