目錄

• 1 非線性回歸Python實現方法
• • 1.1 多項式擬合 polyfit()
  • 1.2 非線性擬合 curve_fit()
• 2 應用實例
• • 2.1 擬合活化能
  • • 2.1.1 實驗方法與原理
    • 2.1.2 結果與討論
    • 2.1.3 附錄：基于Python實現
  • 2.2 非理想反應器停留時間分布擬合
  • • 2.2.1 引言
    • 2.2.2 研究方法
    • 2.2.3 結果與討論

1 非線性回歸Python實現方法

在很多領域中往往需要對實驗數據（或者調查數據）通過某種方程進行擬合，得到方程中的關鍵參數。在統計、物理、化學、生物等領域都是非常常見的方法。
非線性方程是指方程中存在非線性項，例如$e^x$，$\ln x$，$x^3$，$\sin x$等項時，求解過程通過人工計算很難得到準確的解。

1.1 多項式擬合 polyfit()

Python的numpy庫中的函數polyfit()能夠擬合一元多項式
$$y=a_1{x}^m+a_2{x}^{m?1}+...+a_{m}x+a_{m+1}$$numpy.polyfit()專門用于擬合這種形式的多項式，看下面的例子

import numpy as np import matplotlib.pyplot as plt# 生成數據 x = np.linspace(0, 5, 100) y = np.exp(x) # 方差初始容器，輔助內容 i_range = [i for i in range(3,10)] rsds = [] # 擬合并繪制不同項數多項式 for i in i_range:p = np.polyfit(x, y, i, full=True)y_hat = np.polyval(p[0], x)rsds.append(p[1])plt.plot(x, y_hat, lw=2, label=str(i)) # 繪制原始數據 plt.plot(x[::10], y[::10], 'r*', ms=20, label='exp(x)') plt.legend() # 字體大小設置 font = {'size':22} plt.xlabel('x', font) plt.ylabel('y', font) plt.tick_params(labelsize=19) plt.show() # 繪制方差圖 plt.plot(i_range, np.log10(rsds), '-o') plt.xlabel('n', font) plt.ylabel(r'$\lg \sigma^2$', font) plt.tick_params(labelsize=19)

圖1-1多項式擬合

圖1-2 多項式擬合方差隨多項式次數的變化規律
簡單舉個例子，$y=e^x$這個函數也可以擬合成多項式的形式，隨著多項式次數的增加，擬合結果方差逐漸減小。其原理是連續函數的泰勒公式，高數里證明了連續可導的函數總能表示成多項式的形式。
多項式擬合能在一定程度上獲得變化規律，但是也存在著無法忽略的問題，超出原始數據范圍的預測效果差，而且隨著n的增加，有可能出現過擬合的情況（噪音信號被帶入擬合結果）。

---

1.2 非線性擬合 curve_fit()

