非绝热连续搅拌釜反应器:MATLAB文件建模与仿真在Simulink®金宝app

这个例子使用了:

打开脚本

这个例子展示了如何在Simulink®中包含和模拟IDNLGREY模型。金宝app我们使用化学反应系统作为建模基础。示例的第一个建模和识别部分可以在没有Simulink的情况下运行。金宝app

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

在过程工业中遇到的一个相当常见的化学系统是连续搅拌槽反应器(CSTR)。在这里，我们将研究一种夹套非绝热(即非绝热)罐式反应堆，Bequette在1998年由Prentice-Hall出版的《过程动力学:建模、分析和仿真》一书中广泛描述过。假设容器完全混合，发生单次一级放热不可逆反应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建模

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))是反应堆的温度。

同样，利用能量平衡原理(假设反应器内体积恒定)，反应器内每时间单位d T(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) = C_Af(t)进料流A的浓度[kgmol/m^3]。u_2(t) = T_f(t)进料流温度[K]。u_3(t) = T_j(t)护套冷却液温度[K]。

和两个输出信号:

y_1(t) = C_A(t)反应器槽内A的浓度[kgmol/m^3]。y_2(t) = t (t)反应堆温度[K]。

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

在将一些原始参数集中在一起之后，我们最终得到了8个不同的模型参数:

F体积流速(体积/时间)[m^3/h]。固定的。V反应器容积[m^3]。固定的。k_0指数前非热因子[1/h]。免费的。E活化能[kcal/kgmol]。免费的。R玻尔兹曼气体常数[kcal/(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 /(高清* V)) * (x_2 (t) -u_3 (t))

y (t) = x (t) y (t) = x (t)

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

该信息被输入到一个名为cstr_m的MATLAB文件中。M与以下内容。

function [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;35 e6;11850;.．.初始参数。1.98589;-5960;480;145);InitialStates = [8.5695;311.267);初始状态的初始值。Ts = 0;%时间连续系统。nlgr = idnlgrey(文件名，顺序，参数，InitialStates, Ts，“名字”，.．.“搅拌釜式反应器”，.．.“TimeUnit”，“小时”）;

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

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，“固定”， {false false});nlgr。OutputName = {“浓度”.．.% y (1);反应堆罐内A的浓度“反应堆温度。”};%(2)。nglr。OutputUnit = {' kgmol / m ^ 3 '“K”};

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

NLGR = setpar(NLGR，“名字”,{“体积流速(体积/时间)”.．.% F。“反应堆容积”.．.% V。“指数前非热因子”.．.% k_0。“活化能”.．.% E。玻尔兹曼理想气体常数.．.% R。“反应热”.．.% H。热容乘以密度.．.%高清。总传热系数乘以槽面积});.．.%公顷。NLGR = setpar(NLGR，“单位”,{“m ^ 3 /小时”“m ^ 3”1 / h的“千卡/ kgmol”“千卡/ (kgmol * K)”.．.“千卡/ kgmol”“千卡/ (m ^ 3 * K)”“千卡/ (K * h)”});nlgr.Parameters(1)。Fixed = true;%修复F。nlgr.Parameters(2)。Fixed = true;%固定V。nlgr.Parameters(5)。Fixed = true;修复R。nlgr.Parameters(6)。Fixed = true;%修复H。

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

nlgr.Parameters(1)。最小值= 0;% F。nlgr.Parameters(2)。最小值= 0;% V。nlgr.Parameters(3)。最小值= 0;% k_0。nlgr.Parameters(4)。最小值= 0;% E。nlgr.Parameters(5)。最小值= 0;% R。nlgr.Parameters(6)。最大值= 0;% H。nlgr.Parameters(7)。最小值= 0;%高清。nlgr.Parameters(8)。最小值= 0;%公顷。

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

礼物(nlgr);

nlgr = cstr_m定义的连续时间非线性灰盒模型(MATLAB文件):dx/dt = F(t, u(t)， x(t)， p1，…， p8) 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能够充分激发系统，另一方面必须选择“植物友好型”(化学过程必须保持稳定，测试的持续时间应尽可能短，以便对生产影响最小，等等)。[1]中的文章讨论了CSTR输入信号的选择。有理由认为，多正弦输入u_3比多电平伪随机输入信号更有利。在下面的识别实验中，我们将使用两个这样的输入信号，一个用于估计，一个用于验证，它们是通过上述文章的作者提供的MATLAB®输入数据生成工具(GUI)生成的。

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

负载(fullfile (matlabroot“工具箱”，“识别”，“iddemos”，“数据”，“cstrdata”));Ts = 0.1;每小时10个样品!Ze = iddata(y1, u1, 0.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。ExperimentName =“验证”；

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

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

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

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

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

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

初始CSTR模型的性能

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

CLF比较(ze, nlgr);图;比较(zv nlgr);

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

参数估计

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

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

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

opt.Display =“关闭”；opt.SearchOptions.MaxIterations = 5;NLGR = nlgreyest(ze, NLGR, opt);

估计CSTR模型的性能

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

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

图5:CSTR模型的真实输出与模拟输出的比较。

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

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

礼物(nlgr);

nlgr = cstr_m定义的连续时间非线性灰盒模型(MATLAB文件):dx/dt = F(t, u(t)， x(t)， p1，…， p8) 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_param(broot，“modelworkspace”）;cstrws.assignin (“zv”zv);cstrws.assignin (“nlgr”, nlgr);结束

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

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

IDNLGREY块必须在模拟之前进行配置。这是通过输入持有IDNLGREY模型的MATLAB工作空间变量(这里是nlgr)或通过在标记为“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 = evalin(“基地”，“zsim”）;结束zsim。InputName = nlgr.InputName;zsim。InputUnit = nlgr.InputUnit;zsim。OutputName = nlgr.OutputName;zsim。OutputUnit = nlgr.OutputUnit;zsim。TimeUnit =“小时”；

最后通过绘制Simulink得到的仿真结果来总结实例。金宝app

图(“名字”，估金宝app算CSTR模型的Simulink仿真结果）;为I = 1:zsim。纽约次要情节(zsim。Ny, 1, i);情节(zsim。SamplingInstants ze.OutputData (:, i));标题([“输出#”num2str(我)“:”zsim.OutputName{我}]);包含(”）;轴紧；结束包含([zsim。域“(”zsim。TimeUnit“)”]);