groupingCOSMOMargulesall.ipynb 詳細分析報告

1. 總覽

本 Jupyter Notebook 實作了一套物理約束神經網路 (PINN) 預測管線，用於預測氣液平衡 (VLE) 數據。方法結合了 COSMO 分子描述符與 Margules 活度係數模型。Notebook 載入預先訓練好的模型，針對二元混合物的熱力學性質（氣相組成 y 及系統壓力 P）進行預測，共約 17,600 筆數據。

專業背景

• COSMO (COnductor-like Screening MOdel)：一種量子化學方法，產生分子的 sigma-profile（電荷密度分布），在此作為 51 維的分子描述符使用。
• Margules 方程式：用於二元混合物活度係數的熱力學模型，在神經網路中作為物理約束層使用。
• VLE 預測：給定兩個化合物 (A 和 B)、溫度 (T) 及液相組成 (xA)，模型預測氣相組成 (yA) 和系統壓力 (P)。

---

2. Notebook 結構（15 個程式碼儲存格）

儲存格 1 — 初始設定與輔助函數

用途：匯入核心套件、定義檔案路徑、建立三個輔助函數。

關鍵函數：

• training_data(number)：根據化合物 ID 從 CSV 檔讀取 51 點的 sigma-profile。若檔案不存在則回傳零向量。
• name_to_cosmo(COSMO_num)：在對應 CSV (CID ↔ SMILES) 中查找化合物 ID。本質上是一個驗證化合物是否存在的恆等函數。
• train_datasets(fluids, T, x)：組裝 110 維特徵向量，串接兩個 sigma-profile (51×2=102)、溫度 (1)、莫耳分率 xA (1)，以及 6 個佔位零值（P、VpA、VpB、gammaA、gammaB 及備用欄位）。

特徵向量配置（110 維）：

索引	內容
0–50	化合物 A 的 Sigma-profile
51–101	化合物 B 的 Sigma-profile
102	溫度 T (K)
103	莫耳分率 xA
104	yA（目標值，佔位符）
105	ln(P)（壓力）
106	ln(VpA)（A 的蒸氣壓）
107	ln(VpB)（B 的蒸氣壓）
108	gammaA（A 的活度係數）
109	gammaB（B 的活度係數）

儲存格 2 — 已註解的多折預測（MLP V9）

用途：早期完全註解的嘗試，迭代多個交叉驗證折疊（fold），載入各模型並執行預測。使用較簡單的 MLP（非 PINN），採用雙輸入結構 [in_feat, in_phys]。

狀態：完全註解，未使用。

儲存格 3 — 單折 MLP V9 預測

用途：使用簡單 MLP 模型（V9 版本），逐列進行預測，不含 Margules 物理約束層。逐列寫入 CSV。使用 MinMaxScaler 進行特徵正規化，並將預測結果反轉回物理單位。

關鍵邏輯：

• 壓力單位偵測：P > 500 假定為 Pa 單位，否則為 kPa。
• 反正規化：氣相組成使用 y_calc = pred * (max - min) + min；壓力使用 P = exp(pred * (max - min) + min)。

儲存格 4 — 已註解的 Pathlib 版本

用途：另一個已註解的版本，嘗試使用 pathlib.Path 做更乾淨的路徑處理，並對蒸氣壓做更明確的對數轉換。未執行。

儲存格 5 — PINN 模型的匯入區塊

用途：匯入 PINN 架構所需的 TensorFlow/Keras 元件：Model、Input、Dense、Concatenate、Lambda 及 Keras 後端 K。

儲存格 6 — PINN 路徑設定（V8_7）

用途：設定 PINN 模型版本（V8，變體 7）的所有檔案路徑，包括：

• 權重檔案 (.weights.h5)
• 超參數 pickle 檔案 (.pkl)
• 預處理管線 (.joblib)
• 輸入/輸出 CSV 檔案

儲存格 7 — 核心 PINN 架構 + 首次預測嘗試

用途：定義 PINN 模型的核心元件：

