工場の制御室に並ぶ何十本ものトレンドグラフ。流量、温度、圧力、液位、弁の開度——どれも一見、いつもどおりに揺れています。ところが30分後、配管のどこかで小さな異変が起き、数十の計器がじわじわと連動して動き始めます。オペレータは「どの計器が、なぜおかしいのか」を限られた時間で見抜かなければなりません。これがプロセス監視と故障検知・診断（FDD: Fault Detection and Diagnosis）の現場です。
この「化学プラントの異変を多変量センサから当てる」課題を、誰もが同じ土俵で試せるように作られた仮想プラントが Tennessee Eastman Process（TEP） です。Downs と Vogel が1993年に発表して以来、プロセス監視・異常検知・故障診断の論文がほぼ必ず引き合いに出す、この分野の「MNIST」とも言える定番ベンチマークです。
TEP を実データで読み解くと、次の2つがはっきり見えてきます。
- 多変量プロセス監視とは何をする技術なのか — なぜ1本の計器ではなく34本まとめて見るのか、PCA と Hotelling T² が何を測っているのか
- 「検知できる故障」と「ほぼ検知できない故障」の差 — そして運転条件（モード）が変わると、なぜ正常まで誤検知してしまうのか
本記事では、ローカルに展開した Reinartz らの多モード版 TEP（6つの運転モード、各モード2900run × 600ステップ × 34変数、正常＋28種の故障）から本物の波形を描きながら、データ構造・故障の見え方・統計的監視の仕組みを順に追います。最後に PCA と Hotelling T² の管理図を Python で実装し、検知率と多モードの誤報問題を実データで確かめます。
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
Tennessee Eastman Process とは
まず全体像をイメージしましょう。TEP は、実在の化学プラント（Eastman Chemical 社のあるプロセス）を題材に、その動的挙動を忠実に再現したシミュレータです。実機で故障を意図的に起こすのは危険でお金もかかりますが、シミュレータなら何度でも安全に故障を注入し、計器がどう反応するかを記録できます。だからこそ「手法を公平に比べるベンチマーク」として定着しました。
プラントは5つの主要ユニットからなります。気体の原料（A・C・D・E）を反応器（Reactor）に送り込んで発熱反応で製品（G・H）を作り、凝縮器（Condenser）で冷やし、気液分離器（Separator）で気体と液体に分け、未反応ガスはコンプレッサ（Compressor）で昇圧して反応器へ再循環し、液はストリッパ（Stripper）で製品として取り出します。これらは互いに配管で結ばれ、しかも全体が閉ループの自動制御で運転されています。
この図が TEP の本質を凝縮しています。プラントの各所には流量・温度・圧力・液位・弁開度などを測る多数のセンサが付いており、本データではそれを 34個のプロセス変数として記録します。そこへ赤い矢印のように故障（IDV: Disturbance Variable）——原料組成のステップ変化、冷却水入口温度のランダム変動、弁の固着（sticking）など——が注入されます。重要なのは、プラントが閉ループ制御である点です。1か所で異変が起きても、制御系が他の変数を動かして打ち消そうとするため、異常の影響は複数の変数に連動して広がります。だから単独の計器ではなく、34変数の「いつもの関係」が崩れていないかを丸ごと監視する必要があるのです。
古典版の TEP は52変数（41の測定値 XMEAS と11の操作量 XMV）と20種類以上の既知故障 IDV1〜を持ちます。本データはその主要な34変数を使い、故障を28種類に整理した版です。では、その「いつもの揺れ」とは具体的にどんな波形なのか。正常運転を実データで見てみましょう。
正常運転の実波形：揺れているのに「いつもどおり」
下は Mode1（ある運転条件）の正常運転について、代表的な3変数を8run分重ねて描いたものです。1run は600ステップの時系列です。
3枚を見比べると、変数ごとに「揺れ方の個性」があることがわかります。上段（圧力に相当する変数#6、約2800kPa）と下段（ある流量に相当する変数#18）は、制御系の働きで設定値のまわりをランダムに細かく揺れています。これが正常です。一方、中段の変数#8（反応器温度に相当）は約122.9℃でほとんど一直線——実測すると標準偏差はわずか0.057℃で、ほぼ一定に保たれています。これは温度が極めて厳しく制御されていることを意味します。
ここに、多変量監視を難しくする最初のポイントがあります。正常運転でも各変数は常にゆらいでいるので、「値が普段と少し違う」だけでは異常とは言えません。どこまでが正常な揺れで、どこからが異常なのかを、データから定量的に線引きする必要があります。その線引きの主役が、のちほど登場する Hotelling $T^2$ です。
正常の姿がつかめたところで、いよいよ故障が注入されると波形がどう変わるのかを見ていきましょう。
故障を注入すると波形はどう動くか
各 run では、最初しばらくは正常に運転し、途中で故障が注入されます。下は故障クラス13を注入した run（赤）と正常な run（青）を、同じ変数#6（圧力に相当）で重ねたものです。
破線が故障の注入時刻（おおよそステップ29）です。それより前では赤も青も同じように約2800kPaのまわりを細かく揺れているだけで、見分けがつきません。ところが注入後、赤い故障 run は大きく波打ち始め、ときに2650kPa付近まで落ち込み、2870kPa近くまで跳ね上がります。正常の揺れ幅（青）をはっきり超えています。制御系が外乱を抑えようと弁を動かし続けるため、圧力が大きく振動するのです。
ここで「では1変数の振れ幅を見張れば十分では？」という疑問が湧きます。実は故障によっては、どの単独変数もほとんど動かず、変数どうしの関係だけがズレるものがあります。さらに厄介なのは、ほとんど何も変わらない故障が存在することです。検知の難しさが故障ごとにまるで違うことを、次に見てみましょう。
検知が容易な故障と、ほぼ検知不能な故障
TEP が良いベンチマークである理由のひとつは、簡単な故障と極めて難しい故障が混在していることです。下は「全34変数の標準化された逸脱度 $|z|$ の平均」を時系列で描き、容易な故障クラス6（左）と難しい故障クラス15（右）を、それぞれ正常の変動範囲（青帯）と重ねたものです。
左のクラス6では、赤い故障 run の逸脱度が正常の範囲を大きく突き抜けています。これは「誰が見ても異常」とわかる故障です。ところが右のクラス15では、赤い故障 run が青い正常の帯にすっぽり埋もれています。実測すると、クラス15の故障発生後の平均逸脱度は0.74で、正常時の0.726とほとんど変わりません。つまりデータの上では正常と区別がつかないのです。
この「ほぼ検知不能な故障」は偶然ではありません。古典 TEP では IDV3・IDV9・IDV15 が「微小でほとんど検知できない故障」として有名で、多くの手法がここで点を落とします。本データでもクラス3・9・15は検知率が極端に低く、あとで定量的に確かめます。簡単な故障だけで評価すれば、どんな手法も高得点に見えてしまう——TEP はわざと難物を混ぜることで、手法の実力を厳しく測っているのです。
1変数の逸脱では見えない故障があるなら、変数を「まとめて」見るとどうなるのか。34変数の動きを一枚の絵にして眺めてみましょう。
34変数を一望する：標準化ヒートマップ
下は、正常な1run（左）と故障クラス13の1run（右）について、横軸を時刻、縦軸を34個のプロセス変数とし、正常統計で標準化した値 $z$ を色で表したヒートマップです。赤が「普段より高い」、青が「普段より低い」を表します。
左の正常 run はほぼ一様に薄い色——どの変数も標準化値が $\pm 2$ 程度に収まり、際立った逸脱はありません。一方、右の故障 run は破線（注入時刻）を境に、複数の変数の行が同時に濃い赤・青に染まります。これが閉ループ・プロセスの故障の典型的な見え方です。1か所の異変が制御を通じて伝播し、多くの変数が連動して逸脱します。
このヒートマップは、なぜ多変量監視が有効かを直感的に語っています。個々の行（変数）だけ見ると微妙でも、「多くの変数が同時に普段と違う」というパターンは、まとめて見れば明確な異常信号になります。これを1つのスカラー指標に凝縮するのが Hotelling $T^2$ です。その前に、本データならではの難所——運転モードの違い——を押さえておきましょう。
運転モードが変わると「正常」の姿が変わる
同じプラントでも、作る製品のグレードや生産量を変えれば、運転点（操業条件）が変わります。本データはこの多モードを6種類（Mode1〜6）持っています。下は、ある変数（生産量・原料供給量に相当する変数#1）の正常分布を6モードで重ねたヒストグラム（左）と、代表4変数の正常平均を Mode1=1.0 として相対表示した棒グラフ（右）です。
左のヒストグラムが衝撃的です。変数#1の「正常値」は、モードによって約735から約5810まで、7倍以上も違うのです（実測の各モード平均は 3655 / 735 / 5179 / 5811 / 762 / 5810）。どのモードでもプラントは正常に動いているのに、計器が指す値はまったく別物です。右の棒グラフを見ると、変数によって挙動が分かれます。変数#19のように全モードでほぼ同じ（制御で設定値に固定）変数もあれば、変数#1・#24のようにモードごとに大きく変わる変数もあります。
ここに、現実のプロセス監視の核心的な難しさがあります。「正常」はひとつではありません。あるモードで覚えた正常の形を別モードに当てはめると、正常なのに「異常だ」と誤報してしまいます。これがマルチモード・プロセス監視、そしてドメイン適応が取り組む課題です。実際にこの誤報が起きる様子は、Hotelling $T^2$ を組んでから実測します。
ここまで波形と分布を眺めてきました。いよいよ、これらを定量的な「異常スコア」に落とし込む統計手法に入ります。
多変量管理図の考え方：PCA と Hotelling T²
なぜ単純な閾値ではダメなのか
34個の変数それぞれに「上限・下限」を引いて見張る素朴な方法を考えましょう。これには2つの欠点があります。第一に、変数が34個あれば、各々を1%の確率で誤検知するだけで、全体ではほぼ毎回どれかが警報を出してしまいます。第二に、変数どうしの相関を無視しています。たとえば「流量が増えれば温度も上がる」のが正常なら、流量だけ高くて温度が低い状態は、各々の上下限内でも関係としては異常です。
そこで、変数の相関構造ごと正常を学び、そこからの距離を測るのが多変量管理図です。中心となるのがマハラノビス距離、その2乗が Hotelling $T^2$ 統計量です。
Hotelling T² の定義
正常データの平均ベクトルを $\bm{\mu}$、共分散行列を $\bm{\Sigma}$ とします。観測ベクトル $\bm{x}$（34次元）に対する Hotelling $T^2$ は次式で定義されます。
$$ \begin{equation} T^2 = (\bm{x} – \bm{\mu})^{\top} \bm{\Sigma}^{-1} (\bm{x} – \bm{\mu}) \end{equation} $$
これは「正常の中心からの、相関を考慮した距離の2乗」です。$\bm{\Sigma}^{-1}$ が効いているおかげで、ばらつきの大きい方向の逸脱は割り引き、相関を破る方向の逸脱は強く罰します。普段一緒に動く変数がバラバラに動けば $T^2$ は跳ね上がります。
この $T^2$ を、PCA（主成分分析）を使って計算すると見通しが良くなります。正常データを標準化したうえで共分散行列を固有分解し、固有ベクトル $\bm{p}_k$（主成分の方向）と固有値 $\lambda_k$（その方向の分散）を得ます。観測を各主成分に射影したスコアを $t_k = (\bm{x}-\bm{\mu})^{\top}\bm{p}_k$ とすると、$T^2$ は次のように書き換えられます。
$$ \begin{equation} T^2 = \sum_{k=1}^{m} \frac{t_k^2}{\lambda_k} \end{equation} $$
各主成分方向の逸脱を、その方向の分散 $\lambda_k$ で割って正規化し、足し合わせているわけです。式(1)と式(2)は数学的に等価ですが、式(2)は「どの主成分方向で逸脱が起きたか」を分解して見られる利点があります。
T² と SPE：2種類の逸脱を見張る
実務の PCA 監視では、$T^2$ と並んで SPE（Squared Prediction Error, Q統計量） がよく併用されます。両者は監視する「方向」が違います。主成分を上位 $m$ 本だけ残して正常の挙動を表すモデルを作ると、観測ベクトルは「モデルが説明できる成分（主成分が張る部分空間）」と「説明できない残差成分」に分けられます。$T^2$ は前者の中で、正常な相関の範囲を超えた逸脱を測ります。一方 SPE は、残差成分のノルムの2乗
$$ \begin{equation} \text{SPE} = \lVert \bm{x} – \hat{\bm{x}} \rVert^2, \qquad \hat{\bm{x}} = \sum_{k=1}^{m} t_k \bm{p}_k \end{equation} $$
で、「正常では決して現れないはずの新しいパターン（モデル外の方向への飛び出し）」を捉えます。直感的には、$T^2$ は「いつもの動き方の範囲を超えて大きく動いた」、SPE は「いつもと違う動き方を始めた」を検知します。両者を同時に見張ることで、相関構造の崩れと新規パターンの出現の両方を取りこぼさずにすみます。本記事では話を一本化するため全変数版の $T^2$（マハラノビス距離）に絞りますが、実プラントの監視ではこの2枚看板が標準だと覚えておくと、論文の管理図がぐっと読みやすくなります。
管理限界（UCL）の引き方
$T^2$ がどこまで大きければ異常とみなすか、その境界を管理限界（UCL: Upper Control Limit）と呼びます。理論的には、正常がガウス分布なら $T^2$ は自由度（変数の数）の $\chi^2$ 分布に従うので、その上側1%点を UCL に選べます。ただし実データは厳密にはガウスでないので、本記事ではより素直に正常データの $T^2$ の経験的99%点を UCL に採ります。「正常なら100回に1回しか超えない線」を引く、という発想です。
直感がつかめたところで、まず正常 vs 故障を2次元で可視化し、それから $T^2$ の管理図を実装します。
Pythonで触る：データ読み込みと PCA 散布図
まず Mode1 を読み込み、正常データから PCA の方向を学んで、正常と数種の故障を第1・第2主成分の平面に射影してみます。
import numpy as np
import pickle, os
DS = os.path.expanduser("~/workspace/dataset/tep_reinartz")
def load_mode(m):
d = pickle.load(open(f"{DS}/TEPDataset_Mode{m}.pickle", "rb"))
return d["Signals"], d["Labels"] # (2900,600,34), (2900,)
S, y = load_mode(1)
normal = S[y == 0] # 正常run (100,600,34)
NV = S.shape[2]
# 正常統計で標準化（分散ゼロの定数変数は除外）
flat = normal.reshape(-1, NV)
mu, sd = flat.mean(0), flat.std(0)
keep = sd > 1e-6
mu_k, sd_k = flat[:, keep].mean(0), flat[:, keep].std(0)
def zk(X): # 有効変数のみ標準化
return (X[:, keep] - mu_k) / sd_k
# 正常の共分散を固有分解（PCA）
Z = zk(flat)
C = np.cov(Z, rowvar=False)
evals, evecs = np.linalg.eigh(C)
order = np.argsort(evals)[::-1]
evals, evecs = evals[order], evecs[:, order]
print("使用変数数:", keep.sum(), " 除外(定数):", np.where(~keep)[0].tolist())
これを実行すると、34変数のうち3つ（撹拌速度など常に一定の変数）が定数として除外され、31変数で監視モデルを組むことがわかります。
使用変数数: 31 除外(定数): [26, 30, 33]
次に、正常と故障クラス6・1（容易）・15（困難）を主成分平面へ射影します。
import matplotlib.pyplot as plt
P = evecs[:, :2]
pn = zk(normal[:20].reshape(-1, NV)) @ P
plt.scatter(pn[::5, 0], pn[::5, 1], s=6, alpha=.3, label="正常")
for cls in [6, 1, 15]:
Sc = S[y == cls][:20, 300:, :].reshape(-1, NV) # 発症後
pc = zk(Sc) @ P
plt.scatter(pc[::5, 0], pc[::5, 1], s=8, alpha=.4, label=f"故障{cls}")
plt.legend(); plt.show()
この散布図から、故障の「型」の違いがはっきり読み取れます。容易な故障6・1（橙・赤）は正常の青い雲から遠く離れた別の領域に飛び出しています。これらは主成分平面の上ですでに分離できる故障です。一方、困難な故障15（緑）は正常の雲とほぼ完全に重なっています。2次元では見分けようがありません。図4で見た「正常に埋もれる故障」が、PCA の視点でも同じく分離不能であることが確認できます。
散布図は2次元への射影なので情報が落ちます。31次元の情報を保ったまま1つのスコアに凝縮する $T^2$ を、次に実装しましょう。
Pythonで触る：Hotelling T² 管理図で故障を検知する
式(2)に従って $T^2$ を計算し、正常データから UCL を引きます。
rng = np.random.default_rng(0)
def t2(Zx): # 式(2): Σ t_k^2 / λ_k
sc = Zx @ evecs
return (sc ** 2 / evals).sum(1)
# 正常を学習/検証に分け、学習側の99%点をUCLにする
perm = rng.permutation(len(Z))
ntr = int(0.7 * len(Z))
T2_tr = t2(Z[perm[:ntr]])
UCL = np.percentile(T2_tr, 99)
fa = (t2(Z[perm[ntr:]]) > UCL).mean() # 検証正常の誤報率
print(f"UCL(99%)={UCL:.1f} 検証正常の誤報率={fa*100:.2f}%")
実行結果は次のとおりです。
UCL(99%)=51.9 検証正常の誤報率=0.93%
UCL は約52、そして学習に使っていない正常データでの誤報率は0.93%——狙いどおり「正常なら約1%しか超えない線」が引けています。この UCL を、故障 run と正常 run の $T^2$ 時系列に重ねます。
fr = S[y == 13][0] # 故障クラス13の1run
nr = normal[1] # 正常の1run
T2_f, T2_n = t2(zk(fr)), t2(zk(nr))
plt.plot(T2_n, label="正常run")
plt.plot(T2_f, label="故障クラス13 run")
plt.axhline(UCL, ls="--", color="k", label=f"UCL={UCL:.0f}")
plt.yscale("log"); plt.legend(); plt.show()
on = 29 # 故障注入時刻
print("発症後の超過率(故障run): %.0f%%" % (100*(T2_f[on:] > UCL).mean()))
print("超過率(正常run): %.2f%%" % (100*(T2_n > UCL).mean()))
発症後の超過率(故障run): 100%
超過率(正常run): 0.33%
この管理図が、多変量プロセス監視の全体像を一枚で示しています。正常 run（青）の $T^2$ は終始 UCL（破線、約52）の下を這い、超過率はわずか0.33%です。対して故障 run（赤）は、注入時刻（ステップ29）を境に $T^2$ が2桁も跳ね上がり、発症後は100%の時刻で UCL を超え続けます。縦軸が対数であることに注意してください。故障時の $T^2$ は数千〜数万に達し、正常の数十とは桁が違います。1つのスカラーに31変数の情報を凝縮するだけで、これほど明瞭に異常を検知できるのです。
ただし、ここまでの話には大きな前提が隠れています。「UCL を引いた運転モードと、監視するモードが同じ」という前提です。これが崩れると何が起きるか——多モードの誤報問題を実測しましょう。
多モードの誤報問題：別モードの「正常」を異常と誤る
Mode1 の正常で学んだ監視モデル（標準化の $\bm{\mu},\bm{\sigma}$、PCA、UCL）を、そのまま Mode2 の正常に当ててみます。
S2, y2 = load_mode(2)
norm2 = S2[y2 == 0]
T2_m1 = t2(zk(normal[2])) # Mode1 の正常
T2_m2 = t2(zk(norm2[0])) # Mode2 の正常（同じMode1モデルで採点）
plt.plot(T2_m1, label="Mode1 正常run")
plt.plot(T2_m2, label="Mode2 正常run")
plt.axhline(UCL, ls="--", color="k"); plt.yscale("log"); plt.legend(); plt.show()
fa2 = (t2(zk(norm2[:20].reshape(-1, NV))) > UCL).mean()
print("Mode2正常をMode1モデルで採点した誤報率: %.0f%%" % (100*fa2))
Mode2正常をMode1モデルで採点した誤報率: 100%
結果は劇的です。Mode2 はまったく正常に運転しているのに、Mode1 基準の $T^2$ は管理限界をはるか上空で振り切れ、誤報率は100%になります。図6で見たとおり、Mode2 の変数#1の正常値（約735）は Mode1（約3655）と桁違いなので、Mode1 の標準化では巨大な $z$ になり、$T^2$ が爆発するのです。
ここから現実的な教訓が得られます。「正常モデルは運転モードごとに用意する」「モードを跨ぐなら、モードを跨いでも崩れない特徴を学ぶ（ドメイン適応）」——本データが多モード構成になっているのは、まさにこの課題を研究するためです。あるモードで学んだ知識を別モードへ転移できるか、というのが多モード TEP の問いです。
最後に、28種の故障の検知しやすさを一望して、TEP というベンチマークの厳しさを定量化しましょう。
故障クラス別の検知難易度：易しい故障と難物の混在
各故障クラスについて、発症後（ステップ100以降）のサンプルが UCL を超えた割合を「検知率」として棒グラフにします。
det = []
for cls in range(1, 29):
Xc = S[y == cls][:, 100:, :].reshape(-1, NV)
det.append(100 * (t2(zk(Xc)) > UCL).mean())
det = np.array(det)
floor = 100 * (t2(zk(normal[:, 100:, :].reshape(-1, NV))) > UCL).mean()
order_c = np.argsort(det)[::-1]
plt.bar(range(28), det[order_c])
plt.axhline(floor, ls="--", color="k", label=f"正常の誤報水準≈{floor:.1f}%")
plt.legend(); plt.show()
print("易:", [(int(c+1), round(det[c],1)) for c in order_c[:5]])
print("難:", [(int(c+1), round(det[c],1)) for c in order_c[-5:]])
易: [(14, 100.0), (20, 100.0), (2, 100.0), (4, 100.0), (6, 100.0)]
難: [(28, 2.8), (16, 1.8), (15, 1.5), (5, 1.3), (21, 1.0)]
この棒グラフが TEP の難しさを物語ります。左側の緑のバー（クラス14・20・2・4・6・13など）は検知率100%で、Hotelling $T^2$ という古典手法でも完璧に捉えられます。ところが右端の赤いバー（クラス21・5・15・16・28、そして古典 TEP の難物 IDV3・9 もこの近辺）は、検知率がわずか1〜3%——破線で示した正常時の誤報水準（約1%）とほとんど変わりません。統計的にはこれらの故障は「正常と区別がつかない」のです。
この分布こそが、TEP が長年ベンチマークとして使われ続ける理由です。簡単な故障だけなら20年前の手法で十分です。研究者が競っているのは、この右端の難物をどこまで検知率の floor から引き上げられるかです。新しい深層学習手法や関連構造ベースの手法を提案するなら、まずこの「ほぼ検知不能な故障」で素朴な $T^2$ を上回れるかが試金石になります。評価指標そのものの落とし穴については VUS（Volume Under the Surface） や SMDとPSM：時系列異常検知ベンチマーク の記事もあわせてご覧ください。
まとめ
本記事では、化学プロセス異常検知の定番ベンチマーク Tennessee Eastman Process（TEP）を、実データの波形つきで読み解きました。
- 構造：反応器・凝縮器・分離器・コンプレッサ・ストリッパからなる閉ループ制御プラントを、34個のプロセス変数（流量・温度・圧力・液位・弁開度）で監視する。故障は途中で注入され、制御を通じて複数変数に連動して伝播する。
- 正常の姿：各変数は設定値のまわりで常に揺れている。反応器温度のように標準偏差0.057℃で固定される変数もあれば、大きく揺らぐ変数もある。「値が少し違う」だけでは異常と言えない。
- 故障の多様性：28種の故障には、誰が見てもわかるもの（検知率100%）と、データ上で正常と見分けがつかないもの（検知率1〜3%、IDV3・9・15が代表）が混在する。だから手法の実力を厳しく測れる。
- 統計的監視：PCA で相関構造を学び、Hotelling $T^2$ で「正常からの距離」を1スカラーに凝縮する。正常の99%点に管理限界を引けば、容易な故障は発症直後に100%検知できる（正常の誤報は約1%）。
- 多モードの壁：「正常」は運転モードごとに違う。Mode1 で引いた管理限界を Mode2 の正常に当てると誤報率100%。マルチモード監視・ドメイン適応が必要になる。
TEP は「化学プラントの異変を多変量センサから当てる」という現実的な課題を、安全に・再現可能に・公平に試せる貴重なテスト場です。異常検知の新手法を作ったら、まず簡単な故障で動くことを確かめ、次に難物（IDV3・9・15）と多モードの誤報という2つの壁にぶつけてみる——それが TEP の正しい使い方です。
次のステップとして、以下の記事も参考にしてください。