波峰疊在一起，怎麼辦？

1.一個你真的會遇到的問題入門

先講一個情境。理解了這個情境，後面所有的數學才有意義。

你在食品分析實驗室工作，桌上有一瓶市售的機能性飲料。品管主管要你回答一件事： 這瓶飲料裡的咖啡因有沒有超過標示？防腐劑（苯甲酸、己二烯酸）有沒有超過法規上限？ （己二烯酸即 sorbic acid，我國《食品添加物使用範圍及限量暨規格標準》採此名稱，部分教科書譯作「山梨酸」。）

你把樣品打進 HPLC（高效液相層析儀，High-Performance Liquid Chromatography）， 機器跑了十分鐘，吐出一張圖。你期待看到三個漂亮、分得開的波峰， 每個波峰代表一種成分，峰的面積就是它的含量。

結果你看到的是這樣：

在「沒調好」的那張圖裡，咖啡因和苯甲酸的波峰疊在一起。 這不是小事——你沒辦法量出任何一個的正確面積，也就算不出正確的含量。 在檢驗報告上，這代表數據不能用，整批樣品要重跑。

2.HPLC 到底在做什麼入門

如果你把 HPLC 想成一場賽跑，幾乎所有名詞都會瞬間變得直觀。

三個角色

比賽怎麼進行

每一個分子在管柱裡不斷重複兩件事：被黏在顆粒表面上（不動）， 和溶在液體裡被水流帶著走（前進）。

管柱的另一端有一個偵測器（detector），通常是紫外光偵測器。 它像終點線的攝影機，每一瞬間都記錄「現在流過去的液體裡有多少東西」。 把這個訊號對時間畫成一條線，就是層析圖（chromatogram）—— 本教材從頭到尾都在看的那張圖。

3.讀懂一張層析圖入門

一張層析圖只有兩個軸，但藏了三種資訊。看懂這三種，你就看懂了八成。

真實的層析圖長什麼樣實測資料

上面那張是畫出來的示意圖。真實儀器吐出來的，其實沒那麼乾淨。 下面三張是本校實驗室用 Shimadzu HPLC 實際跑出來的資料 （UV 230 nm、10 分鐘、每 0.5 秒取一點，共 1201 個點）， 分析的正是苯甲酸與己二烯酸——第 1 節那個飲料情境的主角。

三件示意圖沒告訴你的事

4.三個一定要記住的數字入門

整篇論文的公式都建立在這三個量上面。它們其實都只是「時間」的不同講法。

t₀ — 死時間（column dead time）

一個完全不黏管柱的分子，從打進去到流出來需要的時間。 它等於管柱裡液體的體積除以流速，代表「不可能比這更快」的下限。 一般分析管柱大約是 1 分鐘左右。在論文所有函數裡它叫 t0，是你必須自己量、自己輸入的參數。

t_R — 滯留時間（retention time）

某個成分的波峰最高點出現的時間。這是你在層析圖上直接讀到的數字。 論文中寫作 tR。

k — 滯留因子（retention factor）

這是最關鍵、也最常被搞混的一個。它的定義是：

5.解析度 R_s：這份教材的主角入門

前面說「波峰要分得開」，但「分得開」不能只憑感覺——需要一個數字。這個數字就是解析度。

解析度（resolution，記作 R_s）衡量相鄰兩個波峰分離得多好。定義是：

這個式子在講一件很單純的事：「兩個峰離多遠」除以「兩個峰本身有多胖」。 分子越大（離得遠）或分母越小（峰越瘦），解析度就越好。

要記住的三個門檻

用真實資料驗算一次 R_s實測資料

公式講完了，來驗算。第 3 節那三張真實層析圖，儀器軟體自己就會報告 R_s。 我們不看它的答案，直接從原始訊號（1201 個數字）自己量峰寬、自己算一次，再對照：

6.那顆能改變一切的旋鈕：有機修飾劑中階

你已經知道要讓 R_s 變大。問題是——你到底能轉哪個旋鈕？

在逆相層析（reversed-phase，就是用 C18 這種油性管柱的主流做法）裡， 最有力、最常用、也最容易調的一個變數，是移動相裡有機溶劑的比例。

移動相是水和有機溶劑（乙腈 acetonitrile 或甲醇 methanol）的混合液。 這個有機溶劑就叫有機修飾劑（organic modifier）， 它佔的體積分率記作 φ（希臘字母 phi）。 在這份教材和論文的所有 R 程式碼裡，它一律寫成 f—— 請記住 f 就是 φ，論文特別註明了這件事。

