Problems: High-Dimensional Data

Problems: High-Dimensional Data #

Get this problem set

Download everything (Topic5.1-High_Dimensional_Data_Problems.zip) — the notebook and tep_sensors.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 High-Dimensional Data. 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 control and monitoring. The plant has a reactor, a condenser, a vapor–liquid separator, a recycle compressor and a stripper, and it converts reactants A, C, D and E into products G and H (plus a byproduct F). The control system records 52 variables every 3 minutes:

columns	what they are
`XMEAS_01` – `XMEAS_22`	continuous measurements: flows, pressures, levels, temperatures
`XMEAS_23` – `XMEAS_36`	on-line gas analyzers (reactor feed, purge), mol%, updated every 6 min
`XMEAS_37` – `XMEAS_41`	product analyzer, mol%, updated every 15 min
`XMV_01` – `XMV_11`	valve positions set by the control system, % open

A few variables come up by name below: XMEAS_01 A feed (kscmh), XMEAS_02 D feed (kg/h), XMEAS_04 A and C feed (kscmh), XMEAS_07 reactor pressure (kPa), XMEAS_09 reactor temperature (°C), XMV_03 A feed valve, XMV_10 reactor cooling water valve and XMV_11 condenser cooling water valve. The full list is in the benchmark’s documentation.

The file data/tep_sensors.csv contains four runs of the plant:

`run`	rows	what happened
`normal`	500	normal operation, 25 h
`fault01`	480	step change in the A/C ratio of the stream-4 feed
`fault04`	480	step change in the reactor cooling water inlet temperature
`fault05`	480	step change in the condenser cooling water inlet temperature

The fault is active from the first sample of each fault run. The fault column holds the fault number (0 for normal), and sample and time_h count samples and hours within each run.

%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
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_sensors.csv')
sensors = [c for c in df.columns if c.startswith('XM')]
normal = df[df['run'] == 'normal']
print(f"{len(df)} rows, {len(sensors)} sensor columns")
print(df['run'].value_counts().to_string())

1940 rows, 52 sensor columns
run
normal     500
fault01    480
fault04    480
fault05    480

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': (5, ['be7c0dfa4aab', 'b95fb28bf6d7', 'd69ec22f1c51', '69cde8e5b11b', 'dc1d1e308586']),
    'q2': (4, ['6321b7e4d08d', 'fc5211c92e16', '51344809cc64', '9eeb07d4ec37', '10a505f337d9', 'e8a782109980', 'baec74376bcf', 'cb7fb7ee3030', '6ba3c3386013', '4e47fe7e4d05', '0f06932e6e7e', '591477fa6e06', '4ec443878ff3', 'ea12c92f77cd', 'c181cc9fbe96', 'b8802875f1fb', 'bec14810b7ec', '2c48cabe8a11', '46447b154aed', '043e32d0eb7e', '420f2168e0d1']),
    'q3': (4, ['24b039797731', '36d4385ef209', 'a7d49d802832', 'c28b76801af1', '4c4c86d28f9b', 'b79281b7fbc0', 'f0a43b788226', '72564d568c08', '8293f5f9f380', 'ae8534a4b2d8', '894b2cf0d7b8', 'c1c104b80444', '4a423ba53119']),
}

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.

Before building any model on 52 sensors, it is worth asking a few basic questions about the sensors themselves. How many of them are just repeating each other? Does a given fault show up in many sensors, or only one? And what happens to the idea of “nearby” points once there are 52 axes instead of two?

A1. How many sensors repeat each other? (10 pts)#

Exercise 1

Using only the normal run, compute the 52 × 52 Pearson correlation matrix of the sensor columns. Count the pairs of distinct sensors with \(|r| > 0.9\). Note that \((i, j)\) and \((j, i)\) are the same pair and should only be counted once, and a sensor paired with itself does not count.

Assign the count to n_pairs_090.

# YOUR CODE HERE
n_pairs_090 = ...

check("q1", n_pairs_090)

A2. A fault with one fingerprint (10 pts)#

Exercise 2