1. **hard_constraint_layer_v17(inputs, **kwargs)** — Margules 物理約束層： - 接收 MLP 原始輸出 [A12, A21, y_mlp, P_delta] 及物理輸入 [xA, lnPsatA, lnPsatB, T] - 將物理輸入從 MinMax 縮放空間反正規化 - 計算 Margules 活度係數：\ln \gamma_A = x_B^2 (A_{12} + 2(A_{21} - A_{12}) x_A) - 透過修正 Raoult 定律計算總壓：P = x_A \gamma_A P_{sat,A} + x_B \gamma_B P_{sat,B} - 回傳 [y_mlp, P_normalized + P_delta, gamma_A, gamma_B]
2. **build_and_load_model(params, fold_idx, scaler)** — 建構較簡單的模型，各層寬度一致（皆為 n_neurons），使用 tanh 啟動函數。
3. 首次預測迴圈 — 載入 fold 1 的模型，執行預測，輸出結果。

儲存格 8 — 基於 PKL 的模型載入（含備用路徑）

用途：從 pickle 檔載入超參數，處理鍵名不匹配（n_neurons vs neurons vs units），並以備用路徑嘗試載入權重。

儲存格 9 — 完整 PINN 管線（V8_7，Fold 1）

用途：自成一體的版本，包含：

• build_pinn_model_from_params() 使用逐層神經元數量 (u0, u1, ...) 及 swish 啟動函數
• 輸出分支：tanh(A12, A21) * 3.0、clip(yA, 0, 1)、P_delta
• 用於正規化驗證的診斷列印語句
• 逐列預測迴圈

儲存格 10 — 幾近相同的變體

用途：與儲存格 9 相同，僅有細微差異：

• 不同的輸出檔名 (_Fix.csv)
• 使用 np.asarray() 修復潛在的 InvalidIndexError
• 無診斷列印

儲存格 11 — 同步預測與即時輸出

用途：相同的預測邏輯，變更如下：

• 版本/折疊改為 fold_idx = 5
• 即時主控台輸出格式化表格
• 累計誤差追蹤 (total_ad_y, total_ard_p)
• 最終摘要統計（AAD-y、AARD-P）

注意：p_args 中有一個 Bug — 'max_gammaB': sc.data_min_[-1] 誤用了 data_min 而非 data_max。

儲存格 12 — 最終執行的預測（17,656 列）

用途：唯一實際執行完成的儲存格，產生主要結果檔案。使用 fold 5，gamma 佔位值設為 1.0（與其他使用 0 的儲存格不同）。此儲存格產生了 17,656 筆預測。

最終結果：AAD-y = 15.98%，AARD-P = 159.43%

儲存格 13 — 批次預測（向量化）

用途：嘗試使用 model.predict() 一次對整個資料集進行批次預測 (batch_size=256)，而非逐列預測。達到相同的 AARD-P 159.43%。

儲存格 14 — Scaler 邊界檢查

用途：重新載入 scaler 並列印邊界值以驗證特徵對齊。使用 grouping_alldata_V3.csv（與其他儲存格的版本不同）。

儲存格 15 — Gamma 訓練記憶檢查

用途：列印 scaler 中 gamma 欄位（索引 108、109）的最小/最大值，揭示訓練資料範圍：gammaA ∈ [0.0001, 45354.2]，gammaB ∈ [0.025, 1083.9]。

---

3. 資料流程圖

輸入 CSV (grouping_alldata_V2.csv)
    │
    ├── comp.A, comp.B → Sigma-profile CSV 檔（各 51 維）
    ├── T, xA → 物理條件
    ├── P_exp, VpA, VpB → 實驗值
    │
    ▼
特徵向量組裝（110 維）
    │  [sigma_A(51) | sigma_B(51) | T | xA | yA | lnP | lnVpA | lnVpB | gammaA | gammaB]
    │
    ▼
