--- jupytext: text_representation: extension: .md format_name: myst kernelspec: display_name: Python 3 language: python name: python3 --- ```{contents} :local: :depth: 2 ``` # Problems: Dimensionality Reduction :::{admonition} Get this problem set :class: seealso {download}`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. ::: :::{admonition} Before you start :class: tip This problem set accompanies {doc}`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](https://github.com/camaramm/tennessee-eastman-profBraatz#continuous-process-measurements). 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. ```{code-cell} ipython3 %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}") ``` 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': (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} :label: pr-eda-tep-tinyeig 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`. ::: ```{code-cell} ipython3 :tags: [skip-execution] # YOUR CODE HERE n_tiny = ... ``` ```{code-cell} ipython3 :tags: [skip-execution] check("q1", n_tiny) ``` ### A2. How many components does normal operation need? (10 pts) :::{exercise} :label: pr-eda-tep-k90 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`. ::: ```{code-cell} ipython3 :tags: [skip-execution] # YOUR CODE HERE k90 = ... ``` ```{code-cell} ipython3 :tags: [skip-execution] 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} :label: pr-eda-tep-spe5 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`. ::: ```{code-cell} ipython3 :tags: [skip-execution] # YOUR CODE HERE det_f5 = ... ``` ```{code-cell} ipython3 :tags: [skip-execution] 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} :label: pr-eda-tep-monitor **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. ::: ```{code-cell} ipython3 :tags: [skip-execution] # 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} :label: pr-eda-tep-hidden **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. :::: ```{code-cell} ipython3 :tags: [skip-execution] # 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](https://doi.org/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](https://doi.org/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](https://doi.org/10.1002/cem.800).