神經網路系統鑑別

系統鑑別（system identification）是只用量測的輸入、輸出資料，建立受控體的動態模型。這一頁以模糊控制章的水槽為例，假裝不知道它的方程式，只記錄「閥開度 → 水位」的資料，比較線性 ARX 模型與神經網路 NARX 模型學得如何。

NARX 模型 Nonlinear AutoRegressive with eXogenous input

離散時間的動態系統，下一刻的輸出由過去的輸出與輸入決定。取過去 $n_y$ 個輸出、$n_u$ 個輸入當作迴歸向量：

$$\hat y(k+1) = f\big(y(k), \ldots, y(k-n_y+1),\ u(k), \ldots, u(k-n_u+1)\big)$$

$f$ 是線性函數時就是傳統的 ARX 模型：

$$\hat y(k+1) = a_1 y(k) + a_2 y(k-1) + b_1 u(k) + b_2 u(k-1) + c$$

係數可以用最小平方法一次解出。$f$ 用類神經網路表示時稱為 NARX，可以描述非線性的關係。本頁兩者都取 $n_y = n_u = 2$，也就是迴歸向量 $\mathbf x(k) = [y(k), y(k-1), u(k), u(k-1)]$。

模型有兩種用法，差別很重要：

一步預測（序列–平行模式，series-parallel）：輸入用量測到的過去輸出。訓練就是這樣做，把它變成一般的迴歸問題，誤差不會累積。
自由模擬（平行模式，parallel）：只給初始值，之後都回授模型自己的預測。拿模型來預測未來或設計控制器時就是這樣用，小誤差會一步步累積放大。

一步預測很準不代表模型好：任何模型只要學到「$y(k+1) \approx y(k)$」，一步預測的誤差就很小。判斷模型好壞要看沒看過的資料上的自由模擬。

實驗設計 Experiment Design

學得好不好，資料比網路更重要：

輸入要涵蓋整個操作範圍：閥開度每 10～60 秒隨機換一個 0～0.62 之間的值（隨機階梯，類似 APRBS 訊號），讓水位跑遍 0～2 m。網路只在資料出現過的範圍內可靠，範圍外是外插，常常錯得離譜。
要有不同的變化快慢：只有很慢的變化學不到動態，只有很快的變化學不到穩態。
訓練與測試分開：訓練用 3000 秒（50 分鐘），測試用另外隨機產生的 600 秒。

網路的兩個實用技巧：輸入先標準化（減平均、除以標準差），以及讓網路輸出變化量 $\Delta y = y(k+1) - y(k)$ 而不是 $y(k+1)$ 本身，也就是 $\hat y(k+1) = y(k) + \text{NN}(\mathbf x(k))$。水位每秒只變幾公分，直接學 $y(k+1)$ 時網路大部分的能力都花在「複製 $y(k)$」上。

範例：水槽的 ARX 與 NARX Example

網路為 4→32→32→1，tanh 激活函數，共 1,249 個參數，用 Adam 對整批資料訓練 3000 次（約 20 秒）：

# 系統鑑別：只用量測的「閥開度 u → 水位 h」資料，學出水槽（同模糊控制章）的動態模型
# 線性 ARX 與神經網路 NARX 比較；兩者都用 y(k+1) = f(y(k), y(k−1), u(k), u(k−1))
import json
import math
from pathlib import Path

import numpy as np
import torch
from torch import nn

def tank(u, h0=0.5):
    """真實系統（只用來產生資料）：dh/dt = 0.5u − 0.05√(2gh)，取樣 1 s，槽高 2 m"""
    h, out = h0, [h0]
    for uk in u:
        for _ in range(10):                                            # 每秒內 RK4 十步
            f = lambda x: 0.5 * uk - 0.05 * math.sqrt(2 * 9.8 * max(x, 0.0))
            k1 = f(h); k2 = f(h + 0.05 * k1); k3 = f(h + 0.05 * k2); k4 = f(h + 0.1 * k3)
            h = min(max(h + 0.1 / 6 * (k1 + 2 * k2 + 2 * k3 + k4), 0.0), 2.0)
        out.append(h)
    return np.array(out)

def steps_input(n, rng, hi=0.62):
    """隨機階梯輸入：每 10～60 s 換一個 0～0.62 的閥開度（涵蓋 0～2 m 的操作範圍）"""
    u = []
    while len(u) < n:
        u += [rng.uniform(0, hi)] * int(rng.integers(10, 61))
    return np.array(u[:n])

rng = np.random.default_rng(0)
u_tr = steps_input(3000, rng); y_tr = tank(u_tr)                        # 訓練資料：50 分鐘
u_te = steps_input(600, rng); y_te = tank(u_te)                         # 測試資料：另外 10 分鐘
X = lambda y, u: np.stack([y[1:-1], y[:-2], u[1:], u[:-1]], 1)          # [y(k), y(k−1), u(k), u(k−1)]
Xtr, Ytr = X(y_tr, u_tr), y_tr[2:]

# 線性 ARX：最小平方
A = np.c_[Xtr, np.ones(len(Xtr))]
theta = np.linalg.lstsq(A, Ytr, rcond=None)[0]
arx = lambda x: np.r_[x, 1.0] @ theta

# 神經網路 NARX：輸入先標準化，輸出水位的變化量 Δy
mu, sd = Xtr.mean(0), Xtr.std(0)
torch.manual_seed(0)
net = nn.Sequential(nn.Linear(4, 32), nn.Tanh(), nn.Linear(32, 32), nn.Tanh(), nn.Linear(32, 1))
opt = torch.optim.Adam(net.parameters(), lr=3e-3)
xt = torch.tensor((Xtr - mu) / sd, dtype=torch.float32)
dt_ = torch.tensor(Ytr - Xtr[:, 0], dtype=torch.float32).unsqueeze(1)
for epoch in range(3000):                                               # 一步預測的誤差（序列–平行模式）
    loss = nn.functional.mse_loss(net(xt), dt_)
    opt.zero_grad(); loss.backward(); opt.step()