MinMaxScaler（來自訓練管線）
    │
    ├── X_f = scaled[:, :104]              → MLP 輸入（sigma-profiles + T + xA）
    └── X_p = scaled[:, [103,106,107,102]] → 物理輸入 [xA, lnVpA, lnVpB, T]
         │
         ▼
    PINN 模型
    ┌──────────────────────────────────────────┐
    │  MLP 骨幹網路（swish，可變層數）           │
    │       ↓                                  │
    │  原始輸出：[A12, A21, yA, P_delta]        │
    │       ↓                                  │
    │  hard_constraint_layer_v17               │
    │  （Margules 方程式 + Raoult 定律）        │
    │       ↓                                  │
    │  輸出：[y_pred, P_pred, γA, γB]           │
    └──────────────────────────────────────────┘
         │
         ▼
    反正規化
         │
         ├── y_cal = pred[0] * (max_y - min_y) + min_y
         └── P_cal = exp(pred[1] * (max_p - min_p) + min_p)
         │
         ▼
    誤差指標（AARD-P、AAD-y）→ 輸出 CSV

---

4. 關鍵問題與程式碼異味

4.1 大量程式碼重複

最嚴重的問題。相同的邏輯在 10 個以上的儲存格中被複製貼上，僅有微小差異：

• hard_constraint_layer_v17 被定義了 6 次（儲存格 7、9、10、11、12、13）
• build_pinn_model_from_params 被定義了 5 次
• get_desc / training_data 被定義了 4 次
• 路徑設定在 8 個儲存格中重複
• 預測迴圈被撰寫了 7 次

這使得 Notebook 極難維護與除錯。修復單一 Bug 必須在多處同時修改。

4.2 硬編碼的 Windows 絕對路徑

所有路徑使用 r"C:\Users\Hung\Desktop\2025_COSMO"。這使得 Notebook：

• 無法在其他電腦上直接執行
• 無法在 macOS/Linux 上執行
• 對目錄結構變更極為脆弱

4.3 隨處可見的魔術數字

欄位索引如 103、105、106、107、[-6]、[-7] 在程式碼中頻繁出現，卻沒有命名常數。110 維特徵向量的配置未被文件化，完全依賴隱含的索引慣例。

4.4 脆弱的壓力單位偵測

使用啟發式規則 P > 500 → Pa，否則 kPa 來自動偵測壓力單位。對於壓力介於 500–1000 kPa 的系統（物理上完全可能），此規則會失效並造成無聲的資料損毀。

4.5 逐列預測迴圈

大多數儲存格在 Python for 迴圈中對 17,600+ 列的每一列呼叫 model.predict()。與批次預測相比，這種方式極度緩慢。儲存格 13 嘗試了批次預測，但僅作為事後補救。

4.6 儲存格 10/11 中的 Bug — gammaB 正規化

'max_gammaB': sc.data_min_[-1]  # Bug：應為 sc.data_max_[-1]

這會導致 gamma_B 的反正規化不正確。

4.7 Gamma 佔位值不一致

• 儲存格 1–11 使用 0 作為特徵向量中 gammaA、gammaB 的佔位值
• 儲存格 12 使用 1.0 作為佔位值
• 正確值取決於 scaler 的預期輸入，這種不一致會導致不同的預測結果

4.8 迴圈中重複讀取 CSV 檔案

training_data() 和 name_to_cosmo() 在每次呼叫時都開啟並讀取 CSV 檔案。對於 17,600 次迭代，這意味著僅 sigma-profile 就有約 35,200 次檔案開啟/關閉操作。

4.9 無進度條或時間預估

使用 print(f"已完成 {i+1}...") 每 500 列輸出一次。使用正規的進度條（如 tqdm）會更具資訊性。

4.10 未清理的註解程式碼

儲存格 2 和 4 是完全被註解的程式碼區塊（數百行），留在 Notebook 中作為「歷史參考」，但增加了雜訊和混亂。

---

5. 改善建議

5.1 重構為乾淨的 Python 模組

將所有可重複使用的邏輯抽取到一個 Python 模組中（例如 cosmo_vle_predictor.py）：

# cosmo_vle_predictor.py

import os
import numpy as np
import pandas as pd
import joblib
import csv
import pickle
import tensorflow as tf
from tensorflow.keras.models import Model
from tensorflow.keras.layers import Input, Dense, Lambda, Concatenate
import tensorflow.keras.backend as K
from pathlib import Path
from dataclasses import dataclass

