昨天結束的時候,我們手上有三樣東西:一個標定好的物理模型、一份誤差基準線,還有一份「物理抓不到什麼」的清單。
今天要用這份清單,把殘差模型建起來。
但在寫任何程式之前,有一件事必須先講清楚,因為它是這整篇最容易出事的地方。
殘差模型會偷走物理模型的方向性
先講我踩的坑。
我第一次做混合模型,結構寫成這樣:
OD_pred = f_物理(所有參數) + g_資料(所有參數 + 情境變數)
看起來很自然:物理模型給基準,殘差模型吃全部的資訊去修正。
結果訓練出來,整體 RMSE 從 0.030 mm 降到 0.011 mm,漂亮。
然後我做了一件事:把訓練好的混合模型拿去算偏微分——也就是 Day 4 那張「每個參數動一單位,外徑動多少」的敏感度表。
溫度的係數變成正的。
物理模型算出來是 −0.004 mm/°C,混合模型算出來是 +0.002 mm/°C。
原因是:g_資料 也吃了溫度。而歷史資料裡,溫度跟操作員的反應習慣糾纏在一起(Day 3 講過的那件事)。殘差模型是純統計的,它不知道什麼叫物理,它只會去擬合殘差裡任何看得到的相關性——包括那個假的、由操作習慣造成的正相關。
於是它學了一個跟物理相反的修正項,量級還比物理項大。
整體誤差降低了,但因果方向被毀掉了。
而這個模型的用途是什麼?是給師傅參數建議。一個方向相反的建議,比沒有建議危險得多。
所以殘差模型的輸入要被限制
修正的做法很簡單,但需要紀律:
物理模型已經描述的變數,不進殘差模型;或進去,但要受約束。
我最後的結構是把變數分三類:
A 類:物理變數。 螺桿轉速、線速、料溫、導體外徑。這些物理模型已經有明確的因果描述,不直接進殘差模型。
B 類:情境變數。 水槽溫度、模具已使用時數、本捲已跑長度、原料批次的 MFI、環境溫濕度、模頭壓力。這些是物理模型沒描述、或描述不足的東西,這才是殘差模型該吃的。
C 類:交互項。 有時候殘差確實跟「溫度 × 料批次」這種交互有關(不同批的料對溫度的反應不同)。這種可以進,但要明確設計、明確記錄理由,不要讓模型自己去挖。
分類的原則是一句話:
殘差模型只回答「同樣的物理參數下,為什麼今天和上週不一樣」,不回答「參數該往哪調」。
後面那個問題,永遠由物理模型回答。
為什麼用梯度提升樹,不用神經網路
這個場景裡,GBDT(LightGBM / XGBoost)幾乎總是比神經網路合適:
資料量級。 殘差建模的有效樣本可能只有幾萬到幾十萬筆(而且穩態段篩完會更少)。這個量級是 GBDT 的主場。
表格型異質特徵。 情境變數單位不同、尺度差很大、有類別型(料批次、模具編號)。樹模型不需要正規化,天然處理得好。
殘差通常是分段的、有閾值的。 「模具超過 200 小時之後開始偏」這種行為,樹很容易表示,MLP 要花很多參數去逼近。
可解釋性。 SHAP 在樹模型上是精確計算的。要跟現場解釋「為什麼建議這樣」,這很重要。
訓練成本。 每次換規格、換料號都要重訓或微調,GBDT 幾秒鐘的事。
神經網路值得考慮的時機,是當你要把「整段時序波形」當輸入(例如用 1D-CNN 抓螺桿脈動的波形特徵)的時候。但那是後面的事,現在不必要。
程式碼
python
import numpy as np
import pandas as pd
import lightgbm as lgb
from sklearn.model_selection import GroupKFold
PHYSICS_VARS = [ # A 類:物理模型負責,不進殘差模型
"screw_rpm", "line_speed", "melt_temp", "d_core",
]
CONTEXT_VARS = [ # B 類:物理沒描述的,殘差模型負責
"water_temp", # 收縮率的變異
"water_level",
"die_hours", # 模具磨損
"run_length", # 本捲已跑長度 → 積料
"hours_since_purge", # 清機後時數
"material_mfi", # 原料批次熔融指數
"material_age_days", # 料的存放天數 (Day 1 師傅提到的)
"ambient_temp",
"ambient_humidity",
"die_pressure", # 壓力異常 = 流道狀態改變
"head_temp_dev", # 模頭溫度與設定值的偏差
]
CATEGORICAL = ["die_id", "material_lot", "machine_id"]
def build_residual_dataset(df, popt, physics_fn):
X_phys = np.vstack([df[c] for c in
["screw_rpm", "line_speed", "d_core", "melt_temp"]])
df = df.copy()
df["od_phys"] = physics_fn(X_phys, *popt)
df["resid_mm"] = (df["od"] - df["od_phys"]) * 1000 # 目標:殘差 (mm)
return df
def train_residual_model(df, n_splits=5):
feats = CONTEXT_VARS + CATEGORICAL
X = df[feats].copy()
for c in CATEGORICAL:
X[c] = X[c].astype("category")
y = df["resid_mm"].values
groups = df["reel_id"].values # 關鍵:同一捲不可跨 fold
params = dict(
objective="regression",
metric="l1",
learning_rate=0.05,
num_leaves=31,
min_data_in_leaf=200, # 設大一點,避免學到單捲的雜訊
feature_fraction=0.8,
bagging_fraction=0.8,
bagging_freq=1,
lambda_l2=1.0,
verbose=-1,
)
oof = np.zeros(len(df))
models = []
gkf = GroupKFold(n_splits=n_splits)
for tr, va in gkf.split(X, y, groups):
ds_tr = lgb.Dataset(X.iloc[tr], y[tr], categorical_feature=CATEGORICAL)
ds_va = lgb.Dataset(X.iloc[va], y[va], reference=ds_tr)
m = lgb.train(params, ds_tr, num_boost_round=2000,
valid_sets=[ds_va],
callbacks=[lgb.early_stopping(100, verbose=False)])
oof[va] = m.predict(X.iloc[va])
models.append(m)
return models, oof, feats
分組交叉驗證:為什麼一定要用 reel_id 分組
這是另一個很容易出事的地方。
長度軸資料是高度自相關的——第 1,250.0 公尺和第 1,250.1 公尺的狀態幾乎一樣。如果隨機切 train/valid,同一捲的資料會同時出現在兩邊,模型只要「記住這捲的偏移量」就能得到極好的驗證分數。
我第一次這樣做,OOF MAE 是 0.004 mm。開心了半天,拿到新的一捲上去測,MAE 是 0.026 mm——幾乎等於沒修。
驗證集必須是模型完全沒見過的捲。 更嚴格的做法是用時間切分(拿最後兩週當測試),因為現實中模型面對的永遠是未來。
驗收:三個必過的關卡
殘差模型訓完,不能只看 RMSE。我固定跑三個檢查。
關卡一:方向性沒有被破壞
python
def check_monotonicity(phys, resid_models, base, feats):
"""
掃描物理參數,確認混合模型的敏感度方向與物理一致
"""
checks = {
"melt_temp": ("neg", np.arange(178, 195, 1.0)),
"line_speed": ("neg", np.arange(2.2, 3.8, 0.1)),
"screw_rpm": ("pos", np.arange(52, 74, 1.0)),
}
for var, (expect, grid) in checks.items():
ods = []
for val in grid:
b = dict(base); b[var] = val
od_p = phys.cold_od({k: b[k] for k in
["screw_rpm", "line_speed", "d_core", "T_melt"]})
# 殘差模型不吃物理變數 → 此處為常數修正
r = np.mean([m.predict(base_ctx_frame(b, feats))[0]
for m in resid_models])
ods.append(od_p * 1000 + r)
slope = np.polyfit(grid, ods, 1)[0]
ok = (slope < 0) if expect == "neg" else (slope > 0)
print(f" {var:<12} slope={slope:+.5f} expect={expect} "
f"{'PASS' if ok else '* FAIL ***'}")
因為 A 類變數不進殘差模型,這個檢查必定通過——這正是那個設計決策的回報。如果你讓物理變數進了殘差模型,這裡就是你會看到 FAIL 的地方。
關卡二:外插測試
拿一個訓練時完全沒出現過的規格(不同的 d_core、不同的目標外徑)去測。
python
def extrapolation_test(df_train, df_unseen, popt, physics_fn, resid_models, feats):
for name, d in [("in-domain", df_train), ("unseen spec", df_unseen)]:
Xp = np.vstack([d["screw_rpm"], d["line_speed"], d["d_core"], d["melt_temp"]])
od_p = physics_fn(Xp, *popt) * 1000
r = np.mean([m.predict(prep(d, feats)) for m in resid_models], axis=0)
for label, pred in [("physics only", od_p), ("hybrid", od_p + r)]:
mae = np.abs(pred - d["od"].values * 1000).mean()
print(f" {name:<12} {label:<13} MAE = {mae:.4f} mm")
理想的結果長這樣:
physics only hybrid
in-domain 0.030 0.011
unseen spec 0.034 0.015
重點不是 hybrid 比較好,而是hybrid 在沒看過的規格上沒有崩掉。
如果你跑純資料模型(把所有變數丟進 LightGBM 直接預測 OD),你會看到 in-domain 0.009(更好),但 unseen spec 0.08 以上(崩掉)。
這個差距就是物理骨架的價值。 它不是為了 in-domain 好看,它是為了 out-of-domain 不出事。
關卡三:殘差的殘差還有沒有結構
把 Day 5 的殘差診斷再跑一次,這次對 y - (物理 + 殘差模型)。
如果還有明顯結構,代表候選特徵沒放全,或者那個效應需要時序特徵(例如 rolling mean、滯後項)才抓得到。
SHAP:把殘差模型翻譯成人話
現場不會看 feature importance,但會看得懂「這捲為什麼偏」。
python
import shap
def explain_reel(model, X_reel, feats, top_k=4):
sv = shap.TreeExplainer(model).shap_values(X_reel)
mean_sv = sv.mean(axis=0)
order = np.argsort(-np.abs(mean_sv))[:top_k]
lines = []
for i in order:
d = mean_sv[i]
lines.append(f"{feats[i]}: {d:+.4f} mm")
return lines
輸出可以直接變成現場的一句話:
這捲物理模型算出來是 5.412 mm,實際偏大 0.021 mm。
主要來自:模具已使用 218 小時(+0.012)、本捲已跑 4,200 公尺(+0.006)、料批 MFI 偏低(+0.004)。
這句話比任何準確率數字都有說服力。 因為它說的是師傅本來就在想的事,只是把它量化了。
而且它天然帶出一個行動:模具該換了。
今天的產出
到這裡,孿生模型的第一個可用版本成形了:
OD_pred = f_物理(製程參數) + g_資料(情境變數)
↑ ↑
給方向與量級 解釋今天為什麼不一樣
可跨線移植 可解釋、可診斷
外插安全 隨資料持續更新
它能回答 Day 2 提到的那種反事實問題:「如果料溫再高兩度,外徑會變成多少」——由 f_物理 回答,而且方向保證正確。
它也能回答「為什麼今天和上週同樣參數跑出來不一樣」——由 g_資料 回答,並且指出是哪個情境變數造成的。
iThome鐵人賽
看影片追技術