7.模型 1：把直覺變成一條公式中階

「φ 越大跑越快」是定性的。要讓電腦幫你算，需要一條定量的式子。

真實的實驗數據告訴我們一件很漂亮的事：如果你把 ln k（滯留因子取自然對數） 對 φ 畫圖，資料點幾乎會落在一條直線上。

這就是論文 Table 2 的模型 1，在 R 程式裡寫作 model=1。 它是整個套件裡最簡單的一個，也是我們整份教材的主角。 （另外還有九個模型，處理 pH、雙溶劑等更複雜的情況——放在 附錄，你現在不需要它們。）

用真實數據看看

下面是論文附帶的真實實驗數據（Data.xlsx 的 i-ret.fit 工作表）。 9 個成分，在 4 種不同的 φ 下各跑一次，管柱死時間 t₀ = 1.4 分鐘。 這批樣品是酚類污染物與芳香烴——像五氯酚、溴酚這些， 正是飲用水和食品安全檢驗會遇到的目標物。

把每個數字換成 k = (t_R − 1.4)/1.4，再取 ln，然後畫圖。 下面這張圖就是結果——九條漂亮的直線：

8.親手最佳化：把旋鈕轉一遍中階

現在你有了九個成分的 c₀ 和 c₁，就能預測任何 φ 下的完整層析圖。 下面這個模擬器把整件事串起來了——請務必動手拖。

電腦是怎麼「找到」最佳解的

你剛剛用手拖滑桿，其實就是在做最佳化——只是很慢。 iopt 做的事一模一樣，只是它每隔 0.005 把 φ 從 0.30 掃到 0.60 （這個步長就是參數 df），一共 61 個點，每個點都算一次最小 R_s， 然後挑出在時間預算內 R_s 最大的那一個。

下面這張圖就是掃描的完整結果。論文的 Fig. 7 畫的就是這種圖：

那個藏在 φ = 0.35 的陷阱

如果你剛剛按過「陷阱 φ = 0.350」，應該會嚇一跳： 明明 φ 更小、時間更長、照理說應該分得更開， R_s 卻從 φ=0.325 的 1.50 暴跌到 0.235。

9.當一個 φ 值救不了你：梯度沖提中階

剛才我們發現：10 分鐘的預算下，不管 φ 怎麼調都達不到基線分離。 這不是運氣不好，而是一個有名字的結構性問題。

一般沖提問題（the general elution problem）

當你的樣品裡同時有很親水的成分和很疏水的成分時，你會陷入兩難：

兩邊都對，也都不夠。沒有任何單一的 φ 能同時滿足兩端。 這就是「一般沖提問題」。

梯度的三個參數

10.第二顆旋鈕：pH中階

φ 和梯度控制的是「所有成分一起變快或變慢」。 但有一顆旋鈕不一樣——它能改變「誰跑得比誰快」，甚至把出峰順序整個翻過來。

回到第 1 節那瓶飲料。三個目標物裡，咖啡因是中性分子， 但苯甲酸和己二烯酸都是弱酸。這個差別在調 φ 時看不出來， 一旦你動 pH，它就會主宰一切。

轉折點在哪裡：pK_a

每個弱酸都有一個專屬的數字 pK_a，代表它「剛好一半解離」時的 pH。 在 pK_a 以下它主要是中性的，以上主要是離子態，轉變集中在前後約兩個 pH 單位之間。

注意這兩個 pK_a 差了 0.56。 就是這 0.56 的差距，讓兩根酸在不同 pH 下「各自以不同的速度提早出來」， 也讓整個出峰順序有機會被改寫。

描述它的公式：模型 7

這就是論文 Table 2 的模型 7（model=7）。 它和模型 1 的角色完全對稱：模型 1 描述 φ 怎麼影響 k，模型 7 描述 pH 怎麼影響 k。 附錄裡的模型 8 和 10 則是把兩者合起來，同時處理 φ 和 pH。

三段不同的出峰順序

把 pH 從 2.5 拖到 7.0，你會經過三個截然不同的區間：

延伸案例：真實的論文怎麼做同一件事

你到目前為止做的（拖一個變數、看最小 R_s、避開陷阱、在時間與解析度之間取捨）， 就是方法開發的核心工作。這裡看一篇真實的研究，它處理的正是同一批成分。