# --- 常數定義 ---
SIGMA_PROFILE_DIM = 51
FEATURE_DIM = 110
MLP_INPUT_DIM = 104
PHYS_INPUT_DIM = 4

# 特徵向量欄位索引
IDX_T = 102
IDX_XA = 103
IDX_YA = 104
IDX_LNP = 105
IDX_LNVPA = 106
IDX_LNVPB = 107
IDX_GAMMA_A = 108
IDX_GAMMA_B = 109


@dataclass
class ModelConfig:
    root: Path
    results_version: int
    version: int
    fold_idx: int

    @property
    def model_train_dir(self) -> Path:
        return self.root / "model_train" / f"results_15_V{self.results_version}_{self.version}_Validation"

    @property
    def pipeline_path(self) -> Path:
        return self.model_train_dir / "model" / "pipelineMLP_COSMO_merged.joblib"

    @property
    def weight_path(self) -> Path:
        return self.model_train_dir / "model" / "h5" / f"PINN_V{self.results_version}_Fold_{self.fold_idx}.weights.h5"

    @property
    def pkl_path(self) -> Path:
        return self.model_train_dir / "model" / "pkl" / f"PINN_V{self.results_version}_Fold_{self.fold_idx}.pkl"

5.2 將 Sigma-Profile 快取於記憶體中

在啟動時一次性載入所有 sigma-profile，而非每次迭代都重新讀取檔案：

class SigmaProfileCache:
    def __init__(self, folder_path: Path):
        self.folder_path = folder_path
        self._cache: dict[str, list[float]] = {}
        self._zero_profile = np.zeros(SIGMA_PROFILE_DIM).tolist()

    def get(self, compound_id) -> list[float]:
        key = str(compound_id)
        if key not in self._cache:
            file_path = self.folder_path / f"{key}-opt-b3lyp.csv"
            if file_path.exists():
                with open(file_path, "r", encoding="utf-8") as f:
                    self._cache[key] = [float(row[1]) for row in csv.reader(f)][:SIGMA_PROFILE_DIM]
            else:
                self._cache[key] = self._zero_profile
        return self._cache[key]

5.3 使用批次預測取代逐列預測

先組裝所有特徵向量，再一次性呼叫預測：