narx = lambda x: x[0] + net(torch.tensor((np.asarray(x) - mu) / sd, dtype=torch.float32)).item()

def one_step(model, y, u):
    return np.array([model(x) for x in X(y, u)])

def free_run(model, y, u):
    """自由模擬（平行模式）：只給前兩個真實水位，之後都拿模型自己的預測回授"""
    p = list(y[:2])
    for k in range(1, len(u)):
        p.append(min(max(model([p[k], p[k - 1], u[k], u[k - 1]]), 0.0), 2.0))
    return np.array(p)

rm = lambda a, b: float(np.sqrt(np.mean((a - b) ** 2)))
print("ARX 係數 [a1 a2 b1 b2 c] =", np.round(theta, 4))
for name, m in (("線性 ARX", arx), ("神經網路 NARX", narx)):
    e1 = rm(one_step(m, y_te, u_te), y_te[2:])
    fr = free_run(m, y_te, u_te)
    print(f"{name}：測試資料的一步預測 RMSE {e1:.4f} m，自由模擬 600 s 的 RMSE {rm(fr, y_te[:len(fr)]):.4f} m")

layers = [[m.weight.tolist(), m.bias.tolist()] for m in net if isinstance(m, nn.Linear)]
fr_n, fr_a = free_run(narx, y_te, u_te), free_run(arx, y_te, u_te)
Path(__file__).with_name("tank_narx.json").write_text(json.dumps({
    "layers": layers, "mu": mu.tolist(), "sd": sd.tolist(), "arx": theta.tolist(),
    "test": {"u": u_te.tolist(), "y": y_te.tolist(), "narx": fr_n.tolist(), "arx": fr_a.tolist()}}), encoding="utf-8")

⬇ 下載 tank_narx.py（Python 3＋PyTorch）

實際執行的輸出：

ARX 係數 [a1 a2 b1 b2 c] = [ 1.8523e+00 -8.5380e-01  4.4240e-01 -4.3850e-01 -2.0000e-04]
線性 ARX：測試資料的一步預測 RMSE 0.0025 m，自由模擬 600 s 的 RMSE 0.1322 m
神經網路 NARX：測試資料的一步預測 RMSE 0.0019 m，自由模擬 600 s 的 RMSE 0.0147 m

結果（下載 tank_narx.json，內含網路權值，下一頁的 MPC 會用到）：

測試資料（600 s）	一步預測 RMSE	自由模擬 RMSE
線性 ARX（5 個係數）	0.25 cm	13.2 cm
神經網路 NARX（1,249 個參數）	0.19 cm	1.5 cm

一步預測兩者差不多，自由模擬差了約 9 倍：

為什麼線性模型不夠 Static Map

水槽的出水量和 $\sqrt h$ 成正比：$\dfrac{dh}{dt} = 0.5u - 0.05\sqrt{2gh}$。閥開度固定時，水位最後停在 $h = (u/0.1)^2 / 2g$，是一條拋物線。把模型的輸入固定、反覆套用到收斂，就得到模型的穩態特性：

ARX 的穩態是 $y = \dfrac{(b_1 + b_2)u + c}{1 - a_1 - a_2}$，只能是直線。而且由範例的係數，$b_1 + b_2 = 0.0039$、$1 - a_1 - a_2 = 0.0015$，穩態增益是兩個很小的數相除，係數的一點誤差就讓穩態差很多。這就是 ARX「一步預測準、自由模擬不準」的原因：每一步的偏差很小，但累積下去就跑到錯的穩態。NARX 學到了拋物線的形狀：水位 1 m 以下誤差約 2 cm 以內，水位越高誤差越大，1.5～1.8 m 時偏高約 5 cm。原因之一是高水位的資料較少：閥開度在 0～0.62 均勻分布，但水位和開度的平方成正比，訓練資料中水位超過 1.5 m 的時間只有 11%，低於 0.5 m 的卻有 46%。誤差雖小，下一頁會看到它仍會造成 MPC 的穩態誤差，需要另外修正。

自由模擬對照 Try It

下方動畫把範例訓練好的 NARX 權值與 ARX 係數匯出到網頁，和真實水槽並行，用同樣的閥開度做自由模擬。兩個模型都看不到真實水位，只能靠自己的預測往前推。可以用隨機輸入，或手動把閥開到很大、很小，看 ARX 在兩端偏多少：

自製動畫（依控制理論數值模擬）

實務注意事項 Practical Notes

直接用自由模擬的誤差訓練：把網路展開好幾步、對多步預測的誤差做反向傳播（類似遞迴網路的時間反向傳播），自由模擬會更準，但訓練較慢、較不穩定。
雜訊：實際的量測有雜訊，迴歸向量中的 $y(k)$ 也含雜訊，一步預測訓練會有偏差；常見的做法是加長預測步數訓練，或改用狀態空間形式的網路。
結構的選擇：$n_y, n_u$ 太小描述不了動態，太大容易過擬合；可以從系統的階數（水槽是一階）開始，用驗證資料比較。
不要外插：網路在資料範圍外沒有根據；用在控制時，要讓控制器知道模型可信的範圍（例如限制 MPC 的搜尋範圍）。

學到模型之後，下一頁用它來做模型預測控制。