11.架好你的 R 環境進階

從這一節開始，你會動手用真正的統計軟體 R 來跑 RChromOptim——這是 Zisi、Pappa-Louisi 與 Nikitas 三位學者在 2020 年發表於《Journal of Chromatography A》的一套工具，把前面幾節教你手算的滯留模型（例如 lnk = c0 − c1·f）、解析度 Rs、峰形，全部包成十個可以直接呼叫的函數。你不需要自己寫程式，只要照著步驟把資料餵給它。

R Console安裝套件
# 一次安裝四個套件，中間用逗號分隔
install.packages(c("optimx", "VGAM", "plot3D", "numDeriv"))

R Console
load("RChromOptim.RData")

R Console檢查
ls()

12.你的資料要長什麼樣子進階

RChromOptim 只吃一種格式的檔案：Tab 字元分隔的純文字檔（.txt），而且永遠只有一列標題（header）。這一節帶你認識四種會用到的檔案長相，之後每個範例都會對照到其中一種。

(a) 恆溶劑滯留資料（isocratic retention data）——給 ifitk、ifitopt 用

第一欄是實驗條件（通常是有機溶劑比例 f），後面每一欄是一個溶質（solute）在該條件下量到的滯留時間 tR（分鐘）。

如果你的滯留模型是 model 7（跟 pH 有關），第一欄就要換成 pH；model 8、9、10 則需要兩欄條件（例如兩種溶劑或 pH＋f 併用）。

(b) 梯度滯留資料（gradient retention data）——給 gfitk 用

前面固定是四欄條件，之後每欄一個溶質。這四欄很容易搞混，務必記住：f1、f2 是梯度起點與終點的有機溶劑比例，t1、t2 是梯度斜坡開始與結束的時間——不是兩種溶劑的意思。

(c) 層析圖原始曲線（chromatogram trace）——給 fitshape、以及 ipred/gpred 的 ic=1 用

就是一條真實的層析圖曲線，永遠只有兩欄：時間、偵測器訊號值。

(d) 溶質參數檔（solute-parameter file）——ifitk/gfitk 的輸出，也是 pred/opt 的輸入

這種檔案是「轉置」的：第一欄是參數名稱，之後每欄是一個溶質。跑完 fitshape 之後，峰形參數（h0, h1, h2, s0, s1）會被附加在同一個檔案下方。

不同模型的參數列名稱不一樣：model 1–6 用 c0, c1, c2；model 7 用 k0, k1, pK；model 8 用 c0, c1, pK, r；model 9 用 c0…c5；model 10 用 c0, c1, c2, c3, pK。

13.十個函數怎麼選進階

十個函數名字看起來很像，但只要抓住命名規則，五秒鐘就能認出每一個在做什麼：字首 i = 恆溶劑（isocratic）、g = 單線性梯度（gradient）、mg = 雙線性梯度（bilinear gradient）；字尾 fitk = 擬合滯留模型、fitshape = 擬合峰形、pred = 預測/模擬層析圖、opt = 找最佳條件。

決策流程：我該用哪一個？

先問自己兩個問題：① 你的實驗是恆溶劑還是梯度？（決定用 i 開頭還是 g/mg 開頭）② 你現在要做的是「求參數」「模擬」還是「找最佳條件」？（決定用 fitk、pred 還是 opt）。

不論恆溶劑或梯度，正規的操作永遠是同一套四步驟流程：

14.完整實戰：恆溶劑最佳化進階

這一節帶你走完整套流程，資料是真實案例：9 個溶質（酚類污染物與芳香烴，跟食品/飲用水安全檢測相關）在 Kinetex 2.6 µm XB-C18 150×4.6 mm 管柱上、乙腈/水（pH 5.7）系統分析，t0 = 1.4 分鐘。

Step 1：準備滯留資料 .txt

把第 12 節格式 (a) 的表格存成 tab 分隔的 .txt（例如 retention.txt），第一欄是 f，A1–A9 是 9 個溶質。

Step 2：擬合滯留模型（ifitk）

R Console擬合 model 1
ifitk(model = 1, t0 = 1.4)
# 執行後跳出檔案選取視窗，用滑鼠點選 retention.txt

model 1 就是你在前面章節學過的 lnk = c0 − c1·f。輸出會是一張表，欄位大致包含：