curve_fit()是scipy.optimize的函數，用于擬合已知形式的函數，接下來還是看個例子。仍然以$y=\exp (x)$為例說明，通過隨機數函數增加生成一些噪音，使噪音范圍為[-0.2, 0.2），然后分別用curve_fit()和polyfit()擬合并對比。用下面的函數擬合，a和n為待擬合參數。
$$y=a\exp (nx)$$
代碼如下：

from scipy.optimize import curve_fit import numpy as np import matplotlib.pyplot as plt from random import random# 生成數據 x = np.linspace(0, 5, 100) y = np.exp(x) # 隨機函數生成噪音 eps = 0.4*np.array([random() for i in range(len(x))]) - 0.2 y = y * (1 + eps)# 擬合函數形式 def f(x, a, n):return a*np.exp(n*x)# 非線性擬合 popt,pcov = curve_fit(f, x, y) y_hat_c = f(x, popt[0], popt[1]) # 多項式擬合 p = np.polyfit(x, y, 6) y_hat_p = np.polyval(p, x) # 繪制圖像 plt.plot(x, y, 'o', label='y') plt.plot(x, y_hat_c, 'r', lw=2, label='curve_fit') plt.plot(x, y_hat_p, 'k--', lw=2, label='polyfit_6') # 設置字體 font = {'size':22} plt.xlabel('x', font) plt.ylabel('y', font) plt.tick_params(labelsize=19) plt.legend(fontsize=20)

圖1-3 非線性擬合+多項式擬合過擬合

圖1-4 非線性擬合
擬合的結果是a=0.9663855 , n=1.01433087，和預期的1, 1很接近，擬合效果還可以。多項式擬合當使用30次方時，就產生了過擬合的情況，偏離了實際的情況，所以多項式擬合還是要慎用，次數取小了誤差大，取大了會過擬合。

---

2 應用實例

2.1 擬合活化能

2.1.1 實驗方法與原理

對于光催化反應甲醛反應，其動力學方程可表示為
$$r=\frac{d{c}_A}{dt}=k_0\exp (?\frac{E}{RT})c_A^m$$
分離變量后得
$$\frac{d{c}_A}{c_A^m}=k_0\exp (?\frac{E}{RT})$$
根據m取值的不同，濃度與時間的變化規律可分為三種形式的函數
$$\begin{gathered} c_A(t)=c_0?k_0\exp (?\frac{E}{RT})t,m=0 \\ \ln c_A(t)=\ln c_0?k_0\exp (?\frac{E}{RT})t,m=1 \\ [c_A(t){]}^{1?m}=c_0^{1?m}?(1?m)k_0\exp (?\frac{E}{RT})t,m=0,1 \end{gathered}$$
通過實驗測得濃度隨時間的變化規律cA(t)，根據不同的動力學模型對數據進行擬合。

圖2-1 實驗濃度變化

2.1.2 結果與討論

從實驗數據初步判斷m≠0，接下來將使用模型2和模型3對實驗數據進行擬合。
$$\ln r=m\ln c_0+\ln k_0?E/RT$$
反應速率的表達式更為統一，對反應速率方程取對數后得到的方程形式更簡潔，因此為減少在擬合過程中不必要的誤差，可使用速率方程進行擬合，計算過程見附錄，擬合結果如下。

圖2-2 速率與濃度變化規律
擬合得到m=1.75978，k0=24.8696，E=15949.7。擬合曲線與實驗數據的對比如下圖所示。

圖2-3 濃度隨時間變化規律
通過非線性回歸方法計算得出甲醛光催化反應的反應級數為1.76，活化能為15.9 kJ/mol。

2.1.3 附錄：基于Python實現

import numpy as np import matplotlib.pyplot as plt from scipy.optimize import curve_fit# 實驗數據 t = [0,10,15,25,35,45,105,115,125,150,165,180,190,210,225,240,255,270] c = [0.29, 0.26, 0.23, 0.2, 0.18, 0.16, 0.11, 0.1, 0.09, 0.09,0.07, 0.06, 0.06, 0.06, 0.06, 0.06, 0.06, 0.06] # 轉換為array t = np.array(t[:-1]) c = np.array(c[:-1]) # 反應速率 r = -np.diff(c)/np.diff(t) T = 298def mfun(c, m, k0, E): ''' 速率方程'''lnr = m*np.log(c) + np.log(k0) - E/2477.6return np.exp(lnr)def mfun2(t, k0, E, m): '''濃度方程，對于模型3'''cm = 0.29**(1-m) - (1-m)*k0*np.exp(-E/2477.6)*treturn cm # 返回值為 c**(1-m)def plot(x1, y1, x2, y2): '''繪制速率-濃度圖'''font = {'size':22}plt.plot(x1, y1, 'o')plt.plot(x2, y2)#plt.ylim([0,0.35])plt.xlabel('c / ppm',font)plt.ylabel('r / (ppm/s)',font)plt.tick_params(labelsize=19)plt.show()def plot2(t1, c1, popt): '''繪制濃度-時間圖'''# 生成數據font = {'size':22}t = np.linspace(t1[0],t1[-1],100)cm = mfun2(t, popt[1], popt[2], popt[0])**(1/(1-popt[0]))# 繪圖plt.plot(t1, c1, 'o')plt.plot(t, cm)plt.ylim([0,0.35])plt.xlabel('t / sec',font)plt.ylabel('c(t) / ppm',font)plt.tick_params(labelsize=19)# 擬合數據，繪制圖像 t1 = t[:-1] p1,pcov1 = curve_fit(mfun, c[:-1], r) # very important cf = np.linspace(c[0],c[-1],100) # 細化數據點 rf = mfun(cf, p1[0], p1[1], p1[2]) # 計算擬合曲線 plot(c[:-1], r, cf, rf) # 繪制速率濃度圖 plot2(t,c,p1) # 繪制濃度時間圖

---

2.2 非理想反應器停留時間分布擬合

2.2.1 引言

非理想反應器內流體存在返混現象，因此流體微元在反應器中的停留時間并不是完全相同的，而存在與反應器特性有關的概率分布，稱為停留時間分布。停留時間分布對于認識反應器傳質過程與評估反應器生產能力等都有重要參考價值。
停留時間分布的測試方法分為兩類：階躍法和脈沖法。階躍法是指在某一時刻突然改變穩定運行的反應器入口的探針物質濃度，通過監測不同時間后出口處探針物質濃度的變化規律，直接得到離散的停留時間分布F(t)。脈沖法是指在某一瞬間，從入口注入一定量探針物質，通過監測出口濃度，直接得到停留時間分布密度E(t)。

圖3-1 階躍法停留時間分布示意圖
停留時間分布函數F(t)與停留時間分布密度函數E(t)體現了反應器的某些工作特點，比如返混程度、平均停留時間等等。對于理想平推流反應器無返混，而理想全混釜反應器返混無窮大，通常非理想反應器的返混程度則是在兩者之間。可見停留時間分布函數對于評估反應器生產能力、設計反應器體積等極為重要。

圖3-2 停留時間分布函數

2.2.2 研究方法

停留時間分布規律通常呈“S”型曲線，F(t)的最大值為1，因此可使用S型函數表示F(t)，其中b，c為擬合參數。將實驗數據通過S函數擬合得到停留時間分布的連續函數。
$$F(t)=\frac{1}{1+b\exp (?ct)}$$
因此，進一步可解算停留時間分布密度函數E(t)，其形式可表示為如下方程：
$$E(t)=\frac{dF(t)}{dt}=\frac{bc\exp (?ct)}{(1+b\exp (?ct){)}^2}$$
平均停留時間則可表示為如下：
$$\bar{t}={\int }_0^{\infty }tE(t)dt=\frac{\ln (b+1)}{c}$$
停留時間分布的方差表示為：
$${\sigma }_t^2={\int }_0^{\infty }(t?\bar{t}{)}^{2}E(t)dt={\int }_0^{\infty }{t}^{2}E(t)dt?{\bar{t}}^2$$

2.2.3 結果與討論

表 1 階躍法測定停留時間分布實驗數據

t/s0152535455565758595

Cout/(kg m-3)	0	0.5	1.0	2.0	4.0	5.5	6.5	7.0	7.7	7.7

（1）通過離散方法計算得到=50.58，σt2=341.54。
（2）通過擬合S型函數得到停留時間分布函數：
$$F(t)=\frac{1}{1+63.5564\exp (?0.0916t)}$$
停留時間分布密度函數：
$$E(t)=\frac{5.8218\exp (?0.0916t)}{(1+63.5564\exp (?0.0916t){)}^2}$$
平均停留時間：
$$\bar{t}=\frac{\ln (63.5564+1)}{0.0916}=45.50$$

停留時間分布方差：
$${\sigma }_t^2=364.60$$

圖3-3 擬合停留時間分布函數

*圖3-4 擬合停留時間分布密度函數
擬合函數的計算結果與離散方法的計算結果相比，數值相差不超過10%，通過圖像可以看到，S函數能較好的擬合出停留時間分布函數。經過多次嘗試發現，S函數擬合返混程度較大的非理想反應器時效果較好，而對接近平推流的反應器擬合效果差。

總結

以上是生活随笔為你收集整理的反应动力学参数拟合与停留时间分布函数——基于Python实现的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。