def predict_batch(model, scaler, feature_matrix: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    scaled = np.asarray(scaler.transform(feature_matrix))
    X_f = scaled[:, :MLP_INPUT_DIM]
    X_p = scaled[:, [IDX_XA, IDX_LNVPA, IDX_LNVPB, IDX_T]]

    preds = model.predict([X_f, X_p], batch_size=256, verbose=0)

    y_cal = preds[:, 0] * (scaler.data_max_[-6] - scaler.data_min_[-6]) + scaler.data_min_[-6]
    ln_p_cal = preds[:, 1] * (scaler.data_max_[-5] - scaler.data_min_[-5]) + scaler.data_min_[-5]
    p_cal = np.exp(ln_p_cal)

    return y_cal, p_cal

5.4 使用設定檔取代硬編碼路徑

將硬編碼路徑替換為 YAML/JSON 設定檔：

# config.yaml
root: ./2025_COSMO
results_version: 8
version: 7
fold_idx: 5
input_file: grouping/data/grouping_alldata_V2.csv
sigma_profile_dir: grouping/data/s-profile-area

5.5 將物理層定義為正式的 Keras Layer

將 Lambda 替換為自訂 tf.keras.layers.Layer，以獲得更好的序列化能力和程式碼清晰度：

class MarguleConstraintLayer(tf.keras.layers.Layer):
    def __init__(self, scaler_params: dict, **kwargs):
        super().__init__(**kwargs)
        self.p = scaler_params

    def call(self, inputs):
        mlp_out, phys_in = inputs
        A12, A21 = mlp_out[:, 0:1], mlp_out[:, 1:2]
        y_mlp, P_delta = mlp_out[:, 2:3], mlp_out[:, 3:4]

        xA = phys_in[:, 0:1] * (self.p["max_xA"] - self.p["min_xA"]) + self.p["min_xA"]
        xB = 1.0 - xA
        lnPsatA = phys_in[:, 1:2] * (self.p["max_VpA"] - self.p["min_VpA"]) + self.p["min_VpA"]
        lnPsatB = phys_in[:, 2:3] * (self.p["max_VpB"] - self.p["min_VpB"]) + self.p["min_VpB"]

        ln_gammaA = tf.square(xB) * (A12 + 2.0 * (A21 - A12) * xA)
        ln_gammaB = tf.square(xA) * (A21 + 2.0 * (A12 - A21) * xB)

        P_phys = xA * tf.exp(ln_gammaA + lnPsatA) + xB * tf.exp(ln_gammaB + lnPsatB)
        P_log = tf.math.log(tf.maximum(P_phys, 1e-7))
        P_norm = (P_log - self.p["min_p"]) / (self.p["max_p"] - self.p["min_p"] + 1e-8)

        return tf.concat([y_mlp, P_norm + P_delta, tf.exp(ln_gammaA), tf.exp(ln_gammaB)], axis=-1)

    def get_config(self):
        config = super().get_config()
        config["scaler_params"] = self.p
        return config

5.6 明確的壓力單位處理

將脆弱的 P > 500 啟發式規則替換為明確的元資料處理：

def convert_pressure(value: float, unit: str) -> float:
    """將壓力轉換為 kPa，需明確指定單位。"""
    conversion = {"Pa": 1e-3, "kPa": 1.0, "bar": 100.0, "atm": 101.325, "mmHg": 0.133322}
    if unit not in conversion:
        raise ValueError(f"未知的壓力單位：{unit}")
    return value * conversion[unit]

或者，在輸入 CSV 中增加一個 P_unit 欄位。

5.7 使用 tqdm 加入進度條

from tqdm import tqdm

for i in tqdm(range(len(df_input)), desc="預測中"):
    ...

5.8 使用型別提示與文件字串

所有函數都應具備清晰的型別提示和文件字串，說明：

• 函數的功能
• 預期的輸入/輸出形狀和單位
• 參數的物理意義

5.9 移除無用程式碼

刪除所有已註解的儲存格（儲存格 2、4），將預測邏輯整合為單一、經過測試的儲存格或函數。

5.10 加入驗證與健全性檢查

def validate_feature_vector(vec: np.ndarray, scaler) -> None:
    """檢查特徵值是否在 scaler 預期的邊界範圍內。"""
    scaled = scaler.transform(vec.reshape(1, -1))
    out_of_range = (scaled < -0.5) | (scaled > 1.5)
    if out_of_range.any():
        bad_cols = np.where(out_of_range[0])[0]
        raise ValueError(f"特徵值超出 scaler 範圍的索引：{bad_cols}")

---

6. 建議的重構 Notebook 結構

乾淨版本的 Notebook 應該只有 5–6 個儲存格：

儲存格	用途
1	設定與匯入（從 YAML 載入設定、定義常數）
2	模型架構定義（單一 hard_constraint_layer、單一 build_model）
3	資料載入與特徵組裝（批次處理，含快取機制）
4	批次預測與反正規化
5	結果輸出（CSV + 摘要統計）
6	（選用）視覺化 / 誤差分析

更好的做法：轉換為 Python 腳本

對於正式的預測管線，建議將 Notebook 轉換為命令列腳本：

python predict_vle.py --config config.yaml --input data.csv --output results.csv

這提供了：

• 版本控制的友善性（避免 JSON Notebook 的差異比較困難）
• 透過設定檔實現可重現性
• 易於整合至自動化工作流程

---

7. 效能考量

面向	現狀	建議
預測方式	逐列預測（17,600 次呼叫 model.predict）	批次預測（1 次呼叫，batch_size=256）
檔案 I/O	每次迭代讀取 sigma-profile	啟動時快取於記憶體中
CSV 寫入	逐列追加 mode='a'	收集結果至串列，最後一次性寫入
模型載入	每個儲存格都從頭重建	載入一次，重複使用
預期加速	—	僅預測部分即可達 50–100 倍

---

8. 已發現的 Bug

#	位置	Bug 描述	影響
1	儲存格 10/11 p_args	'max_gammaB': sc.data_min_[-1] 應為 sc.data_max_[-1]	gamma_B 反正規化不正確
2	儲存格 12	Gamma 佔位值設為 1.0，而其他儲存格使用 0	特徵縮放不一致
3	儲存格 1 name_to_cosmo	硬編碼範圍 range(0, 1960) — 對超過 1960 個化合物的資料集會無聲失敗	遺漏化合物 ID
4	儲存格 3	縮放前未對 P_exp 做對數轉換 — 與後續儲存格的做法不同	特徵與訓練時不對齊
5	儲存格 12 p_args	缺少 'min_gammaA'、'max_gammaA'、'min_gammaB'、'max_gammaB' 鍵值	hard_constraint_layer_v17 將使用 kwargs.get() 回傳 None，導致執行時錯誤或無聲的不正確輸出

---

9. 總結

本 Notebook 是一個可運作但高度實驗性的研究程式碼，經歷了多次反覆試錯的演進。它成功實作了基於 PINN 的 VLE 預測系統，結合了 COSMO 分子描述符與 Margules 活度係數物理約束。然而，它存在嚴重的程式碼重複問題（約 10 個幾乎相同的核心邏輯副本）、硬編碼路徑、魔術數字，以及物理單位和佔位值的不一致處理。

最具影響力的改善方向為：

1. 整合所有定義為單一模組（消除約 90% 的重複程式碼）
2. 使用批次預測取代逐列預測（約 50-100 倍加速）
3. 將 sigma-profile 快取於記憶體中（約 2 倍 I/O 加速）
4. 修復 gamma 正規化 Bug（物理約束層中的錯誤）
5. 將硬編碼路徑替換為相對路徑或設定檔

預測結果（AAD-y = 15.98%，AARD-P = 159.43%）確實存在嚴重問題，尤其是壓力預測。讓我從幾個層面分析可能的原因和建議。

---

** 預測結果不好的原因！！

一、從輸出結果觀察到的模式

從儲存格 12 的 17,656 筆輸出中，我觀察到幾個明顯的問題模式：

1. 某些系統的壓力預測偏差極大

例如系統 3.0-501.0，AARD-P 高達 1600%，模型預測的壓力比實驗值高出 10-20 倍。而系統 7.0-1105 的 y_cal 經常卡在 0.0001，與 y_exp 完全不符。

2. y_cal 經常飽和在邊界值

大量預測結果的 y_cal = 0.9999 或 0.0001，這代表模型的輸出被 clip_by_value(v[:, 2:3], 0.0, 1.0) 截斷了——模型的原始輸出已超出合理範圍。

3. 某些化合物的預測全面崩潰

尤其是涉及化合物 ID 1921、1911 等的系統，幾乎所有預測都極度偏差。

---

二、根本原因分析

原因 1：Sigma-Profile 缺失（零向量問題）

這可能是最嚴重的問題。程式碼中有這樣的邏輯：

if not os.path.exists(file_full_path):

return [np.zeros(51).tolist()] # 回傳全零向量

當某個化合物的 sigma-profile 檔案不存在時，模型接收到的是 51 個零值。這等於告訴模型「這個分子沒有電荷分布」，在物理上是完全沒有意義的。模型從未在訓練中看過這種輸入，自然會產生垃圾輸出。

建議：

• 統計有多少化合物缺少 sigma-profile 檔案，將這些資料點從預測集中排除
• 或在結果中標記這些資料點，不計入整體誤差統計
• 長期方案：補算缺失化合物的 COSMO 計算

原因 2：Margules 模型的物理限制

Margules 方程式是最簡單的活度係數模型，只有兩個參數（A12、A21）。它只適用於中等偏離理想的系統。對於以下類型的系統，Margules 模型根本不足以描述：

• 強烈非理想系統（如水-有機溶劑）
• 存在液-液分相的系統
• 具有氫鍵或強極性交互作用的系統

AARD-P 159% 的結果暗示，預測資料集中可能包含大量 Margules 模型無法處理的系統。

建議：

• 考慮升級為 Wilson、NRTL 或 UNIQUAC 模型作為物理約束層，這些模型對非理想系統有更好的描述能力
• 分析哪些「類型」的系統誤差最大（例如：按化合物極性分類、按溫度範圍分類），找出模型的弱點

原因 3：訓練與預測的特徵不一致

我在程式碼中發現了幾個訓練/預測不一致的地方：

a) Gamma 佔位值不一致