• 參數估計值（c0、c1）——每個溶質各一組。
• 標準差（standard deviation）——這個參數估得穩不穩，數字越小代表估計越可靠。
• p 值（p-value）——白話說：p 值小（通常小於 0.05）表示這個參數確實有在起作用，不是雜訊湊出來的巧合；p 值大則代表這個參數可能可有可無。
• SEE（standard error of estimate，估計標準誤）——代表模型算出來的 tR 跟實際量到的 tR 差多少，SEE 越小，代表模型跟實驗數據貼合得越好。
• ssr（sum of squared residuals，殘差平方和）——所有誤差平方加總，是計算 SEE 的原始材料，越小越好。

本例用 9 個溶質實際擬合出來的結果（model 1）：

ifitk 會另外跳出存檔視窗，讓你把這張參數表存成「格式 (d)」的溶質參數檔（例如 param.txt）——之後每一步都要用到它。

Step 3：擬合峰形（fitshape）

R Console擬合峰形
fitshape(pd = 1)
# 跳出視窗選一個或多個真實層析圖檔（格式 c）
# 圖上會提示你用滑鼠點兩下，標出峰的範圍

跑完後，h0、h1、h2（峰高相關）與 s0、s1（峰寬相關）這五列參數，會被附加到你在 Step 2 存的 param.txt 下方，讓這個檔案變成完整的六列以上參數檔。

Step 4：驗證模型（ipred）

R Console模擬並比對真實層析圖
ipred(model = 1, t0 = 1.4, fp = 0.45)
# 先選 param.txt（含峰形列），再選要比對的真實層析圖檔

fp = 0.45 代表你要模擬「f = 0.45」這個條件下的層析圖。因為 ic=1 是預設值，畫面會把「模型模擬出來的曲線」跟「你手上真實量到的曲線」疊在一起。如果兩條線對得起來，代表模型跟峰形參數都可信，可以放心進入最佳化；如果對不起來，要回頭檢查 Step 2、3 是不是哪裡出了問題。

Step 5：找最佳分離條件（iopt）

R Console最佳化
iopt(model = 1, t0 = 1.4, tmin = 2, tmax = 20,
     fmin = 0.3, fmax = 0.6, df = 0.005)
# 選 param.txt；tmax=20 代表你能接受的最長分析時間是 20 分鐘

tmax 是整個最佳化裡最關鍵的參數——它是你的「時間預算」：不管 Rs 理論上能拉多高，分析時間都不能超過這個上限。df 則是搜尋 f 值時的網格步長，越小搜得越細但越慢。

R Consoleiopt 輸出（數值為本教材以同一演算法重算）
Table for manual selection of the optimum f value
     f      Rs   tRmax
 0.300   2.824   38.35   <- Rs 很好，但遠超過 tmax
 0.305   2.555   36.94
 ......................
 0.350   0.235   26.46   <- 共沖提陷阱（見第 8 節）
 ......................
 0.410   2.455   17.13   <- 最佳解
 ......................
 0.595   0.167    5.14
 0.600   0.002    5.00   <- 全部擠在一起

Optimum f    :   0.410
Optimum Rs   :   2.455
Optimum tR values:
3.53  4.39  6.26  7.07  7.70  9.41  10.21  14.23  17.13

白話翻譯這個結果：在 f = 0.410（也就是 41% 有機溶劑）這個條件下，9 個溶質裡分得最差的那一對相鄰峰的 Rs = 2.455（遠超過安全門檻 1.5，分離良好），而整批跑完最慢的溶質要 tRmax = 17.13 分鐘，還在你設定的 20 分鐘預算之內。

時間預算 vs. 解析度：同一組資料，換一個 tmax 會怎樣

只改 tmax，其他都不變，看看結果怎麼變化：

15.一步到位：ifitopt進階

如果你只是想快速看看「這批資料大概能不能分開、大概要多久」，不想一步步跑 ifitk → fitshape → iopt 三個函數，ifitopt 把這三步合併成一次呼叫——但代價是它只支援恆溶劑、且只支援 model 1、2、3（比較簡單的滯留模型）。

R Consoleifitopt 一次做完
ifitopt(model = 1, t0 = 1.4, blc = c(1, 3), pd = 1,
        tmax = 20, fmin = 0.3, fmax = 0.6)
# 依序跳出視窗：先選滯留資料 .txt，再選 blc 指定列數對應的層析圖檔

什麼時候用 ifitopt，什麼時候乖乖走三步驟

