非绝热连续搅拌槽式反应器的建模

连续搅拌槽式反应器(CSTR)是过程工业中较为常见的一种化学系统。在这里，我们将研究在Bequette的《过程动力学:建模、分析和模拟》(Prentice-Hall出版，1998)一书中广泛描述的加套非绝热(即非绝热)罐式反应器。假设容器完全混合，发生单一的一级放热和不可逆反应，a -> B。容器和周围冷却夹套的示意图显示在绘图窗口中。注意这是一个草图;在现实中，冷却剂流动，例如，通常围绕着整个反应堆外壳，而不仅仅是底部。

图1:CSTR原理图。

更先进的控制方法需要CSTR模型。试剂A的进口流以恒定的速率f进料。搅拌后，最终产物以与试剂A进料槽相同的速率流出容器(因此反应器槽中的体积V保持恒定)。该控制策略要求控制夹套温度u_3(t)，以使试剂A y_1(t)的浓度保持在期望的水平，而不受进料流浓度和温度(输入分别为u_1(t)和u_2(t))的干扰。由于y_2(t)罐内的温度在反应器运行过程中可能发生显著变化，因此，也需要确保该过程变量保持在合理的范围内。

对装运箱建模

CSTR系统使用基本的会计和节能原理建模。容器内试剂A浓度每时间单位d C_A(t)/dt (= d y_1(t)/dt)的变化可表示为:

d C_A (t ) -------- = F / V * (C_Af (t) -C_A (t) - r (t) dt

其中第一项表示由于试剂A在进口流和容器中的浓度不同而引起的浓度变化，第二项表示由于容器中的化学反应而引起的浓度变化(反应速率)。单位体积的反应速率用阿伦尼乌斯速率定律描述:

r (t) = k_0 * exp (e / (r * t (t))) * C_A (t)

它表明化学反应的速率随绝对温度呈指数增长。k_0在这里是一个未知的非热常数，E是活化能，R玻尔兹曼理想气体常数，T(T) (= y_2(T))是反应堆的温度。

同样，利用能量平衡原理(假设反应器内体积恒定)，反应器内每时间单位dt (T)/dt的温度变化可建模为:

d T (T ) ------ = F / V (T_f (t) - t (t)) - (H / c_p *ρ)* r (t) - (U *) / (c_p *ρ* V) * (t (t) -T_j (t) / dt

其中第一项和第三项描述了由于进料流温度T_f(t)和夹套冷却剂温度T_j(t)分别与反应堆温度不同而发生的变化。第二项是容器内化学反应对反应堆温度的影响。在这个方程中，H是反应热参数，c_p是热容项，rho是密度项，U是总传热系数，a是换热面积(冷却剂/容器面积)。

总之，CSTR有三个输入信号:

u_1（t）=入口进料流中A的C_Af（t）浓度[kgmol/m^3]。u_2（t）=t_f（t）入口进料流温度[K]。u_3（t）=t_j（t）夹套冷却剂温度[K]。

和两个输出信号:

y_1（t）=反应器罐中A的C_A（t）浓度[kgmol/m^3]。y_2（t）=t（t）反应器温度[K]。

这也是自然模型状态，即y_1(t) = x_1(t)和y_2(t) = x_2(t)

将一些原始参数集中在一起后，我们得到了八个不同的模型参数：

F体积流量(体积/时间)[m^3/h]。固定的。V反应器体积[m^3]。固定的。k_0指前非热因子[1/h]。免费的。E活化能[千卡/kgmol]。免费的。R玻尔兹曼气体常数[千卡/(kgmol*K)]。固定的。 H Heat of reaction [kcal/kgmol]. Fixed. HD = c_p*rho Heat capacity times density [kcal/(m^3*K)]. Free. HA = U*A Overall heat transfer coefficient times tank area [kcal/(K*h)] Free.

这里指定了四个参数是自由的(即，待估计的)。在实际应用中，还可以通过实验确定指前非热因子k_0和活化能E。这将简化识别任务，只考虑两个未知数:热容c_p和总传热系数U(分别由HD和HA推导，即已知rho和A)。

通过上述引入的符号，获得了CSTR的以下状态空间表示。

d x_1（t）--=F/V*（u_1（t）-x_1（t））-k_0*exp（-E/（R*x_2（t）））*x_1（t）dt

d x_2（t）-------------=F/V（u_2（t）-x_2（t））-（H/HD）*k_0*exp（-E/（R*x_2（t））*x_1（t）dt-（HA/（HD*V））*（x_2（t）-u 3（t））

y_1（t）=x_1（t）y_2（t）=x_2（t）

创建非绝热连续搅拌槽式反应器IDNLGREY对象

该信息输入到名为cstr_m.m的MATLAB文件中，包含以下内容。

函数[dx，y]=cstr_m（t，x，u，F，V，k_0，E，R，H，HD，HA，varargin）%cstr_m非绝热连续搅拌釜式反应器（cstr）。

%输出方程。y = [x (1);．．．反应器中A物质的浓度%。x(2)……%反应堆温度。］;

%状态方程。dx = [F / V * (u (1) - x (1)) -k_0 * exp (e / (R * x (2))) * x (1);．．．F / V * (u (2) - x(2))——(H / HD) * k_0 * exp (e / (R * x (2))) * x (1) - (HA /(高清* V)) * (x (2) - u(3))…］;

接下来创建一个反映建模情况的IDNLGREY对象。

文件名=“cstr_m”；%描述模型结构的文件。订单= [2 3 2];模型订单[ny nu nx]。参数=[1；1；35e6；11850；...%初始参数。1.98589;-5960;480;145);InitialStates = (8.5695;311.267);%初始状态的初始值。t = 0;%时间连续系统。nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts，“名字”，...搅拌釜反应器的，...“TimeUnit”，“小时”）;

CSTR模型结构的输入、状态和输出是使用SET和setit方法指定的。我们还指定默认情况下应该估计初始状态。

nlgr。InputName = {进料流中A的浓度...%u（1）。进料流温度...% u(2)。夹克冷却液温度的};% u(3)。nlgr.InputUnit={‘kgmol/m^3’“K”“K”};nlgr = setinit (nlgr,“名字”, {反应堆储罐中A的浓度...% x(1)。“反应堆温度”}）;...% x(2)。nlgr = setinit (nlgr,“单位”, {‘kgmol/m^3’“K”}）;nlgr = setinit (nlgr,“固定的”,{假假});nlgr。OutputName = {“集中注意力”...% y (1);反应池中A的浓度“反应堆温度。”};%(2)。nglr。OutputUnit = {‘kgmol/m^3’“K”};

定义了CSTR模型结构的参数，并指定F、V、R和H固定。

nlgr=setpar（nlgr，“名字”, {“体积流率（体积/时间）”...%F。“反应堆体积”...%五,。“Pre-exponential非热能的因素”...%k_0。“活化能”...%E。玻耳兹曼理想气体常数...% R。反应热的...%H。“热容乘以密度”...%高清。‘总传热系数乘以储罐面积’}）;...%公顷。nlgr=setpar（nlgr，“单位”, {“m ^ 3 /小时”“m ^ 3”1 / h的“千卡/千克摩尔”“千卡/ (kgmol * K)”...“千卡/千克摩尔”“千卡/ (m ^ 3 * K)”“千卡/ (K * h)”})；nlgr.参数（1）.固定=真；%修正F。nlgr.Parameters(2)。固定= true;% V修复。nlgr.参数（5）.Fixed=true；%修复R。nlgr.Parameters(6)。固定= true;%固定H。

通过物理推理，我们也知道除了反应热参数(总是负的，因为反应是放热的)外，其他参数都是正的。让我们将这些(有点粗糙的)知识整合到我们的CSTR模型结构中:

nlgr.Parameters(1)。最小值= 0;%F。nlgr.Parameters(2)。最小值= 0;%五,。nlgr.Parameters(3)。最小值= 0;%k_0。nlgr.Parameters(4)。最小值= 0;%E。nlgr.Parameters(5)。最小值= 0;% R。nlgr.Parameters(6)。最大= 0;%H。nlgr.参数（7）。最小值=0；%高清。nlgr.参数（8）。最小值=0；%公顷。

接下来通过PRESENT命令得到输入的CSTR模型结构的概要:

礼物(nlgr);

nlgr = 'cstr_m'定义的连续时间非线性灰箱模型(MATLAB文件):dx/dt = F(t, u(t)， x(t)， p1，…y(t) = H(t, u(t)， x(t)， p1，…，p8) + e(t) with 3 input(s), 2 state(s), 2 output(s), and 4 free parameter(s) (out of 8). Inputs: u(1) Concentration of A in inlet feed stream(t) [kgmol/m^3] u(2) Inlet feed stream temperature(t) [K] u(3) Jacket coolant temperature(t) [K] States: Initial value x(1) Concentration of A in reactor tank(t) [kgmol/m^3] xinit@exp1 8.5695 (estimated) in [-Inf, Inf] x(2) Reactor temperature(t) [K] xinit@exp1 311.267 (estimated) in [-Inf, Inf] Outputs: y(1) A Concentration(t) y(2) Reactor temp.(t) Parameters: Value p1 Volumetric flow rate (volume/time) [m^3/h] 1 (fixed) in [0, Inf] p2 Volume in reactor [m^3] 1 (fixed) in [0, Inf] p3 Pre-exponential nonthermal factor [1/h] 3.5e+07 (estimated) in [0, Inf] p4 Activation energy [kcal/kgmol] 11850 (estimated) in [0, Inf] p5 Boltzmann's ideal gas constant [kcal/(kgmo..] 1.98589 (fixed) in [0, Inf] p6 Heat of reaction [kcal/kgmol] -5960 (fixed) in [-Inf, 0] p7 Heat capacity times density [kcal/(m^3*K)] 480 (estimated) in [0, Inf] p8 Overall heat transfer coefficient times tank area [kcal/(K*h)] 145 (estimated) in [0, Inf] Name: Stirred tank reactor Status: Created by direct construction or transformation. Not estimated. More information in model's "Report" property.

输入输出数据

许多非线性系统的系统辨识实验设计通常比线性系统要复杂得多。CSTR也是如此，一方面希望可控输入u_3能够充分激励系统，另一方面必须选择“工厂友好型”（化学过程必须保持稳定，试验持续时间应尽可能短，以便对生产的影响最小，等等）讨论CSTR输入信号的选择。有人认为多正弦输入u_3有利于多级伪随机输入信号，原因有几个。在下面的识别实验中，我们将使用两个这样的输入信号，一个用于估计，一个用于验证，它们是通过MATLA生成的B®输入数据生成工具（GUI），由上述文章的作者提供。

我们加载CSTR数据，并将其放在两个不同的IDDATA对象中，ze用于估计，zv用于验证:

负载(fullfile (matlabroot“工具箱”，“识别”，“iddemos”，“数据”，“cstrdata”））;t = 0.1;每小时% 10个样品!y = 1, y = 1, y = 1“名字”，“估计数据”）;泽。InputName = nlgr.InputName;泽。InputUnit = nlgr.InputUnit;泽。OutputName = nlgr.OutputName;泽。OutputUnit = nlgr.OutputUnit;泽。Tstart = 0; ze.TimeUnit =“小时”；泽。ExperimentName =“估计”；zv=iddata（y2，u2，0.1，“名字”，验证数据的）;zv。InputName = nlgr.InputName;zv。InputUnit = nlgr.InputUnit;zv。OutputName = nlgr.OutputName;zv。OutputUnit = nlgr.OutputUnit;zv。Tstart = 0; zv.TimeUnit =“小时”；zv.name=“验证”；

估计数据集ze的输入和输出显示在两个图中。

图(“名字”，[ze.Name“：输入数据”]);对于我= 1:泽。ν次要情节(泽。ν,1,);情节(泽。SamplingInstants ze.InputData (:, i));标题([“输入#”num2str(我)“:”ze.InputName{我}]);包含(”)；轴线紧；结束包含([泽。域“(”泽。TimeUnit“)”]);

图2:向CSTR输入的估计数据集。

图(“名字”，[ze.Name“：输出数据”]);对于我= 1:泽。纽约次要情节(泽。纽约,1,我);情节(泽。SamplingInstants ze.OutputData (:, i));标题([“输出#”num2str(我)“:”ze.OutputName{我}]);包含(”)；轴线紧；结束包含([泽。域“(”泽。TimeUnit“)”]);

图3:从CSTR输出的估计数据集。

之前进行识别实验,我们应该提到生成的输入信号迫使装运箱的输出显示很多反应堆非线性整体(80度左右的温度变化,导致一些反应堆的“点火”现象明显)。虽然这激发了反应堆(通常从识别的角度来看是好的)，但这可能不是工程师们想要操作一个真实的反应堆的方式，尤其是一个像这个一样放热的反应堆。使用布劳恩等人(2002)所描述的指导方针，人们可以在实际进行实验之前重新设计实验。在这种情况下，尝试减少实验的持续时间并使用周期较短的多正弦输入信号将是有趣的。其目的是减少可控输入信号的低频含量，从而减少电抗器输出的变化。

初始CSTR模型的性能

初始CSTR模型有多好？让我们通过使用ze和zv的输入信号模拟初始模型来研究这一点，并将计算输出与ze和zv中包含的真实输出（通过使用其他参数模拟上述IDNLGREY模型并添加一些噪声获得）分别进行比较。

clf比较(泽nlgr);图;比较(zv nlgr);

图4:初始CSTR模型的真实输出和模拟输出的比较。

参数估计

初始CSTR模型的真实输出和模拟输出之间的一致性良好。为了进一步改进它，我们使用估计数据集ze估计了四个自由模型参数以及模型的初始状态向量。我们指示NLGREYEST显示迭代信息，并执行最多25次搜索迭代。

选择= nlgreyestOptions (“显示”，“开”）;opt.SearchOptions.MaxIterations = 25;NLGR = nlgreyest(ze, NLGR, opt);

如果需要，可以通过第二次调用NLGREYEST（或PEM）继续搜索。这一次NLGREYEST被指示不显示任何迭代信息，并且最多只执行5次迭代。

opt.Display =“关闭”；opt.SearchOptions.MaxIterations=5；nlgr=NLGRYEST（ze，nlgr，opt）；

估计CSTR模型的性能

为了评估估计模型的性能，我们再次使用COMPARE:

比较(泽nlgr);图;比较(zv nlgr);

图5:估计的CSTR模型的真实输出和模拟输出之间的比较。

目视检查立即显示，估计模型的输出接近真实的输出，无论是对ze和zv。这种改进对于验证数据集尤其重要，其中，对于两个模型输出，模型拟合值从负值分别增加到大约70%和99%。

PRESENT将返回关于估计的CSTR模型的进一步信息:

礼物(nlgr);

nlgr = 'cstr_m'定义的连续时间非线性灰箱模型(MATLAB文件):dx/dt = F(t, u(t)， x(t)， p1，…y(t) = H(t, u(t)， x(t)， p1，…，p8) + e(t) with 3 input(s), 2 state(s), 2 output(s), and 4 free parameter(s) (out of 8). Inputs: u(1) Concentration of A in inlet feed stream(t) [kgmol/m^3] u(2) Inlet feed stream temperature(t) [K] u(3) Jacket coolant temperature(t) [K] States: Initial value x(1) Concentration of A in reactor tank(t) [kgmol/m^3] xinit@exp1 8.62914 (estimated) in [-Inf, Inf] x(2) Reactor temperature(t) [K] xinit@exp1 311.215 (estimated) in [-Inf, Inf] Outputs: y(1) A Concentration(t) y(2) Reactor temp.(t) Parameters: Value Standard Deviation p1 Volumetric flow rate (volume/time) [m^3/h] 1 0 (fixed) in [0, Inf] p2 Volume in reactor [m^3] 1 0 (fixed) in [0, Inf] p3 Pre-exponential nonthermal factor [1/h] 3.55889e+07 17548.3 (estimated) in [0, Inf] p4 Activation energy [kcal/kgmol] 11853.9 0.0703052 (estimated) in [0, Inf] p5 Boltzmann's ideal gas constant [kcal/(kgmo..] 1.98589 0 (fixed) in [0, Inf] p6 Heat of reaction [kcal/kgmol] -5960 0 (fixed) in [-Inf, 0] p7 Heat capacity times density [kcal/(m^3*K)] 500.71 0.194139 (estimated) in [0, Inf] p8 Overall heat transfer coefficient times tank area [kcal/(K*h)] 150.127 0.0697455 (estimated) in [0, Inf] Name: Stirred tank reactor Status: Termination condition: Maximum number of iterations or number of function evaluations reached.. Number of iterations: 6, Number of function evaluations: 7 Estimated using Solver: ode45; Search: lsqnonlin on time domain data "Estimation data". Fit to estimation data: [71.36;99.18]% FPE: 0.02736, MSE: 0.6684 More information in model's "Report" property.

在Simulink中对估计CSTR模型进行仿真金宝app

还可以导入IDNLGREY模型并在Simulink中使用。金宝appSimu金宝applink模型“cstr_sim”导入验证数据集zv，并将其数据传递给Simulink IDNLGREY模型块，该模型块在模拟时产生输出，与使用的输入信号一起存储在MATLAB工作空间中的IDDATA对象zsim中。(下面的最后五行代码用于确保向Simulink模型提供适当的模型输入。金宝app只有当idnlgreydemo9在函数中运行时才需要这样做，因为无法保证对zv和nlgr的访问。)

open_system (“cstr_sim”）;如果~ evalin (“基地”，的存在(“zv”、“var”))cstrws=get_参数（bdroot，“modelworkspace”）;cstrws.assignin (“zv”zv);cstrws.assignin (“nlgr”, nlgr);结束

图6:金宝appSimulink模型包含估计的CSTR模型。

通用IDNLGREY Simulink库块金宝app可以在标准系统识别Simulink库中找到，可以复制到任何Simulink模型中使用。例如，在CSTR的情况下，它可以很好地用于闭环控制安排。

在模拟IDNLGREY块之前，必须先配置IDNLGREY块。这可以通过输入持有IDNLGREY模型(这里是nlgr)的MATLAB工作区变量来实现，或者通过在标签为“IDNLGREY模型”的参数编辑框中使用IDNLGREY构造函数定义一个合适的IDNLGREY模型对象来实现。这里还可以指定要使用的初始状态向量(默认值是指定IDNLGREY对象的内部存储的初始状态向量)。

IDNLGREY模型对象存储SIM、PREDICT、NLGREYEST等所使用的微分/差分方程求解器的设置属性(参见nlgr。SimulationOptions用于控制模型的仿真选项)。在Si金宝appmulink中，这些设置总是被重写，以便使用Simulink指定解算器的选项。如果IDNLGREY模型对象指定并使用不同的求解器，那么在Simulink中得到的仿真结果可能与在MATLAB中使用IDNLGREY/SIM得到的仿真结果不同。金宝app示例“idnlgreydemo10”提供了一个示例来说明这一点。

解决了求解器选项后，我们接下来可以执行cstr_sim Simulink模型的命令提示符仿真(这里将对估计的CSTR模型进行仿真)。金宝app(如果idnlgreydemo9在函数中运行，则需要evalin调用来检索zsim。)

sim卡(“cstr_sim”）;如果~存在(“zsim”，“var”) zsim =“基地”，“zsim”）;结束zsim。InputName = nlgr.InputName;zsim。InputUnit = nlgr.InputUnit;zsim。OutputName = nlgr.OutputName;zsim。OutputUnit = nlgr.OutputUnit;时间单位=“小时”；

最后，我们通过绘制Simulink获得的仿真结果来总结示例。金宝app

图(“名字”，‘金宝app估计CSTR模型的Simulink仿真结果’）;对于我= 1:zsim。纽约次要情节(zsim。纽约,1,我);情节(zsim。SamplingInstants ze.OutputData (:, i));标题([“输出#”num2str(我)“:”zsim.OutputName{我}]);包含(”)；轴线紧；结束xlabel（[zsim.Domain“(”时间单位“)”]);

图7:在Simulink中模拟估计的CSTR模型得到的输出。金宝app

结论

本教程涵盖了非绝热连续搅拌槽式反应器的建模和识别。特别是，它说明了如何在Simulink中导入和使用IDNLGREY模型。金宝app

参考文献

[1] Braun, M.W.， R. ortizo - mojica和D.E. Rivera，“最小波峰因子多正弦信号在非线性过程系统的‘植物友好’识别中的应用”，刊于控制工程实践，第10卷，第3期，2002年，301-313页。