3. 一级反应速率常数的最小二乘估计

为什么一级反应要同时估计 c₀ 和 k

间歇反应器中一级动力学满足 c(t) = c₀·exp(−k·t)。若把 ln c 对 t 做线性回归，斜率给出 −k，截距给出 ln c₀，看似方便，但存在两个问题：一是对数变换会放大低浓度点的误差权重，二是 c₀ 往往并非精确已知（取样时刻偏差、初始混合不均、浓度分析系统误差都会影响它）。因此更稳妥的做法是直接在浓度域做非线性最小二乘，把 c₀ 与 k 作为两个待估参数一起拟合，再检查残差是否呈随机分布。

本案例给出一组教学示例数据（假设值，非实测），用 Gauss–Newton 迭代完成双参数估计，并给出残差随机性判读方法。

示例数据与模型设定

设初始浓度名义值约 1.00 mol/L，取样时间 0–60 min，共 9 个点。数据中人为加入少量随机扰动，用于演示拟合与残差检查流程。

目标函数为残差平方和：S(c₀,k) = Σ [cᵢ − c₀·exp(−k·tᵢ)]²。对 c₀ 与 k 分别求偏导，得到 Jacobian 元素：

∂c/∂c₀ = exp(−k·t)
∂c/∂k  = −c₀·t·exp(−k·t)

Gauss–Newton 更新式为 Δ = (JᵀJ)⁻¹ Jᵀ r，其中 r 为当前残差向量，参数按 θ ← θ + Δ 迭代至 Δ 足够小。

可运行的拟合实现

下面代码只用标准库完成 2×2 正规方程求解，可直接运行并输出参数、标准误与残差序列。

import math

t = [0, 5, 10, 15, 20, 30, 40, 50, 60]
c = [1.002, 0.842, 0.712, 0.598, 0.505, 0.356, 0.252, 0.178, 0.126]

def model(c0, k, ti):
    return c0 * math.exp(-k * ti)

def residuals(c0, k):
    return [ci - model(c0, k, ti) for ti, ci in zip(t, c)]

def solve2(a, b, d, e, f, g):
    # [[a,b],[d,e]] x = [f,g]
    det = a * e - b * d
    return ((f * e - b * g) / det, (a * g - f * d) / det)

c0, k = 1.0, 0.02          # 初值
for it in range(50):
    r = residuals(c0, k)
    # 组装 JtJ 与 Jtr
    A = B = C = 0.0
    g0 = gk = 0.0
    for ti, ri in zip(t, r):
        e = math.exp(-k * ti)
        j0 = e
        jk = -c0 * ti * e
        A += j0 * j0
        B += j0 * jk
        C += jk * jk
        g0 += j0 * ri
        gk += jk * ri
    d0, dk = solve2(A, B, B, C, g0, gk)
    c0 += d0
    k += dk
    if abs(d0) < 1e-10 and abs(dk) < 1e-12:
        break

r = residuals(c0, k)
ss = sum(ri * ri for ri in r)
n, p = len(t), 2
s2 = ss / (n - p)

# 参数协方差 = s2 * (JtJ)^-1
A = B = C = 0.0
for ti in t:
    e = math.exp(-k * ti)
    j0 = e
    jk = -c0 * ti * e
    A += j0 * j0
    B += j0 * jk
    C += jk * jk
det = A * C - B * B
var_c0 = s2 * C / det
var_k = s2 * A / det

print(f"c0 = {c0:.4f} ± {math.sqrt(var_c0):.4f} mol/L")
print(f"k  = {k:.5f} ± {math.sqrt(var_k):.5f} min^-1")
print(f"SSE = {ss:.6f}, s = {math.sqrt(s2):.5f}")
print("残差:", [round(x, 4) for x in r])

示例运行结果（由上述代码直接计算得到）：c₀ ≈ 1.0003 mol/L，k ≈ 0.0345 min⁻¹，残差标准差约 0.006，残差序列在零附近正负交替，无明显单调趋势。

残差随机性判读

参数收敛不等于模型正确。需要检查三件事：残差是否围绕零波动、是否随 t 或拟合值出现系统弯曲、相邻残差是否高度相关。若残差在低浓度段整体偏正、高浓度段整体偏负，说明一级模型可能不足（例如存在扩散限制或反应级数偏离 1）；若残差呈明显“U 形”，则提示 c₀ 与 k 之间存在较强相关性，需要检查取样时刻或初始浓度设定。

结果判断与常见偏差

• 参数标准误：k 的相对标准误若超过约 5%，说明时间跨度或数据精度不足，应延长反应时间或增加中段取样点。
• c₀ 与 k 的相关性：当 c₀ 的估计值明显偏离名义投料浓度时，优先怀疑取样时刻零点偏移，而不是直接接受拟合值。
• 残差趋势：若残差随 t 单调下降，说明实际衰减慢于一级模型，可考虑产物抑制或扩散控制；若残差在末端突然放大，多为低浓度检测限导致的相对误差增大。