--- jupytext: text_representation: extension: .md format_name: myst kernelspec: display_name: Python 3 language: python name: python3 --- ```{contents} :local: :depth: 2 ``` # Problems: High-Dimensional Data :::{admonition} Get this problem set :class: seealso {download}`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. ::: :::{admonition} Before you start :class: tip This problem set accompanies {doc}`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](https://github.com/camaramm/tennessee-eastman-profBraatz#continuous-process-measurements). 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. ```{code-cell} ipython3 %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()) ``` 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. ```{code-cell} ipython3 :tags: [skip-execution] # 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} :label: pr-eda-tep-redund 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`. ::: ```{code-cell} ipython3 :tags: [skip-execution] # YOUR CODE HERE n_pairs_090 = ... ``` ```{code-cell} ipython3 :tags: [skip-execution] check("q1", n_pairs_090) ``` ### A2. A fault with one fingerprint (10 pts) :::{exercise} :label: pr-eda-tep-onesensor 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`. ::: ```{code-cell} ipython3 :tags: [skip-execution] # YOUR CODE HERE r_top_f4 = ... ``` ```{code-cell} ipython3 :tags: [skip-execution] 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} :label: pr-eda-tep-nnratio **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`. ::: ```{code-cell} ipython3 :tags: [skip-execution] # YOUR CODE HERE nn_ratio = ... ``` ```{code-cell} ipython3 :tags: [skip-execution] 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} :label: pr-eda-tep-panels **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. ::: ```{code-cell} ipython3 :tags: [skip-execution] # 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} :label: pr-eda-tep-panel **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. :::: ```{code-cell} ipython3 :tags: [skip-execution] # YOUR CODE HERE ``` --- ## Additional Reading 1. 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](https://doi.org/10.1016/0098-1354(93)80018-I) — the process and its 20 disturbances. 2. 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](https://doi.org/10.1007/978-1-4471-0347-9) — the source of these data files, and the standard reference for monitoring them. 3. [The Braatz-group TEP data and simulation code](https://github.com/camaramm/tennessee-eastman-profBraatz) — variable lists and the full set of 21 fault runs.