神經網路模型預測控制

模型預測控制（model predictive control，MPC）每個取樣時刻都用模型預測未來，以數值最佳化找出最好的控制輸入。它不要求模型是線性的，所以很適合搭配上一頁學到的 NARX 網路。MPC 從 1970 年代起就是石化工業的主力控制方法，近年隨計算能力提升，也用在汽車、無人機與機器人上。

滾動時域 Receding Horizon

在時刻 $k$，從目前量測到的狀態出發，對未來 $N$ 步的輸入序列 $u(k), \ldots, u(k+N-1)$ 解一個最佳化問題：

$$\min_{u(k), \ldots, u(k+N-1)} \; J = \sum_{j=1}^{N} \big(r - \hat y(k+j)\big)^2 + \lambda \big(u(k) - u(k-1)\big)^2$$ $$\text{subject to}\quad \hat y(k+j+1) = f\big(\hat y(k+j), \hat y(k+j-1), u(k+j), u(k+j-1)\big), \qquad 0 \le u \le 1$$

第一項要水位接近設定值 $r$，第二項懲罰閥開度的劇烈變動。解出整個序列，但只執行第一步 $u(k)$；下一個取樣時刻量測新的水位，整個問題重解一次。因為每一步都用新的量測重新出發，模型的誤差與外來干擾會被回授修正，這就是 MPC 和「離線算好一條開迴路軌跡」的差別。

MPC 的優點：

限制條件（閥開度 0～1、水位上下限）直接寫在最佳化問題裡，不需要另外處理積分飽和。
非線性模型可以直接用，不必先線性化。
多輸入多輸出系統的寫法和單輸入一樣。

代價是每個取樣時刻都要解一個最佳化問題，計算量大，而且非線性模型的最佳化可能有多個區域最小值。

怎麼解：輸入分段 Move Blocking

預測 $N = 15$ 秒，有 15 個未知的輸入。常用的簡化是輸入分段（move blocking）：前 $m = 3$ 秒用 $u_1$、之後都用 $u_2$，未知數只剩兩個。兩個變數可以直接格點搜尋，保證找到全域最小值：

粗搜：$u_1, u_2$ 各取 0, 0.05, …, 1，共 21 × 21 = 441 組。
細搜：在最好的一組附近，各取間隔 0.005 的 11 個值，共 121 組。

每組都要用網路往前推 15 步，每秒共約 8,400 次網路計算。用 PyTorch 把 441 組排成一批同時計算，一個取樣時刻不到 10 毫秒。變數多時，則改用梯度法：網路可微分，用自動微分求 $\partial J / \partial u$，再交給最佳化器。

只用 $u_1$ 一段（整個區間閥開度固定）也可以，但試過的結果比 PI 還差：要快速加水時，固定的開度只能在「開大、之後會超越」與「開小、上升很慢」之間妥協。兩段讓 MPC 可以「先開大、再收回」。

偏差修正 Offset-Free MPC

模型一定有誤差。如果模型在某個水位預測水會繼續上升，實際上卻已經平衡，MPC 就會為了「抵消」不存在的上升而把閥關小一點，水位停在設定值下方，形成穩態誤差。範例的 NARX 在 1.5 m 平衡時，每秒會錯估 0.3 cm 的上升，15 秒的預測累積起來約 4 cm，結果 MPC 的穩態誤差約 −3 cm。

標準的解法是把誤差當成一個固定的干擾：模型和真實系統並行，用同樣的輸入自由模擬，得到模型輸出 $y_m(k)$，兩者的差

$$d(k) = y(k) - y_m(k)$$

就是模型目前的偏差。預測時假設它在未來保持不變，成本改為 $\sum (r - \hat y(k+j) - d(k))^2$，預測則從 $y_m$ 出發。到穩態時，模型輸出不再變化、預測曲線是平的，成本最小的條件就變成 $y_m + d = y = r$，所以不論模型在這個水位偏多少，都沒有穩態誤差。作用和 PI 的積分項相同。

範例：控制水槽水位 Example

水位設定值 0.5 → 1.5 → 0.8 → 1.8 m，各維持 100 秒。比較六種控制器：

PI：在 $h = 1$ m 線性化設計，$K_p = 0.578$、$K_i = 0.08$（閉迴路兩極點都在 −0.2），閥飽和時停止積分。
MPC 搭配四種模型：真實的水槽方程式（理想情況）、NARX、ARX；後兩者再各加偏差修正。

# 用學到的模型做模型預測控制（MPC）：控制水槽水位（同模糊控制章），和 PI 控制比較
# 需要先執行 tank_narx.py 產生 tank_narx.json（裡面有 NARX 網路與 ARX 係數）
import json
import math
from pathlib import Path

import torch

torch.set_default_dtype(torch.float64)
HERE = Path(__file__).parent
d = json.loads((HERE / "tank_narx.json").read_text(encoding="utf-8"))
W = [(torch.tensor(w), torch.tensor(b)) for w, b in d["layers"]]
mu, sd, theta = torch.tensor(d["mu"]), torch.tensor(d["sd"]), torch.tensor(d["arx"])

