目錄

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

1 非線性回歸Python實(shí)現(xiàn)方法

在很多領(lǐng)域中往往需要對(duì)實(shí)驗(yàn)數(shù)據(jù)（或者調(diào)查數(shù)據(jù)）通過某種方程進(jìn)行擬合，得到方程中的關(guān)鍵參數(shù)。在統(tǒng)計(jì)、物理、化學(xué)、生物等領(lǐng)域都是非常常見的方法。
非線性方程是指方程中存在非線性項(xiàng)，例如$e^x$，$\ln x$，$x^3$，$\sin x$等項(xiàng)時(shí)，求解過程通過人工計(jì)算很難得到準(zhǔn)確的解。

1.1 多項(xiàng)式擬合 polyfit()

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

import numpy as np import matplotlib.pyplot as plt# 生成數(shù)據(jù) x = np.linspace(0, 5, 100) y = np.exp(x) # 方差初始容器，輔助內(nèi)容 i_range = [i for i in range(3,10)] rsds = [] # 擬合并繪制不同項(xiàng)數(shù)多項(xiàng)式 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)) # 繪制原始數(shù)據(jù) plt.plot(x[::10], y[::10], 'r*', ms=20, label='exp(x)') plt.legend() # 字體大小設(shè)置 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多項(xiàng)式擬合

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

---

1.2 非線性擬合 curve_fit()

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

from scipy.optimize import curve_fit import numpy as np import matplotlib.pyplot as plt from random import random# 生成數(shù)據(jù) x = np.linspace(0, 5, 100) y = np.exp(x) # 隨機(jī)函數(shù)生成噪音 eps = 0.4*np.array([random() for i in range(len(x))]) - 0.2 y = y * (1 + eps)# 擬合函數(shù)形式 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]) # 多項(xiàng)式擬合 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') # 設(shè)置字體 font = {'size':22} plt.xlabel('x', font) plt.ylabel('y', font) plt.tick_params(labelsize=19) plt.legend(fontsize=20)

圖1-3 非線性擬合+多項(xiàng)式擬合過擬合

圖1-4 非線性擬合
擬合的結(jié)果是a=0.9663855 , n=1.01433087，和預(yù)期的1, 1很接近，擬合效果還可以。多項(xiàng)式擬合當(dāng)使用30次方時(shí)，就產(chǎn)生了過擬合的情況，偏離了實(shí)際的情況，所以多項(xiàng)式擬合還是要慎用，次數(shù)取小了誤差大，取大了會(huì)過擬合。

---

2 應(yīng)用實(shí)例

2.1 擬合活化能

2.1.1 實(shí)驗(yàn)方法與原理

對(duì)于光催化反應(yīng)甲醛反應(yīng)，其動(dòng)力學(xué)方程可表示為
$$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})$$
根據(jù)m取值的不同，濃度與時(shí)間的變化規(guī)律可分為三種形式的函數(shù)
$$\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}$$
通過實(shí)驗(yàn)測(cè)得濃度隨時(shí)間的變化規(guī)律cA(t)，根據(jù)不同的動(dòng)力學(xué)模型對(duì)數(shù)據(jù)進(jìn)行擬合。

圖2-1 實(shí)驗(yàn)濃度變化

2.1.2 結(jié)果與討論

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

圖2-2 速率與濃度變化規(guī)律
擬合得到m=1.75978，k0=24.8696，E=15949.7。擬合曲線與實(shí)驗(yàn)數(shù)據(jù)的對(duì)比如下圖所示。

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

2.1.3 附錄：基于Python實(shí)現(xiàn)

import numpy as np import matplotlib.pyplot as plt from scipy.optimize import curve_fit# 實(shí)驗(yàn)數(shù)據(jù) 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] # 轉(zhuǎn)換為array t = np.array(t[:-1]) c = np.array(c[:-1]) # 反應(yīng)速率 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): '''濃度方程，對(duì)于模型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): '''繪制濃度-時(shí)間圖'''# 生成數(shù)據(jù)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)# 擬合數(shù)據(jù)，繪制圖像 t1 = t[:-1] p1,pcov1 = curve_fit(mfun, c[:-1], r) # very important cf = np.linspace(c[0],c[-1],100) # 細(xì)化數(shù)據(jù)點(diǎn) rf = mfun(cf, p1[0], p1[1], p1[2]) # 計(jì)算擬合曲線 plot(c[:-1], r, cf, rf) # 繪制速率濃度圖 plot2(t,c,p1) # 繪制濃度時(shí)間圖

---

2.2 非理想反應(yīng)器停留時(shí)間分布擬合

2.2.1 引言

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

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

圖3-2 停留時(shí)間分布函數(shù)

2.2.2 研究方法

停留時(shí)間分布規(guī)律通常呈“S”型曲線，F(t)的最大值為1，因此可使用S型函數(shù)表示F(t)，其中b，c為擬合參數(shù)。將實(shí)驗(yàn)數(shù)據(jù)通過S函數(shù)擬合得到停留時(shí)間分布的連續(xù)函數(shù)。
$$F(t)=\frac{1}{1+b\exp (?ct)}$$
因此，進(jìn)一步可解算停留時(shí)間分布密度函數(shù)E(t)，其形式可表示為如下方程：
$$E(t)=\frac{dF(t)}{dt}=\frac{bc\exp (?ct)}{(1+b\exp (?ct){)}^2}$$
平均停留時(shí)間則可表示為如下：
$$\bar{t}={\int }_0^{\infty }tE(t)dt=\frac{\ln (b+1)}{c}$$
停留時(shí)間分布的方差表示為：
$${\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 結(jié)果與討論

表 1 階躍法測(cè)定停留時(shí)間分布實(shí)驗(yàn)數(shù)據(jù)

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

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

圖3-3 擬合停留時(shí)間分布函數(shù)

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

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。