16.現代 R 版：tidyverse + ggplot2進階

論文的套件寫於 2019 年，是 base R 風格的互動式程式。 R 語言在這幾年變化很大——本節用 R 4.6.1 + tidyverse 2.0 + ggplot2 4.0 把同一套方法重寫一次， 並附上可直接下載執行的完整腳本。

改了哪些地方

核心：用 nest + map 一次擬合所有溶質

這是整份改寫最能展現差別的地方。原程式要用迴圈跑過每個溶質、把結果塞進矩陣； tidyverse 的作法是把每個溶質的資料收進一格裡（nest），再對每一格套用同一個函式（map）。

R逐溶質擬合模型 1：ln k = c0 − c1·f
fit_model1 <- function(data, t0 = T0) {
  data |>
    mutate(lnk = log((tR - t0) / t0)) |>      # 先換成 ln k
    nest(.by = solute) |>                        # 每個溶質收成一列
    mutate(
      fit    = map(data, \(d) lm(lnk ~ f, data = d)),   # 每列各跑一次迴歸
      tidy   = map(fit, broom::tidy),
      glance = map(fit, broom::glance)
    ) |>
    mutate(
      c0 = map_dbl(tidy, \(x) x$estimate[x$term == "(Intercept)"]),
      c1 = map_dbl(tidy, \(x) -x$estimate[x$term == "f"]),  # 斜率取負才是 c1
      r_squared = map_dbl(glance, "r.squared"),
      .keep = "unused"
    )
}

解析度與最佳化：把迴圈換成資料表

lead() 讓「相鄰兩峰相減」變成一行；掃描最佳化則是先把每個 φ 的完整結果算出來收進表格， 再從表格裡挑出符合時間預算、R_s 最大的那一列。

R相鄰峰解析度 + 掃描最佳化（對應 iopt）
# Rs = 2(tR2 - tR1) / (w1 + w2)，只算相鄰對
resolution_table <- function(peaks) {
  peaks |>
    arrange(tR) |>
    mutate(pair = paste(solute, lead(solute), sep = "/"),
           Rs   = 2 * (lead(tR) - tR) / (w + lead(w))) |>
    filter(!is.na(Rs)) |>
    select(pair, Rs)
}

optimise_isocratic <- function(params, f_range = c(0.30, 0.60), df = 0.005,
                               t_max = 20, t0 = T0) {
  scan <- tibble(f = seq(f_range[1], f_range[2], by = df)) |>
    mutate(
      peaks  = map(f, \(x) predict_isocratic(params, x, t0)),
      worst  = map(peaks, min_resolution),
      Rs     = map_dbl(worst, "Rs"),
      tR_max = map_dbl(peaks, \(p) max(p$tR)),
      feasible = tR_max <= t_max          # 時間預算
    )
  best <- scan |> filter(feasible) |> slice_max(Rs, n = 1)
  list(scan = scan, best = best)
}

梯度：數值解基本梯度方程式

第 9 節那條積分式 ∫dt/(t₀k(t)) = 1， 在這裡用累積和解決——把時間切成小格，每格前進一點點，累積到 1 就出管柱。 完全向量化，沒有迴圈。

R梯度沖提的滯留時間
predict_gradient_one <- function(c0, c1, f_start, f_end, tG, tD, t0 = T0, dt = 0.002) {
  tt  <- seq(0, 400, by = dt)
  phi <- case_when(                            # 管柱入口看到的組成（含延遲 tD）
    tt <= tD              ~ f_start,
    (tt - tD) / tG >= 1   ~ f_end,
    .default              = f_start + (f_end - f_start) * (tt - tD) / tG
  )
  k    <- exp(c0 - c1 * phi)
  prog <- cumsum(dt / (t0 * k))                 # 已走完的管柱比例
  i    <- which(prog >= 1)[1]
  before <- if (i > 1) prog[i - 1] else 0
  te   <- tt[i] + dt * (1 - before) / (prog[i] - before)   # 格內線性內插
  tibble(tR = te + t0, k_elute = k[i])
}

R退化檢查
iso   <- predict_isocratic(params, 0.45)
degen <- predict_gradient(params, 0.45, 0.45, tG = 15)
stopifnot(all(abs(sort(iso$tR) - sort(degen$tR)) < 1e-4))

執行結果