def true_model(y, ym1, u, um1):
    """真實的水槽（取樣 1 s，每秒 RK4 十步）；可一次算一批"""
    h = y.clone()
    f = lambda x: 0.5 * u - 0.05 * torch.sqrt(2 * 9.8 * x.clamp(min=0))
    for _ in range(10):
        k1 = f(h); k2 = f(h + 0.05 * k1); k3 = f(h + 0.05 * k2); k4 = f(h + 0.1 * k3)
        h = (h + 0.1 / 6 * (k1 + 2 * k2 + 2 * k3 + k4)).clamp(0, 2)
    return h

def narx_model(y, ym1, u, um1):
    h = (torch.stack([y, ym1, u, um1], -1) - mu) / sd
    for k, (w, b) in enumerate(W):
        h = h @ w.T + b
        if k < len(W) - 1:
            h = torch.tanh(h)
    return y + h[..., 0]

def arx_model(y, ym1, u, um1):
    return theta[0] * y + theta[1] * ym1 + theta[2] * u + theta[3] * um1 + theta[4]

def mpc(model, y, ym1, u_prev, r, N=15, m=3, lam=0.5, d=0.0):
    """閥開度在預測區間內分兩段：前 m 秒 u1、之後 u2。先在 0～1 間隔 0.05 試 21×21 組，再在最好的附近間隔 0.005 細搜。
    成本 J = Σ (r − ŷ − d)² + λ(u1 − u_前)²，只執行 u1（滾動時域）；d 為偏差修正量（量測值 − 模型輸出）"""
    def best_of(cand):
        u1, u2 = cand[:, 0], cand[:, 1]
        yk, ykm1, up = torch.full_like(u1, y), torch.full_like(u1, ym1), torch.full_like(u1, u_prev)
        J = lam * (u1 - u_prev) ** 2
        for j in range(N):
            u = u1 if j < m else u2
            yn = model(yk, ykm1, u, up).clamp(0, 2)
            J = J + (r - yn - d) ** 2
            ykm1, yk, up = yk, yn, u
        return cand[int(J.argmin())]
    g = torch.arange(21) * 0.05
    c = best_of(torch.cartesian_prod(g, g))
    off = (torch.arange(11) - 5) * 0.005
    return float(best_of(torch.cartesian_prod((c[0] + off).clamp(0, 1), (c[1] + off).clamp(0, 1)))[0])

ref = [0.5] * 20 + [1.5] * 100 + [0.8] * 100 + [1.8] * 100           # 設定值：0.5 → 1.5 → 0.8 → 1.8 m
u_ss = lambda h: 0.1 * math.sqrt(2 * 9.8 * h)                           # 維持水位 h 所需的閥開度

MODELS = {"MPC（真實模型）": true_model, "MPC（NARX）": narx_model, "MPC（ARX）": arx_model,
          "MPC（NARX＋偏差修正）": narx_model, "MPC（ARX＋偏差修正）": arx_model}

def run(controller):
    y, u_prev, us = [0.5, 0.5], u_ss(0.5), []
    I = u_ss(0.5) / 0.08
    ym = [0.5, 0.5]                                                     # 偏差修正：模型和真實系統並行，用同樣的輸入往前跑
    for k, r in enumerate(ref):
        if controller == "PI":                                          # PI：在 h = 1 m 線性化設計，兩極點都在 −0.2
            e = r - y[-1]
            u_raw = 0.578 * e + 0.08 * (I + e)
            u = min(max(u_raw, 0.0), 1.0)
            if u == u_raw:
                I += e                                                  # 閥飽和時停止積分（防積分飽和）
        elif "偏差修正" in controller:
            d = y[-1] - ym[-1]                                          # 把模型的誤差當成一個固定的干擾
            u = mpc(MODELS[controller], ym[-1], ym[-2], u_prev, r, d=d)
        else:
            u = mpc(MODELS[controller], y[-1], y[-2], u_prev, r)
        ym.append(float(MODELS.get(controller, true_model)(torch.tensor(ym[-1]), torch.tensor(ym[-2]), torch.tensor(u), torch.tensor(u_prev))))
        y.append(float(true_model(torch.tensor(y[-1]), None, torch.tensor(u), None)))
        us.append(u); u_prev = u
    return y[1:], us

out = {"ref": ref}
for name in ("PI", "MPC（真實模型）", "MPC（NARX）", "MPC（ARX）", "MPC（NARX＋偏差修正）", "MPC（ARX＋偏差修正）"):
    y, us = run(name)
    iae = sum(abs(a - b) for a, b in zip(y[1:], ref))
    over = max(y[k + 1] - ref[k] for k in list(range(20, 120)) + list(range(220, 320)))
    out[name] = {"y": y, "u": us}
    ss = [y[k] - ref[k - 1] for k in (120, 220, 320)]                   # 各段結束時的誤差（穩態誤差）
    print(f"{name:16s}：IAE {iae:5.2f} m·s，往上調時最大超越 {max(over, 0) * 100:4.1f} cm，"
          f"各段結束時的誤差 {", ".join(f"{e * 100:+.1f}" for e in ss)} cm")
