昨天寫了一個押出物理模型。今天要做一件很多人會跳過的事:
先量它有多不準,再決定要不要補。
跳過這一步直接接機器學習,你會得到一個看起來很準、但你不知道為什麼準的東西。而在製造現場,「不知道為什麼準」等於「隨時可能不知道為什麼不準」。
標定:那三個係數從哪來
Day 4 的模型裡有幾個待定係數:
係數 物理意義 大概範圍
k_screw 每 rpm 的產出量 (kg/hr) 由螺桿幾何決定,試車量
beta_swell die swell 的溫度修正 負值,量級 1e-3
shrink_rate 冷卻線收縮率 PE 0.015~0.03,PVC 較小
beta_rho 熔體密度溫度係數 查材料表,可固定
beta_rho 我直接查表固定住,不參與擬合——能從手冊查到的,就不要浪費資料去學。 這是混合模型的一個原則:自由度花在真的不知道的地方。
剩下三個用歷史資料擬合。
為什麼不能用一般迴歸
有人會問,為什麼不直接做多元線性迴歸把 rpm、線速、溫度丟進去。
因為 Day 4 那條式子是非線性的:
D_od = sqrt( d_core² + 4ṁ / (π ρ v) )
外徑對線速不是線性關係,是反比開根號。用線性迴歸擬合,在訓練範圍內看起來還好,一離開就崩。而且擬合出來的係數沒有物理意義,不能拿去驗證、不能跨產線移植。
所以要用非線性最小平方,直接擬合物理模型的參數。
資料準備
用 Day 3 對齊好的長度軸資料。但有幾件事要先做:
只取穩態段。 開機升速、換規格、停機重啟這些過渡段,物理模型(基於穩態質量守恆)根本不適用。硬塞進去只會把係數帶歪。
判斷穩態的簡單做法:滑動窗口內,線速和螺桿轉速的變異係數都低於閾值。
涵蓋足夠的操作範圍。 如果歷史資料裡線速永遠是 3.0,你擬不出線速的效應。實務上如果資料太集中,值得安排一次小型 DOE(設計實驗)——刻意在安全範圍內掃幾組參數,跑個幾百公尺。這幾百公尺的料是最便宜的一次學費,因為它換來的是一組可以重複使用的係數。
按規格分層。 不同料號、不同模具配置,k_screw 和 shrink_rate 都不一樣。一開始先針對單一主力規格做,不要一次想吃全部。
程式碼:標定與診斷
python
import numpy as np
import pandas as pd
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
def tag_steady_state(df, cols=("line_speed", "screw_rpm"),
win=200, cv_thresh=0.01):
"""
在長度軸資料上滑動,標記穩態段。
win: 窗口長度(格數),grid=0.1m 時 200 格 = 20 公尺
"""
mask = pd.Series(True, index=df.index)
for c in cols:
roll = df[c].rolling(win, center=True)
cv = roll.std() / (roll.mean().abs() + 1e-9)
mask &= (cv < cv_thresh)
# 線速為 0 (停機) 一律排除
mask &= (df["line_speed"] > 0.1)
return mask.fillna(False)
RHO_REF, T_REF, BETA_RHO = 940.0, 200.0, -6.5e-4 # 查表固定
def physics_od(X, k_screw, beta_swell, shrink):
"""
X: (4, n) array -> screw_rpm, line_speed, d_core(m), T_melt(°C)
回傳冷態外徑 (m)
"""
rpm, v, d_core, T = X
m_dot = k_screw * rpm / 3600.0 # kg/s
rho = RHO_REF * (1 + BETA_RHO * (T - T_REF))
A_ann = m_dot / (rho * v)
od_hot = np.sqrt(d_core**2 + 4 * A_ann / np.pi)
od_hot *= (1 + beta_swell * (T - T_REF))
return od_hot * (1 - shrink)
def calibrate(df):
X = np.vstack([df["screw_rpm"], df["line_speed"],
df["d_core"], df["melt_temp"]])
y = df["od"].values
p0 = [1.80, -1.2e-3, 0.020]
bounds = ([0.5, -1e-2, 0.000],
[5.0, 1e-2, 0.060]) # 物理上合理的範圍
popt, pcov = curve_fit(physics_od, X, y, p0=p0,
bounds=bounds, maxfev=20000)
perr = np.sqrt(np.diag(pcov))
names = ["k_screw", "beta_swell", "shrink_rate"]
for n, v, e in zip(names, popt, perr):
print(f" {n:<12} = {v:+.5g} ± {e:.2g}")
return popt
def report(y_true, y_pred, tol_mm=0.05):
err = (y_pred - y_true) * 1000 # mm
print(f" bias = {err.mean():+.4f} mm")
print(f" MAE = {np.abs(err).mean():.4f} mm")
print(f" RMSE = {np.sqrt((err**2).mean()):.4f} mm")
print(f" P95 |err| = {np.percentile(np.abs(err),95):.4f} mm")
print(f" within ±{tol_mm}mm = {(np.abs(err)<tol_mm).mean()*100:.1f}%")
return err
標定完之後,重點來了
擬合出來會得到一組係數和一個 RMSE。假設是 0.03 mm,規格公差是 ±0.05 mm。
看起來不錯。但誤差大小不是重點,誤差的結構才是。
一個好的物理模型,殘差應該是隨機的白噪聲。如果殘差有結構——隨某個變數系統性地變化、有趨勢、有週期——那代表有物理沒被抓到。
這是我認為整個標定流程裡最有價值的一步:
python
def diagnose_residual(df, resid_mm, candidates):
"""
candidates: 可能與殘差相關的變數名稱 list
殘差與某變數相關 → 該變數的物理效應沒被模型描述
"""
print("\n [residual correlation]")
rows = []
for c in candidates:
if c not in df: continue
r = np.corrcoef(df[c].values, resid_mm)[0, 1]
rows.append((c, r))
for c, r in sorted(rows, key=lambda x: -abs(x[1])):
flag = " <-- 有結構" if abs(r) > 0.2 else ""
print(f" {c:<18} r = {r:+.3f}{flag}")
return rows
def plot_residual(df, resid_mm, by):
"""殘差對某變數作圖,看有沒有非線性結構"""
fig, ax = plt.subplots(figsize=(7, 4))
ax.scatter(df[by], resid_mm, s=3, alpha=0.25)
ax.axhline(0, color="k", lw=1)
# 分箱平均,讓趨勢看得出來
b = pd.qcut(df[by], 20, duplicates="drop")
m = pd.Series(resid_mm).groupby(b.values, observed=True).mean()
ax.plot([i.mid for i in m.index], m.values, "r-o", ms=4, lw=1.5)
ax.set_xlabel(by); ax.set_ylabel("residual (mm)")
ax.set_title(f"Residual vs {by}")
return fig
殘差可能長出來的幾種樣子
我把幾種常見的結構和它的物理含意整理一下:
殘差隨 run_length(本捲已跑長度)單調上升
→ 通常是模具積料或模頭溫度漂移。押出機跑久了,模頭流道會結焦積料,有效流道變窄。這件事物理模型完全沒有描述。
殘差隨 hours_since_die_change 上升
→ 模具磨損。這是一個非常有價值的訊號,因為它可以直接轉成預測性維護:當殘差的長期趨勢超過閾值,就該換模了。用不著等品質超規。
殘差與 water_temp 明顯相關
→ 收縮率不是常數,它跟冷卻歷程有關。Day 4 我把 shrink_rate 當固定值,這是個已知的簡化。
殘差在不同 material_mfi(熔融指數)批次間有階梯
→ 料的批次差異。不同批的黏度不同,die swell 和實際流量都會變。這解釋了 Day 1 那位師傅說的「這批料放太久了」。
殘差呈現週期性震盪
→ 這個要小心。可能是螺桿的脈動(surging)、可能是上游絞線節距的週期性起伏(Day 2 那捲三公里的根因就是這個)、也可能是收線盤的張力波動。作法是對殘差做 FFT,看主頻對應到哪個轉動元件的週期。
殘差在低線速區偏大
→ 低速時熱傳導時間變長,冷卻更完全,穩態假設也更容易被破壞。可能需要分區間建模。
這一步真正的產出
很多人以為標定的產出是「一組係數」。
不是。標定真正的產出是這三樣:
一、一個有物理意義、可以跨線移植的基準模型。 k_screw 換一台機器要重標,但 shrink_rate 是材料的性質,可以帶著走。
二、一份「物理模型抓不到什麼」的清單。 這份清單直接決定了明天殘差模型的輸入特徵該放什麼。這比盲目把 50 個欄位丟進去做特徵選擇有效率得多,而且每個特徵都說得出理由。
三、一個誤差基準線。 之後任何 AI 模型,都要跟這個基準比。如果一個深度學習模型只比物理模型好 5%,但完全不能解釋、需要三年資料、換規格就要重訓——那它不值得。
我看過太多專案直接從「丟 XGBoost」開始,跳過這一步。結果是有一個 R² 0.92 的模型,但沒有人能回答「為什麼今天不準」。
一個誠實的記錄:我失敗過的地方
第一次做這件事,我把所有資料(含過渡段)丟進去擬合,得到的 k_screw 比合理值小了將近 20%。
原因是開機升速段佔了相當比例的資料,那段的質量守恆根本不成立(螺桿裡的料還沒穩定),模型為了同時遷就穩態和過渡態,把係數往中間拉。
我沒發現,因為 RMSE 看起來還可以。是後來拿這個係數去算另一個規格,算出來的外徑差了 0.3mm,才回頭查。
教訓是:RMSE 好不代表模型對。 一定要用物理常識檢查係數本身——k_screw 對應的產出量換算成 kg/hr,跟押出機的銘牌規格對得起來嗎?shrink_rate 2% 對 PE 合理嗎?
係數必須自己說得通,不能只靠指標。
iThome鐵人賽
看影片追技術