終端機Rscript rchromoptim_modern.R
== 模型 1 擬合結果 ==
  solute    c0    c1 r_squared  sigma   p_c1     n
1 A1      3.14  6.64     0.998 0.0333 0.0011     4
2 A2      3.46  6.60     0.998 0.0298 0.0009     4
3 A3      3.67  5.91     0.999 0.0201 0.0005     4
  ...
9 A9      5.60  7.77     0.997 0.0433 0.0014     4

c0 範圍 3.1404 - 5.6030
c1 範圍 5.9113 - 8.5366
論文 i-optim 工作表公布值：c0 3.1404-5.6030，c1 5.9113-8.5366   <- 完全一致

== 最佳化 (t_max = 20) ==
 f = 0.410  Rs = 2.453  tR_max = 17.13  最難分 = A4/A5
 t_max = 15 -> f = 0.430  Rs = 2.027
 t_max = 10 -> f = 0.490  Rs = 0.914

== 梯度 (0.35->0.70, tG=11, tD=0.7) ==
 Rs = 1.633 (A4/A5)  tR_max = 9.99
 檢查通過：f_start == f_end 時梯度解退化為恆溶劑解（誤差 < 1e-4 min）

ggplot2 的輸出

三張圖用 patchwork 的 / 與 | 運算子拼起來，一行程式碼：

R組圖並存檔
p <- (plot_model1(data, params) | plot_optimisation(opt)) /
      plot_chromatogram(best_peaks, "最佳條件 f = 0.410")

ggsave("demo.png", p, width = 12, height = 8, dpi = 130)

還能再快嗎：從「跑得動」到「跑得動大的」

上面的腳本跑示範資料只要零點幾秒，看起來沒有優化的必要。 但真正的最佳化很少只掃一個變數——論文的 iopt 對模型 8、10 要掃 φ × pH 的二維網格， gopt 要掃 t_G × pH。條件數從幾十跳到幾千， 原本沒感覺的成本就會變成好幾分鐘的等待。

優化一：把 61 個資料框變成 1 個矩陣

原本的寫法對每一個 φ 都建一個 tibble，再從中取最小 R_s。 61 個 φ 就配置 61 次記憶體、跑 61 次排序與資料框運算。 但這件事的本質只是一個矩陣：9 個溶質 × 61 個條件。

R優化後：整張表一次算完
scan_isocratic <- function(params, f, t0 = T0, shape = SHAPE) {
  tR  <- t0 * (1 + exp(params$c0 - outer(params$c1, f)))   # n_solute x n_f，一行算完
  ord <- apply(tR, 2, order)                               # 每欄的沖提順序
  S   <- apply(tR, 2, sort)
  W   <- 4 * (shape$s0 + shape$s1 * S) / sqrt(2)
  n   <- nrow(S)
  Rs  <- 2 * (S[-1, , drop = FALSE] - S[-n, , drop = FALSE]) /   # 矩陣位移相減
             (W[-1, , drop = FALSE] + W[-n, , drop = FALSE])
  tibble(f = f, Rs = apply(Rs, 2, min), tR_max = S[n, ], ...)
}

outer(c1, f) 一次產生所有「c₁ × φ」的組合； S[-1,] - S[-n,] 就是「每一欄裡相鄰兩個峰相減」。 迴圈消失了，因為問題本來就是矩陣形狀的。

優化二：梯度的三個浪費

原本每個溶質各呼叫一次數值求解器，裡面藏了三筆重複成本：

• 梯度曲線 φ(t) 被重算 9 次，但它對所有溶質完全相同——只算一次就好。
• 每個溶質都配置 200,001 個時間格點（固定算到 400 分鐘）， 但實際上所有成分 10 分鐘內就出完了。改成由 φ_max 的恆溶劑滯留時間估一個合理上限， 不夠再放大，格點數降到約 1 萬。
• k 逐溶質算，可以合併成一個 (溶質 × 時間) 矩陣。

實測結果

注意第二列：掃描越細，優化的效益越大（42× → 152×）。 這正是重點——優化真正的價值不在於讓示範跑得快一點， 而在於讓原本做不到的規模變得做得到。 Python 版做了完全相同的兩項優化。

R新舊對照
o <- opt_old(p)$best;  nw <- optimise_isocratic(p)$best
cat(abs(o$Rs - nw$Rs))              # -> 0.00e+00
cat(max(abs(sort(grad_old(...)$tR) - sort(predict_gradient(...)$tR))))
                                       # -> 0.00e+00 min