訓練時 scaler 看到的 gamma 範圍是 gammaA ∈ [0.0001, 45354]，但預測時：

• 部分儲存格使用 gammaA = gammaB = 0
• 儲存格 12 使用 gammaA = gammaB = 1.0

雖然這些索引 (108, 109) 不直接輸入模型，但如果 scaler 不是逐欄獨立的 MinMaxScaler，就可能影響其他特徵的縮放。

b) yA 佔位值

索引 104 (yA) 在預測時設為 0，但訓練時有真實值。同樣的考量適用。

建議：

• 確認 scaler 是逐欄獨立的 MinMaxScaler（如果是，則此問題影響較小）
• 若有影響，使用訓練集的中位數作為佔位值，而非 0 或 1

原因 4：P_delta 的約束不足

模型輸出中的 P_delta 是一個無約束的修正項，用來補償 Margules 物理模型的誤差：

P_norm_final = P_norm_phys + P_delta

如果 P_delta 沒有正則化，模型可能在訓練時過度依賴 P_delta 來擬合訓練數據，導致物理參數（A12、A21）學得不好。在新數據上，P_delta 的泛化能力差，物理參數又不準確，預測就會崩潰。

建議：

• 在訓練損失函數中加入 P_delta 的 L2 正則化，強迫模型優先使用物理模型
• 例如：loss = MSE_y + MSE_P + lambda * ||P_delta||^2
• lambda 的值需要調整，過大會限制模型能力，過小則無效

