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实现方法

在很多领域中往往需要对实验数据（或者调查数据）通过某种方程进行拟合，得到方程中的关键参数。在统计、物理、化学、生物等领域都是非常常见的方法。
非线性方程是指方程中存在非线性项，例如exe^xex，ln⁡x\ln xlnx，x3x^3x3，sin⁡x\sin xsinx等项时，求解过程通过人工计算很难得到准确的解。

1.1 多项式拟合 polyfit()

Python的numpy库中的函数polyfit()能够拟合一元多项式
y=a1xm+a2xm−1+...+amx+am+1y = a_1 x^m + a_2 x^{m-1} +...+ a_m x + a_{m+1} y=a1xm+a2xm−1+...+amx+am+1numpy.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=exy=e^xy=ex这个函数也可以拟合成多项式的形式，随着多项式次数的增加，拟合结果方差逐渐减小。其原理是连续函数的泰勒公式，高数里证明了连续可导的函数总能表示成多项式的形式。
多项式拟合能在一定程度上获得变化规律，但是也存在着无法忽略的问题，超出原始数据范围的预测效果差，而且随着n的增加，有可能出现过拟合的情况（噪音信号被带入拟合结果）。

1.2 非线性拟合 curve_fit()

curve_fit()是scipy.optimize的函数，用于拟合已知形式的函数，接下来还是看个例子。仍然以y=exp⁡(x)y=\exp(x)y=exp(x)为例说明，通过随机数函数增加生成一些噪音，使噪音范围为[-0.2, 0.2），然后分别用curve_fit()和polyfit()拟合并对比。用下面的函数拟合，a和n为待拟合参数。
y=aexp⁡(nx)y = a\exp(nx)y=aexp(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=dcAdt=k0exp⁡(−ERT)cAmr = \frac{dc_A}{dt} = k_0\exp(-\frac{E}{RT})c_A^mr=dtdcA=k0exp(−RTE)cAm
分离变量后得
dcAcAm=k0exp⁡(−ERT)\frac{dc_A}{c_A^m} = k_0\exp(-\frac{E}{RT})cAmdcA=k0exp(−RTE)
根据m取值的不同，浓度与时间的变化规律可分为三种形式的函数
cA(t)=c0−k0exp⁡(−ERT)t,m=0ln⁡cA(t)=ln⁡c0−k0exp⁡(−ERT)t,m=1[cA(t)]1−m=c01−m−(1−m)k0exp⁡(−ERT)t,m≠0,1c_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\neq 0,1 cA(t)=c0−k0exp(−RTE)t, m=0lncA(t)=lnc0−k0exp(−RTE)t, m=1[cA(t)]1−m=c01−m−(1−m)k0exp(−RTE)t, m=0,1
通过实验测得浓度随时间的变化规律cA(t)，根据不同的动力学模型对数据进行拟合。

图2-1 实验浓度变化

2.1.2 结果与讨论

从实验数据初步判断m≠0，接下来将使用模型2和模型3对实验数据进行拟合。
ln⁡r=mln⁡c0+ln⁡k0−E/RT\ln r = m\ln c_0 + \ln k_0 - E/RTlnr=mlnc0+lnk0−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)=11+bexp⁡(−ct)F(t) = \frac{1}{1+b\exp(-ct)}F(t)=1+bexp(−ct)1
因此，进一步可解算停留时间分布密度函数E(t)，其形式可表示为如下方程：
E(t)=dF(t)dt=bcexp⁡(−ct)(1+bexp⁡(−ct))2E(t) = \frac{dF(t)}{dt} = \frac{bc\exp(-ct)}{(1+b\exp(-ct))^2}E(t)=dtdF(t)=(1+bexp(−ct))2bcexp(−ct)
平均停留时间则可表示为如下：
tˉ=∫0∞tE(t)dt=ln⁡(b+1)c\bar t = \int_0^\infty tE(t) dt = \frac{\ln (b+1)}{c}tˉ=∫0∞tE(t)dt=cln(b+1)
停留时间分布的方差表示为：
σt2=∫0∞(t−tˉ)2E(t)dt=∫0∞t2E(t)dt−tˉ2\sigma^2_t = \int_0^\infty (t-\bar t)^2E(t)dt = \int_0^\infty t^2E(t)dt - \bar t^2σt2=∫0∞(t−tˉ)2E(t)dt=∫0∞t2E(t)dt−tˉ2

2.2.3 结果与讨论

表 1 阶跃法测定停留时间分布实验数据

t/s	0	15	25	35	45	55	65	75	85	95
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)=11+63.5564exp⁡(−0.0916t)F(t) = \frac{1}{1+63.5564\exp(-0.0916t)}F(t)=1+63.5564exp(−0.0916t)1
停留时间分布密度函数：
E(t)=5.8218exp⁡(−0.0916t)(1+63.5564exp⁡(−0.0916t))2E(t) = \frac{5.8218\exp(-0.0916t)}{(1+63.5564\exp(-0.0916t))^2}E(t)=(1+63.5564exp(−0.0916t))25.8218exp(−0.0916t)
平均停留时间：
tˉ=ln⁡(63.5564+1)0.0916=45.50\bar t = \frac{\ln (63.5564+1)}{0.0916} = 45.50tˉ=0.0916ln(63.5564+1)=45.50