17.Python 版：pandas + numpy + scipy進階

同一套方法，用 Python 再寫一次。這不只是「換個語言」—— 兩個獨立實作跑出同一組數字，本身就是最強的正確性證據。

同一件事，兩種語言怎麼寫

把兩份腳本並排看，是很有效的學習方式——你會發現思路完全一樣，只是詞彙不同。

擬合：groupby 取代 nest

Python逐溶質擬合模型 1
def fit_model1(data: pd.DataFrame, t0: float = T0) -> pd.DataFrame:
    """以模型 1 擬合滯留資料：ln k = c0 - c1 * f。"""
    if not (data["tR"] > t0).all():
        raise ValueError("所有 tR 必須大於 t0，否則 ln k 無定義")

    rows = []
    for solute, grp in data.groupby("solute", sort=False, observed=True):
        f   = grp["f"].to_numpy(dtype=float)
        lnk = np.log((grp["tR"].to_numpy(dtype=float) - t0) / t0)
        res = stats.linregress(f, lnk)
        rows.append(dict(
            solute=solute,
            c0=res.intercept,
            c1=-res.slope,          # 斜率取負才是 c1
            r_squared=res.rvalue**2,
            p_c1=res.pvalue,
        ))
    return pd.DataFrame(rows)

梯度：用 searchsorted 找沖提點

Python 版用 np.searchsorted() 在累積陣列裡二分搜尋，比逐格比對更快也更乾淨：

Python數值解基本梯度方程式
tt  = np.arange(0, t_limit + dt, dt)
phi = np.where(                                # 管柱入口看到的組成（含延遲 tD）
    tt <= tD,
    f_start,
    np.where((tt - tD) / tG >= 1,
             f_end,
             f_start + (f_end - f_start) * (tt - tD) / tG),
)
k    = np.exp(c0 - c1 * phi)
prog = np.cumsum(dt / (t0 * k))               # 已走完的管柱比例

idx    = np.searchsorted(prog, 1.0)            # 第一個 >= 1 的位置
before = prog[idx - 1] if idx > 0 else 0.0
te     = tt[idx] + dt * (1 - before) / (prog[idx] - before)
return te + t0, k[idx]

執行結果：和 R 版一模一樣

終端機python rchromoptim_modern.py
== 模型 1 擬合結果 ==
solute     c0     c1  r_squared  sigma   p_c1  n
    A1 3.1404 6.6413     0.9977 0.0333 0.0011  4
    A2 3.4649 6.6043     0.9981 0.0298 0.0009  4
    ...
    A9 5.6030 7.7662     0.9972 0.0433 0.0014  4

c0 範圍 3.1404 - 5.6030
c1 範圍 5.9113 - 8.5366
論文 i-optim 工作表公布值：c0 3.1404-5.6030，c1 5.9113-8.5366

== 最佳化 (t_max = 20) ==
 f = 0.410  Rs = 2.4535  tR_max = 17.13  最難分 = A4/A5
 t_max = 15 -> f = 0.430  Rs = 2.027
 t_max = 10 -> f = 0.490  Rs = 0.914

== f = 0.350 的共沖提檢查 ==
 最小 Rs = 0.2343（A7/A6）           <- 第 8 節那個陷阱

== 梯度 (0.35->0.70, tG=11, tD=0.7) ==
 Rs = 1.6330 (A4/A5)  tR_max = 9.99

檢查通過：f_start == f_end 時梯度解退化為恆溶劑解
         （最大誤差 8.67e-13 min < 1e-4）

18.互動筆記本：Quarto 與 Jupyter進階

第 16、17 節的腳本是「跑完看結果」。這一節提供兩份邊讀邊改的筆記本—— 文字說明和程式碼交錯排列，你可以隨時改一個數字、按下執行，馬上看到差別。

筆記本裡有什麼

兩份的結構刻意做成一樣，方便對照：

一個做中文教材一定會遇到的坑

R自動挑一個系統上存在的中文字型
cjk_font <- local({
  cands <- c("Microsoft JhengHei", "Noto Sans TC", "PingFang TC",
              "Heiti TC", "Noto Sans CJK TC", "SimSun")
  avail <- if (requireNamespace("systemfonts", quietly = TRUE))
    systemfonts::system_fonts()$family else character(0)
  hit <- cands[cands %in% avail]
  if (length(hit)) hit[1] else ""       # 找不到就用預設值
})
theme_set(theme_minimal(base_family = cjk_font))

