化学プラントの反応器を監視しているとします。オペレーターが流量の設定値（setpoint）を上げれば、温度も圧力もそれに追従して動きます。これは正常です。ところが温度センサだけを見て「閾値を超えたら異常」とすると、操作のたびにアラートが鳴り、本物の異常はその誤報の海に埋もれます。
問題の根は、前回のサーベイ記事で見た「操作で動いた分を異常としない」をどう実現するか、です。
この記事が扱うのは、サーベイの (あ) 条件付き のうち、制御工学にもっとも素直な系統——残差ベース異常検知です。アイデアは一言でいえば「制御入力 $u$ を使って観測 $y$ を予測し、予測しきれなかった残差 $r=y-\hat{y}$ だけで異常を測る」。操作で説明できる変動は予測に吸収されて残差から消え、説明できない逸脱だけが残ります。
この考え方は、制御・信号処理の世界で モデルベースFDI（Fault Detection and Isolation、故障検知・分離） として40年以上磨かれてきました。Isermann、Gertler、Frank らが体系化した古典であり、いまも航空エンジン・化学プラント・上下水道・自動車の診断で現役です。本記事では、その核心を数式とコードで腹落ちさせます。
この手法が武器になる場面は明確です。
- 制御ループのあるプロセス監視: setpoint やバルブ開度が頻繁に変わるプラントで、操作起因の変動と故障を分離する
- 回転機械・エンジンの劣化検知: 回転数・負荷指令で大きく動く信号から、ゆっくり進む劣化だけを取り出す
- 制御系・ICS のセキュリティ監視: アクチュエータ指令とセンサ応答の整合性を見て、攻撃・故障を検出する
この記事の内容
- 残差の作り方の基本 — ARX（外生入力つき自己回帰）と最小二乗同定
- モデルベースFDIの3系統 — パラメータ推定・オブザーバ・パリティ空間（状態消去の幾何）
- 故障分離（FDIの I）— 構造化残差と故障シグネチャ行列で「どの故障か」を当てる
- 残差をどう評価するか — 白色性検定、CUSUM（逐次尤度比としての正体）、しきい値
- 閉ループ（フィードバック制御下）の難しさ — 操作量 $u$ も見て制御労力の増加を捉える
- 非線形への拡張 — NARXニューラルネット、grey-box RNN
- NumPy実装 — ARXを同定し、生信号ベースが操作で誤報する一方、残差CUSUMが劣化を捉える様子を実測
前提知識
1. 残差の作り方 — ARXで「操作への正常応答」を学ぶ
まず、システムが制御入力 $u$ にどう応答するかをモデル化します。線形でよく使われるのが ARX（AutoRegressive with eXogenous input、外生入力つき自己回帰） モデルです。
$$ \begin{equation} \hat{y}_t = \sum_{k=1}^{n_a} a_k\, y_{t-k} + \sum_{k=1}^{n_b} b_k\, u_{t-k} \end{equation} $$
「いまの出力は、過去の出力（自己回帰項 $a_k$）と過去の入力（外生項 $b_k$）の重み付き和で決まる」という素直なモデルです。$y$ は出力（センサ）、$u$ は制御入力（操作・コマンド）です。
係数 $\{a_k, b_k\}$ は、正常運転データから最小二乗で同定します。観測を回帰ベクトル $\varphi_t = [y_{t-1}, \dots, u_{t-1}, \dots]^\top$ とパラメータ $\theta = [a_1, \dots, b_1, \dots]^\top$ で書くと、
$$ \begin{equation} \hat{y}_t = \varphi_t^\top \theta, \qquad \hat{\theta} = \arg\min_\theta \sum_t \left(y_t – \varphi_t^\top \theta\right)^2 \end{equation} $$
これは通常の線形回帰なので閉形式（正規方程式）で解けます。こうして得たモデルは「この操作 $u$ をしたら、正常ならこの $y$ が出るはず」という予測器になります。
そして異常スコアの素になるのが残差（residual）です。
$$ \begin{equation} r_t = y_t – \hat{y}_t = y_t – \varphi_t^\top \hat{\theta} \end{equation} $$
正常運転中は、$\hat{y}_t$ が $y_t$ をよく説明するので $r_t$ は小さなランダム雑音になります。操作 $u$ がどれだけ大きく変動しても、それはモデルが予測するので残差には現れません——ここが「操作で動いた分を消す」の正体です。故障や劣化が起きると、モデルの予測と実際がずれ、残差が構造的に偏ります。
ARXは線形ですが、考え方は非線形にもそのまま拡張できます（後述のNARX）。次節では、残差を作る古典的な3つの流儀を整理します。
2. モデルベースFDIの3系統 — 残差をどう作るか
「正常モデルの予測と観測の差で異常を測る」という発想は、制御工学で3つの代表的な実装に結晶しています。
| 系統 | 残差の作り方 | 代表 |
|---|---|---|
| パラメータ推定 | モデル係数 $\theta$ を逐次推定し、正常値からのずれを残差に | Isermann |
| 状態推定（オブザーバ） | オブザーバ/カルマンで状態を推定し、イノベーション（予測誤差）を残差に | Frank |
| パリティ空間 | 入出力の有限窓から状態を消去し、冗長関係の破れを残差に | Gertler, Chow–Willsky |
前節のARXは「パラメータ推定」系統の一例です。3系統は実は無関係ではなく、残差を作るという一点で深くつながっています。実際、オブザーバ法とパリティ空間法は数学的に等価であることが知られています（Chow–Willsky, 1984）。まずもっとも実務で使われるオブザーバベースを掘り下げ、次にその等価物であるパリティ空間の幾何を見ます。
システムを状態空間で書きます。
$$ \begin{equation} x_{t+1} = A x_t + B u_t + w_t, \qquad y_t = C x_t + v_t \end{equation} $$
オブザーバ（ルーエンバーガ観測器、あるいはカルマンフィルタ）は、入力 $u$ と出力 $y$ から状態を推定します。
$$ \begin{equation} \hat{x}_{t+1} = A \hat{x}_t + B u_t + L\,(y_t – C\hat{x}_t) \end{equation} $$
ここで残差（イノベーション）は予測出力と実測の差です。
$$ \begin{equation} r_t = y_t – C\hat{x}_t \end{equation} $$
ゲイン $L$ で残差を観測器に戻し、推定をなめらかに補正します。正常なら $r_t$ は白色雑音（無相関なランダム信号）になり、故障が入ると残差の平均や分散、相関構造が偏ります。カルマンフィルタを使う場合、このイノベーションは理論上白色になるよう設計されており、イノベーションの白色性が崩れたら異常という判定が自然に導けます。
2.1 パリティ空間 — 状態を「消去」して冗長関係の破れを見る
オブザーバは状態 $\hat{x}$ を推定して残差を作りますが、パリティ空間法は発想が逆で、状態 $x$ を消去します。過去 $s$ ステップの出力と入力をまとめてスタックすると、状態空間モデルから次の関係が出ます。
$$ \begin{equation} Y = H_o\, x_{t-s} + H_u\, U + H_f\, f \end{equation} $$
ここで $Y=[y_{t-s},\dots,y_t]^\top$、$U=[u_{t-s},\dots,u_t]^\top$ を窓で積んだベクトル、$H_o$ は可観測性行列（状態が出力に効く係数）、$H_u$ は入力が出力に効く係数、$H_f$ は故障が出力に効く係数です。問題は未知の状態 $x_{t-s}$ が邪魔なこと。そこで $H_o$ の左零空間にあるベクトル $w$（$w^\top H_o = 0$ を満たす）を掛けると、状態の項が消えます。
$$ \begin{equation} r = w^\top (Y – H_u U) = w^\top H_f\, f \end{equation} $$
これがパリティ残差です。正常（$f=0$）なら $r=0$、故障があれば $r=w^\top H_f f \neq 0$。状態を推定せずに、入出力の冗長関係（パリティ関係）の破れだけを取り出します。
幾何的に見るとさらに腹落ちします。正常運転中の観測ベクトル $z=[U;Y]$ は、入出力の冗長関係が張る部分空間の上にぴったり乗ります。パリティ方向 $w$ はその部分空間と直交する方向で、正常な $z$ を $w$ に射影すると必ずゼロになる——これが「残差ゼロ」の正体です。
故障が起きると観測ベクトルが部分空間から外れ、その外れた分が $w$ への射影、すなわち残差として現れます。オブザーバ法（状態を推定）とパリティ空間法（状態を消去）は、見た目はまるで違うのに、最終的に同じ残差を生む——これがChow–Willskyらが示した等価性です。実装上は、低次の残差生成器を設計したいときはパリティ法で次数を決め、オブザーバ法でゲインを選ぶ、といった使い分けがされます。
2.2 故障分離（FDI の I）— 構造化残差で「どの故障か」を切り分ける
ここまでは「異常があるか（Detection）」だけでした。FDI の I＝Isolation（分離） は、さらに「どの部品の故障か」を当てます。鍵は、残差を1本でなく複数本作り、それぞれを特定の故障群にだけ感度を持つよう設計することです（構造化残差, Gertler）。
左の故障シグネチャ行列は「どの残差がどの故障に反応するか」を 0/1 で並べた設計図です。たとえば残差 $r_1$ はセンサBとアクチュエータの故障に反応するが、センサAの故障には反応しない（その故障成分が $w_1^\top H_f$ で消えるよう $w_1$ を選ぶ）、というように設計します。すると故障 $f_2$ が起きたとき、残差の発火パターンは $(r_1,r_2,r_3)=(1,0,1)$ という故障コードになり、これはシグネチャ行列の $f_2$ 列にしか一致しません。コードを読めば「故障はセンサBだ」と分離できます。
Gertlerが整理した設計条件は明快です——全故障を検出するには、どの故障の列もすべて0であってはならない（どれかの残差が必ず反応する）。故障を分離するには、すべての故障の列（コード）が互いに異なる必要があります。コードが重複する故障同士は区別できません。残差ベース異常検知が単なる「異常／正常」を超えて診断（diagnosis）に踏み込めるのは、この構造化残差の設計があるからです。
3系統はいずれも「残差を作る」ところまでが仕事です。作った残差をどう評価して異常と判定するかが次のテーマです。
3. 残差の評価 — 白色性検定とCUSUM
残差 $r_t$ を手にしたら、「これは正常な雑音か、それとも故障の兆候か」を判定します。2つの代表的な道具を見ます。
3.1 白色性 — 残差に構造が残っていないか
正常なら残差は白色、すなわち時間的に無相関なはずです。残差の自己相関関数（ACF）を見て、ラグ0以外で相関が立っていれば「予測しきれていない＝モデルと実機がずれている＝故障」と判断できます。
左（正常）はラグ0以外ほぼゼロで白色。右（異常）はラグが立ち、残差に構造が残っています。Ljung–Box 検定などで定量化できます。
3.2 CUSUM — 微小な平均シフトを積み上げる
故障は最初から大きく出るとは限りません。ゆっくり進む劣化では、残差の平均がわずかにずれます。一点ずつ見ると雑音に埋もれて見えませんが、累積和（CUSUM, Cumulative Sum） で積み上げると確実に検出できます。
$$ \begin{equation} S_t = \max\!\left(0,\ S_{t-1} + (r_t – k)\right) \end{equation} $$
$k$ は許容するドリフト（参照値）で、$S_t$ が閾値 $h$ を超えたら異常とします。
上段の残差では時刻380からの微小シフトはほとんど見えませんが、下段のCUSUMは積み上げで閾値を超え、明確に検出します。「小さいけれど持続するずれ」を捉えるのがCUSUMの真価です。
なぜこの単純な累積和が最適なのでしょうか。CUSUM は実は逐次確率比検定（SPRT）の片側版として導かれます。正常分布 $p_0(r)$（平均 $\mu_0$）と故障分布 $p_1(r)$（平均 $\mu_1=\mu_0+\delta$）を考えると、各時刻が故障寄りか正常寄りかは対数尤度比 $s_t = \log\frac{p_1(r_t)}{p_0(r_t)}$ で測れます。ガウスの場合これは $s_t \propto (r_t – k)$、ただし参照値 $k=\frac{\mu_0+\mu_1}{2}$ になり、CUSUM の更新式の中身 $(r_t-k)$ そのものです。
$$ \begin{equation} s_t = \log\frac{p_1(r_t)}{p_0(r_t)} = \frac{\delta}{\sigma^2}\left(r_t – \frac{\mu_0+\mu_1}{2}\right) \end{equation} $$
左図の2つのガウスが正常と故障の残差分布で、参照値 $k$ はちょうど2つの交点に当たります。右図が肝心です。尤度比をそのまま足し続けると（灰色）正常区間でも乱歩のように下へドリフトしてしまい、いつ故障が始まったか判別できません。CUSUM は $\max(0,\cdot)$ で負になったら0にリセットするため（赤）、正常区間ではゼロ付近に張り付き、故障が始まった瞬間から積み上がり始めます。このリセットこそが「正常がどれだけ続いても誤検知せず、故障が来た瞬間から素早く反応する」性質を生みます。閾値 $h$ と参照値 $k$ は、検出したい最小シフト幅 $\delta$ と許容誤報率から設計できます。
これで残差を作り評価する一通りが揃いました。しかし実システムには厄介な落とし穴があります——フィードバック制御です。
4. 閉ループの難しさ — 制御が異常を隠す
実際のプラントの多くは、目標値に追従するようフィードバック制御がかかっています。ここに残差ベース監視の最大の難所があります。
故障でプラントの振る舞いが変わっても、制御器がそれを打ち消すように操作量 $u$ を調整してしまうため、出力 $y$ には異常が現れにくくなります。出力だけ監視していると、制御がうまく隠した故障を見逃します。
対策は、出力 $y$ だけでなく操作量 $u$ も一緒に監視することです。制御器が「いつもより頑張って $u$ を動かしている」こと自体が異常の兆候になります。これを定量的に確かめましょう。1次プラントをPI制御で目標値に追従させ、途中でプラントの実効ゲインを $1.0\to0.55$ に落とす故障を入れます。
上段の出力 $y$ は故障の前も後も目標値にぴたりと追従しており、$y$ だけ見ても故障開始（赤線）はまったく分かりません——制御がきれいに隠しています。ところが下段を見ると、同じ目標値を保つのに必要な制御労力 $|u|/|r^*|$ が故障後に明確に段差を作って増えます（実測で $0.78\to1.12$、約1.4倍）。ゲインが落ちた分、制御器は同じ出力を出すのに $u$ をより大きく振らねばならないからです。出力は正常に見えても、その裏で制御器が余計に働いている——この「努力の増加」が故障の指紋になります。
実装としては、残差を $y$ と $u$ の同時モデルで作る、あるいは $u$ の統計量（振幅・制御労力）も併せて見る——制御系の異常検知では「センサとアクチュエータ（指令）をペアで見る」のが鉄則です。これは ICS セキュリティ（SWaT のような水処理テストベッドでセンサとアクチュエータ状態の整合性を見る研究）とも通じます。
次は、線形ARXでは表せない非線形システムへの拡張です。
5. 非線形への拡張 — NARXニューラルネットとgrey-box
実システムの多くは非線形です。線形ARXの予測精度が頭打ちなら、予測器を非線形にします。
NARX（Nonlinear ARX）ニューラルネットは、過去の入出力 $u_{t-1}, \dots, y_{t-1}, \dots$ を入力に、$\hat{y}_t$ を出力する非線形関数をニューラルネットで学習します。LSTM や TCN（Temporal Convolutional Network）を予測器にすれば、長い時間依存も扱えます。残差の作り方・評価の仕方は線形のときと同じで、予測器だけ差し替わるイメージです。実際、上下水道プロセスのような大規模産業システムを、データ駆動のNARXで監視する研究も報告されています。
grey-box（灰色箱）RNNは、既知の物理（線形ARXや微分方程式）とニューラルネットを組み合わせ、「物理で説明できる分は物理、残りをNNで補う」構成です。エンジン故障診断で物理ベースのgrey-box RNNで残差を生成する研究があり、純粋なブラックボックスより解釈性と外挿性に優れます。
出典: D. Jung, “Residual Generation Using Physically-Based Grey-Box Recurrent Neural Networks For Engine Fault Diagnosis,” arXiv:2008.04644, 2020, Fig. 2.
論文の原図は、本記事で繰り返し述べた残差生成の骨格をそのまま示しています。実システム（入力 $u_t$、故障 $f_t$）が出力 $y_t$ を生み、並走するモデル $\dot{\hat{x}}=g(\hat{x},u),\ \hat{y}=h(\hat{x},u)$ が予測 $\hat{y}_t$ を出し、その差 $r_t = y_t – \hat{y}_t$ が残差です。予測器が線形ARXでも非線形NNでも、この「観測とモデル予測の差を取る」骨格は不変だ——という本記事の主張を、論文の図が端的に裏づけています。
ポイントは、予測器がどれだけ高度でも、異常検知の枠組みは「残差を作って評価する」で一貫していること。ここまでの理屈を、最後に実データで確かめます。
6. NumPy実装 — 操作で動く信号から劣化だけを取り出す
ARXを同定し、「生信号を見ると操作で誤報するが、残差を見ると劣化を捉えられる」ことを実測します。
設定：2次の正常プラントを、ステップ状に切り替わる制御入力 $u$（setpoint風）で駆動します。後半でプロセスゲインが緩やかにドリフトする劣化（操作では説明できない故障）を入れます。
import numpy as np
rng = np.random.default_rng(0)
N = 4000
# 制御入力(操作): ステップ状に切り替わる setpoint 風
u = np.zeros(N); t = 0
while t < N:
L = rng.integers(120, 300); u[t:t+L] = rng.uniform(-1.5, 1.5); t += L
# 2次の正常プラント + 後半でゲインがドリフト(劣化)
a1, a2, b1, b2 = 1.2, -0.45, 0.5, 0.25
fault_start = 2600
yf = np.zeros(N); drift = np.zeros(N)
drift[fault_start:] = np.linspace(0, 0.6, N - fault_start)
for k in range(2, N):
gain = b1 + (drift[k] if k >= fault_start else 0.0)
yf[k] = a1*yf[k-1] + a2*yf[k-2] + gain*u[k-1] + b2*u[k-2] + rng.normal(0, 0.05)
正常区間（前半）だけでARX係数を最小二乗同定します。故障区間のデータは同定に使いません（正常モデルを学ぶのが目的）。
idx = np.arange(2, fault_start - 100)
Phi_tr = np.stack([yf[idx-1], yf[idx-2], u[idx-1], u[idx-2]], 1)
theta, *_ = np.linalg.lstsq(Phi_tr, yf[idx], rcond=None) # [a1,a2,b1,b2] を推定
# 全区間の1ステップ予測残差
kk = np.arange(2, N)
Phi_all = np.stack([yf[kk-1], yf[kk-2], u[kk-1], u[kk-2]], 1)
resid = yf[kk] - Phi_all @ theta
label = (kk >= fault_start).astype(int)
3つの検知器を比べます。(1) 生信号の偏差（移動平均からのずれ）、(2) 残差の大きさ、(3) 残差のCUSUM。
win = 50
ma = np.convolve(yf[kk], np.ones(win)/win, mode="same")
raw_score = np.abs(yf[kk] - ma) # 生信号ベース
res_score = np.convolve(np.abs(resid), np.ones(win)/win, mode="same") # 残差の大きさ
sd = resid[:fault_start-200].std() # 残差CUSUM
S = np.zeros(len(resid))
for i in range(1, len(resid)):
S[i] = max(0, S[i-1] + abs(resid[i]) - sd) # 参照値=正常残差の標準偏差
def auc(score, y):
order = np.argsort(score); ranks = np.empty(len(score)); ranks[order] = np.arange(1, len(score)+1)
npos = y.sum(); nneg = len(y) - npos
return (ranks[y==1].sum() - npos*(npos+1)/2) / (npos*nneg)
for name, s in [("生信号の偏差", raw_score), ("残差の大きさ", res_score), ("残差CUSUM", S)]:
print(f"{name}: AUC={auc(s, label):.3f}")
実行結果（seed固定）：
生信号の偏差: AUC=0.528
残差の大きさ: AUC=0.955
残差CUSUM: AUC=0.962
上図がすべてを物語ります。最上段の観測 $y$（青）は制御入力 $u$（灰）のたびに大きく動き、故障開始（赤線）がどこなのか目視ではまったく分かりません。中段の生信号ベースのスコアは、操作のたびに山ができて誤報だらけ。最下段の残差は、操作で動いた分がきれいに消え、故障開始後に偏りが出始めます。
数字でも明快です。
- 生信号の偏差＝0.528：ほぼランダム（0.5）。操作による変動と故障による変動を区別できず、検知器として機能していません。
- 残差の大きさ＝0.955：操作分を予測で消したことで、劣化がくっきり残差に出ます。
- 残差CUSUM＝0.962：微小なドリフトを積み上げ、さらに安定して検出します。右図の通り、故障開始後にCUSUMが明確に立ち上がります。
教訓は明快です。同じ信号・同じ閾値でも、「操作で動いた分を予測で差し引いたか」だけで、検知器はランダム同然から実用レベルへ変わる。残差ベースの神髄は、高度なモデルではなく「操作を予測に織り込む」という設計思想そのものにあります。
7. まとめ
制御入力・コマンドで動くシステムの異常検知を、モデルベースFDIの残差という古典的かつ強力な枠組みで見てきました。
- 核心：制御入力 $u$ で観測 $y$ を予測し、残差 $r=y-\hat{y}$ で異常を測る。操作で説明できる分は残差から消える。
- 残差の作り方：ARX/NARX（パラメータ推定）、オブザーバ／カルマン（状態推定）、パリティ空間（状態を消去して冗長関係の破れを見る）。3系統は等価で、いずれも「正常モデルの予測と観測の差」。
- 故障分離：残差を構造化し、故障シグネチャ行列の発火パターン（故障コード）で「どの故障か」を切り分ける（FDIの I）。
- 残差の評価：白色性検定（構造が残っていないか）と CUSUM（逐次尤度比検定の片側版。微小な持続シフトを検出）。
- 落とし穴：閉ループ制御は故障を隠すため、出力 $y$ と操作量 $u$（制御労力）をペアで監視する。
- 拡張：非線形は NARX ニューラルネットや grey-box RNN で予測器を強化。枠組みは不変。
残差ベースは「制御入力が観測できる」ときの第一選択です。一方、制御入力がコマンド（離散イベント）として与えられ、アクチュエータの状態遷移そのものを見る系統や、運用文脈との適合性を見る系統もあります。次の記事ではそれらに進みます。
主な参考文献
- R. Isermann, “Fault-Diagnosis Systems: An Introduction from Fault Detection to Fault Tolerance,” Springer, 2006.
- J. Gertler, “Fault Detection and Diagnosis in Engineering Systems,” Marcel Dekker, 1998.
- P. M. Frank, “Fault diagnosis in dynamic systems using analytical and knowledge-based redundancy,” Automatica, 26(3), 459–474, 1990.
- E. Y. Chow and A. S. Willsky, “Analytical redundancy and the design of robust failure detection systems,” IEEE Trans. Automatic Control, 29(7), 603–614, 1984.（パリティ空間）
- E. S. Page, “Continuous Inspection Schemes,” Biometrika, 41, 100–115, 1954.（CUSUM）
- Monitoring large-scale industrial systems with data-driven NARX, Computers & Chemical Engineering 系, 2025.
- Residual Generation Using Physically-Based Grey-Box Recurrent Neural Networks, arXiv:2008.04644, 2020.