停留时间分布方差：
σt2=364.60\sigma_t^2 = 364.60σt2=364.60

图3-3 拟合停留时间分布函数

*图3-4 拟合停留时间分布密度函数
拟合函数的计算结果与离散方法的计算结果相比，数值相差不超过10%，通过图像可以看到，S函数能较好的拟合出停留时间分布函数。经过多次尝试发现，S函数拟合返混程度较大的非理想反应器时效果较好，而对接近平推流的反应器拟合效果差。

反应动力学参数拟合与停留时间分布函数——基于Python实现相关推荐

基于python的AD-census立体匹配算法实现
文章目录前言一.AD-census是什么? 1.代价计算 2.代价聚合 3.视差优化 4.视差后处理二.基于python的AD-census立体匹配算法实现前言 AD-Census算法来自 ...
python 物理实验_基于Python和梯度下降算法的物理实验数据一元线性拟合方法
基于 Python 和梯度下降算法的物理实验数据一元线性拟合方法关毅铬 ; 程敏熙 [期刊名称] < <物理通报> > [年 ( 卷 ), 期] 2019(000)010 ...
基于Python的随机森林（RF）回归与多种模型超参数自动优化方法
本文详细介绍基于Python的随机森林(Random Forest)回归算法代码与模型超参数(包括决策树个数与最大深度.最小分离样本数.最小叶子节点样本数.最大分离特征数等等)自动优化代码. ...
python深度神经网络量化_基于Python建立深度神经网络！你学会了嘛？
原标题:基于Python建立深度神经网络!你学会了嘛? 图1 神经网络构造的例子(符号说明:上标[l]表示与第l层:上标(i)表示第i个例子:下标i表示矢量第i项) 单层神经网络图2 单层神经网络示 ...
基于Python的卷积神经网络和特征提取
基于Python的卷积神经网络和特征提取发表于2015-08-27 21:39| 4577次阅读| 来源blog.christianperone.com/| 13 条评论| 作者Christi ...
ML之XGBoost：XGBoost参数调优的优秀外文翻译—《XGBoost中的参数调优完整指南(带python中的代码)》(四)
ML之XGBoost:XGBoost参数调优的优秀外文翻译-<XGBoost中的参数调优完整指南(带python中的代码)>(四) 目录 Step 3: Tune gamma步骤3:伽马微 ...
ML之XGBoost：XGBoost参数调优的优秀外文翻译—《XGBoost中的参数调优完整指南(带python中的代码)》(二)
ML之XGBoost:XGBoost参数调优的优秀外文翻译-<XGBoost中的参数调优完整指南(带python中的代码)>(二) 目录 2. xgboost参数/XGBoost Para ...
TF学习：Tensorflow基础案例、经典案例集合——基于python编程代码的实现
TF学习:Tensorflow基础案例.经典案例集合--基于python编程代码的实现目录 Tensorflow的使用入门 1.TF:使用Tensorflow输出一句话 2.TF实现加法 3.TF实 ...
基于Python的特征自动化选择：两行代码完成特征工程
本文介绍一个特征选择神器:特征选择器是用于减少机器学习数据集的维数的工具,可以傻瓜式地进行特征选择,两行代码即可搞定!! 来源:Will Koehrsen 代码整理及注释翻译:黄海广代码和数据下载地 ...

反应动力学参数拟合与停留时间分布函数——基于Python实现

目录