19.十個滯留模型完整對照表

大部分食品分析的工作，其實只會用到模型 1 或模型 3。只有當你遇到「可解離的化合物」（例如有機酸、胺類——像乙酸、檸檬酸、生物胺）、或是需要控制 pH 值來調整分離效果時，才需要翻到這張完整清單。

20.函數速查表

RChromOptim 套件一共有 10 個函數，命名有規律可循，看懂前綴和字根就能猜出功能。

21.常見問題排除

這些是使用 RChromOptim 時真實會遇到的狀況，坦白列出來，照著檢查通常就能解決。

R安裝相依套件
install.packages(c("optimx","VGAM","plot3D","numDeriv"))

22.中英名詞對照與延伸閱讀

整理課程中出現過的關鍵術語中英對照，以及原始論文的補充資源。

延伸閱讀

本教材食品分析情境所依據的文獻

這門課在整條學習路徑上的位置

教材中的實測層析資料

本教材使用的軟體與版本

授權、版權與引用方式

你學會了什麼

• 能讀懂層析圖的基本結構，認出峰、基線與時間軸代表的意義
• 會計算滯留因子 k，並理解它與 t0、tR 的關係
• 能判斷兩個相鄰峰的 Rs 是否達到基線分離的標準
• 看得懂 ln k 對 f 的直線圖在表達什麼，以及曲率代表的意義
• 知道什麼情況該從恆溶劑沖提改用梯度沖提
• 能準備符合欄位格式要求的 R 輸入檔（滯留資料、層析圖資料）
• 能跑完 ifitk → fitshape → ipred → iopt 這條完整分析流程
• 能解讀最佳化功能輸出的最佳條件、Rs 與層析圖模擬結果
• 理解分析時間與解析度之間必須取捨，不存在「又快又完美分離」的萬能解
• 知道遇到不同的化合物或分離情境時，該去哪裡查十個滯留模型該用哪一個

23.附錄：pH 實驗指引（可印出帶進實驗室）

第 10 節的 pH 模擬器是算出來的。這份指引讓你用真實的儀器去驗證它—— 三支標準品、五種移動相、一個下午。

1. 器材與試劑

2. 標準品配製

3. 移動相配方

五種移動相，差別只有緩衝液的 pH，其餘完全相同——這是本實驗唯一該變的變數。

4. 儀器條件與操作順序

5. 資料記錄表

6. 預期結果與判讀

7. 資料分析：擬合模型 7

把每個 pH 下的 t_R 換成 k = (t_R − t₀)/t₀， 再擬合第 10 節那條式子，就能得到每個化合物的 k₀、k₁ 與 pK。

R模型 7 擬合（已實測驗證）
library(tidyverse)

# 模型 7：k = (k0 + k1*10^(pH-pK)) / (1 + 10^(pH-pK))
fit_model7 <- function(pH, tR, t0) {
  k <- (tR - t0) / t0
  m <- nls(k ~ (k0 + k1 * 10^(pH - pK)) / (1 + 10^(pH - pK)),
           start = list(k0 = max(k), k1 = min(k), pK = median(pH)),
           control = nls.control(maxiter = 200, warnOnly = TRUE))
  tibble(k0 = coef(m)[["k0"]], k1 = coef(m)[["k1"]],
         pK = coef(m)[["pK"]], rss = sum(residuals(m)^2))
}

# 換成你自己的數字
pH  <- c(3.02, 4.03, 4.47, 5.01, 5.98)   # 實測 pH，不是目標值
tR  <- c(5.37, 5.10, 4.42, 3.55, 2.71)   # 苯甲酸的滯留時間
fit_model7(pH, tR, t0 = 1.20)

8. 報告該回答的問題

• 三個預測分別成立嗎？用你的數據說明，不要只寫「符合預期」。
• 你量到的交叉 pH 是多少？和第 10 節模擬器的 4.45 / 5.42 差多少？為什麼會差？
• 擬合出的 pK 與文獻值差多少？這個差距是誤差還是有系統性的原因？
• 如果要你為這三個成分挑一個 pH 寫進 SOP，你選哪個？理由是什麼 （回顧第 10 節思考題三：耐用性 vs 速度）。
• 咖啡因在這個實驗裡扮演什麼角色？如果拿掉它，實驗的說服力會少掉什麼？