(HERE / "tank_mpc.json").write_text(json.dumps(out, ensure_ascii=False), encoding="utf-8")

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

實際執行的輸出：

PI              ：IAE  6.33 m·s，往上調時最大超越  1.6 cm，各段結束時的誤差 +0.0, -0.0, +0.0 cm
MPC（真實模型）       ：IAE  5.29 m·s，往上調時最大超越  0.7 cm，各段結束時的誤差 +0.6, +0.0, -0.1 cm
MPC（NARX）       ：IAE  9.02 m·s，往上調時最大超越  0.0 cm，各段結束時的誤差 -2.9, -0.4, -2.6 cm
MPC（ARX）        ：IAE  4.70 m·s，往上調時最大超越  0.2 cm，各段結束時的誤差 -0.3, -0.4, +0.0 cm
MPC（NARX＋偏差修正）  ：IAE  5.83 m·s，往上調時最大超越  0.7 cm，各段結束時的誤差 +0.6, +0.1, -0.0 cm
MPC（ARX＋偏差修正）   ：IAE  5.57 m·s，往上調時最大超越  5.3 cm，各段結束時的誤差 +0.4, -0.2, +0.3 cm

需要先執行上一頁的 tank_narx.py 產生 tank_narx.json（或直接下載 tank_narx.json）。

控制器	IAE（m·s）	往上調時最大超越	各段結束時的誤差
PI	6.33	1.6 cm	0
模糊控制（模糊控制章的水槽 FIS）	8.69	1.0 cm	−1.0 cm
MPC（真實模型）	5.29	0.7 cm	+0.6 cm 以內
MPC（NARX）	9.02	0	−2.9、−0.4、−2.6 cm
MPC（NARX＋偏差修正）	5.83	0.7 cm	+0.6 cm 以內
MPC（ARX）	4.70	0.2 cm	−0.4 cm 以內
MPC（ARX＋偏差修正）	5.57	5.3 cm	+0.4 cm 以內

從這個比較可以看到：

MPC 本身確實比 PI 好：用真實模型時，上升更快（到達 90% 的時間 4 秒，PI 要 5～7 秒），超越不到 PI 的一半。PI 只在 1 m 附近設計，MPC 則在每個水位都依模型決定該開多大。
學到的模型＋偏差修正，幾乎達到真實模型的效果（IAE 5.83 對 5.29）。沒有修正的 NARX-MPC 有約 3 cm 的穩態誤差，原因如上一節所述。
ARX-MPC 的 IAE 反而最小，但不代表 ARX 是好模型。沒有偏差修正時，MPC 每一步都從量測的水位重新預測，只用到模型未來 15 秒的預測；ARX 每一步的偏差很小（在 1.5 m 平衡時只錯估每秒 0.03 cm），15 秒內還算準，這是運氣好的情形，它的穩態增益其實是錯的（見上一頁）。加上偏差修正後，需要模型自由模擬，ARX 自由模擬的誤差達 13 cm，修正量 $d$ 很大又一直在變，反而造成 5.3 cm 的超越。偏差修正假設的是「模型大致對、只差一個慢慢變的偏差」，模型太差時就不成立。
模糊控制不需要任何模型，但留有約 1 cm 的穩態誤差，這是水槽 FIS 規則不對稱造成的。
真實模型與 NARX＋修正的 MPC 在各段結束時還有約 +0.5 cm 的誤差，原因是細搜的格點：閥開度間隔 0.005 在 1.5 m 附近相當於約 2.7 cm 的穩態水位，MPC 在 0.540 與 0.545 兩個格點間來回切換（平衡需要 0.542）。縮小格點或改用連續最佳化即可改善。
設定值往下時（1.5 → 0.8 m），關閥後水從出水口流掉只要幾秒，差別在何時重新開閥：MPC 預先算出水位會降到哪裡，先全關一秒再逐步打開，幾乎沒有下衝；PI 沒有全關、反應較慢，下衝 1.3 cm；模糊控制全關太久，下衝約 8 cm。

即時的神經網路 MPC Try It

下方動畫在網頁裡即時執行同樣的 MPC：每個取樣時刻用範例訓練好的 NARX 做 562 組的格點搜尋（每秒約 8,400 次網路計算）。綠點線是這一刻 MPC 預測的未來 15 秒水位，可以看到每一步都在修正預測。移動設定值滑桿，比較有無偏差修正與 PI：

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

實務上的神經網路 MPC In Practice

計算時間：變數多、取樣快時，每步解最佳化太慢。常見做法是離線解很多次 MPC，再訓練一個網路直接模仿 MPC 的輸出（approximate / explicit MPC），線上只要一次網路計算，可以放進微控制器。
安全：網路模型在資料範圍外不可靠，最佳化器卻會專門找模型「以為很好」的地方，常常就是外插區。要限制搜尋範圍、加上狀態限制，必要時保留傳統控制器作為後備。
和強化學習的關係：MPC 用模型「每一步重新規劃」，強化學習則把策略學成一個固定的函數。結合兩者的 model-based RL（例如先學模型、再用 MPC 規劃）在機器人上很常見。