Combine the rows of the normal and fault04 runs, and make a target that is 1 on the fault04 rows and 0 on the normal rows. Then compute the Pearson correlation between this target and each of the 52 sensors, similar to the feature–target correlations in the chapter.

Find the sensor whose correlation with the target has the largest magnitude, and assign that correlation (including its sign) to r_top_f4.

# YOUR CODE HERE
r_top_f4 = ...

check("q2", r_top_f4)

A3. What “nearby” means in 52 dimensions (10 pts)#

The chapter’s mean_nn_distance measures how far each point is from its nearest neighbor. On its own this number doesn’t have a natural scale, so we will compare it to the mean pairwise distance (the average distance over all pairs of distinct points):

\[ \rho = \frac{\text{mean nearest-neighbor distance}}{\text{mean pairwise distance}} \]

If \(\rho\) is near 0 then every point has close neighbors compared to the size of the cloud, and as \(\rho\) approaches 1 the nearest point is not much closer than a typical one.

Exercise 3

Standardize the 52 sensors of the normal run (subtract each column’s mean and divide by its standard deviation), then compute \(\rho\) for these 500 standardized points using the Euclidean distance over all 52 columns.

Assign it to nn_ratio.

# YOUR CODE HERE
nn_ratio = ...

check("q3", nn_ratio)

Part B — Visualization (35 pts)#

Part A gave us three numbers about the sensors, but numbers alone don’t tell us why. In this part we will use the standard EDA plots from the chapter to look at what the individual sensors look like, how they are tied to each other, and how the three faults look in the sensors that carry them.

Exercise 4

1. Plot the distributions. Plot a histogram of each of the 52 sensors in the normal run in a single grid of panels, with each panel titled by its column name. Report the number of distinct values that XMEAS_09 takes, and explain in one sentence why its histogram looks the way it does (the Setup table may help).

2. Plot the correlation structure. Plot the 52 × 52 correlation matrix of the normal run as a heatmap, with a diverging color scale fixed from −1 to 1. Using your A1 pairs, give one sentence with a physical reason why one group of strongly correlated sensors moves together.

3. Plot the faults. Apply your A2 method to fault01 and fault05 as well, and report the top sensor and its correlation for each of the three faults. Make a sns.pairplot of these three sensors over all four runs, colored by run. In one to two sentences, say which faults a single sensor separates from normal operation and which it does not.

Label all axes.

# YOUR CODE HERE

Part C — Open Ended (35 pts)#

Part B showed that the information in the 52 sensors is spread very unevenly. Some sensors repeat each other almost exactly, faults 1 and 4 each show up clearly in a single sensor, and fault 5 barely shows up anywhere. A3 also showed that with all 52 axes, normal samples are not very close to each other either. This raises a practical question: if we build an alarm for an operator, how many sensors should it watch? An operator dashboard only has room for a handful of sensors, but is that a compromise, or is it actually better than watching everything?

We will answer this using the simplest possible alarm rule, which only needs the summary statistics from the chapter. The range alarm learns each sensor’s minimum and maximum from the first 250 samples of the normal run, and a sample raises an alarm if any sensor on the panel falls outside its learned range. We can score a panel in two ways:

the false-alarm rate, the fraction of the last 250 normal samples that raise an alarm;
the detection rate for each fault, the fraction of that fault run’s samples that raise an alarm.

Exercise 5

1. Score the full panel. Implement the range alarm and report the false-alarm rate and the three detection rates for a panel that watches all 52 sensors.

2. Build a short panel. Choose a panel of at most six sensors. Use the fault–sensor correlations from Part B to pick sensors that carry each fault, and the A1 pairs to avoid spending two slots on sensors that repeat each other. Report your panel and its four rates.

3. Choose a panel. In two sentences, say which panel you would give the operator and why, citing your rates, and use A3 to explain the difference in false-alarm rate between the two panels.

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

# YOUR CODE HERE

Additional Reading #

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 — the process and its 20 disturbances.
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 — the source of these data files, and the standard reference for monitoring them.
The Braatz-group TEP data and simulation code — variable lists and the full set of 21 fault runs.

Problems: High-Dimensional Data

Contents

Problems: High-Dimensional Data#

Setup#