原因 5：訓練集與預測集的分布差異

預測集 grouping_alldata_V2.csv 可能包含訓練集未涵蓋的：

• 化合物對（新的二元系統）
• 溫度範圍（外推到更高或更低溫度）
• 組成範圍（接近純組分的極端 xA）

建議：

• 比較訓練集與預測集的化合物 ID 重疊率
• 檢查溫度和壓力的分布是否一致
• 對「外推」的數據點單獨評估誤差

---

三、具體行動建議（按優先順序）

優先級	行動	預期效果
1	統計並排除缺少 sigma-profile 的化合物，重新計算誤差	可能大幅降低 AARD-P
2	分析誤差最大的前 20 個系統，找出共同特徵	定位問題根源
3	對 P_delta 加入正則化並重新訓練	改善泛化能力
4	將 Margules 升級為 NRTL 或 Wilson 模型	處理非理想系統
5	加入「可信度評估」：若 scaled 特徵超出 [0, 1] 範圍，標記為不可靠預測	避免錯誤結論
6	使用 K-fold 交叉驗證的平均結果，而非單一 fold	更穩健的評估

---

四、快速診斷腳本建議

如果你想快速定位問題，可以在 Notebook 中加入以下診斷邏輯（思路如下）：

a) 統計缺失的 sigma-profile：

• 遍歷所有 comp.A 和 comp.B，檢查 {id}-opt-b3lyp.csv 是否存在
• 統計有多少筆資料受影響

b) 分組誤差分析：

• 按二元系統 (comp.A, comp.B) 分組計算平均 AARD-P
• 找出誤差最大的前 20 個系統
• 檢查這些系統是否有缺失的 sigma-profile

c) 特徵範圍檢查：

• 對預測集的每個特徵，檢查 scaled 值是否在 [0, 1] 範圍內
• 超出範圍的特徵代表模型在外推

這些診斷結果會幫助你判斷是資料問題（缺失 sigma-profile、分布差異）還是模型問題（Margules 限制、過擬合）。如果需要我協助實作這些診斷腳本，可以切換到 Agent 模式。

診斷腳本位於 grouping/code/diagnostic_prediction.py，輸出檔案在 grouping/code/diagnostic_output/ 資料夾中。gamma 欄位的負值問題是最值得優先調查的方向——需要確認這些值是 gamma 還是 ln(gamma)，因為這會直接影響物理約束層的計算。