matlab求解微分方程组例题(matlab中ode45函数的用法)

matlab求解微分方程组例题(matlab中ode45函数的用法)

  摘要:连续搅拌反应釜(CSTR)在生产过程中得到了广泛应用。因其在实际生产过程中会受到许多不利因素的影响,不易实现面向性能的控制。以连续搅拌反应釜为对象,采用常规PID控制,为了达到实时修改模型参数,动态显示控制曲线和变量数值的目的,设计了GUI人机界面,使得用户可以方便、实时地对CSTR控制系统进行监控。

  0引言

  连续搅拌反应釜(Continuous Stirred Tank Reactor, CSTR)作为一类化学反应器,由于其成本低、热交换能力强和产品质量稳定等特点,成为生产聚合物的核心设备,在化工、发酵、生物制药、石油生产等工业生产过程中得到了广泛的应用[1]。其示意图如图1所示。

  在连续搅拌反应釜系统中通过控制其内部的工艺参数(如温度、浓度等)的稳定,来保证反应的正常进行,其控制质量直接影响到生产的效益和质量指标。连续搅拌反应釜的对象是高度非线性的化学反应系统,一般用一组非线性常微分方程来描述该反应釜的数学模型[2]。本实验通过采用单回路控制系统,在回路中运用基于四阶五级RungeKutta的PID控制算法,通过整定PID参数实现对被控变量的控制。同时在控制系统中设置扰动,模拟现场控制中所遇到的干扰。为了能够实时监控工艺参数,本设计采用GUI人机界面,动态显示被控对象及状态变量的控制结果。

  1CSTR模型

  实验针对的化学反应是由环戊二烯(组分A)生成主产品环戊烯(组分B)和副产品二环戊二烯(组分D)以及由环戊烯继续反应生成的副产品环戊酮(组分C)。化学反应方程如下:

  Ak1Bk2C,2Ak3D。基于能量守恒定律,该CSTR模型的理想动态特性可以由以下非线性微分方程组描述:

  如上各参数中,CA为反应器中物质A的浓度,CB反应器中物质B的浓度,CA0为A的进料浓度,K1、K2、K3为3个化学反应的反应速率,V为物料A进料体积流量,VR为反应器体积,T为反应器温度,Tk为冷却剂温度,T0为反应器入口温度,ΔHRAB、ΔHRBC、ΔHRAD分别为K1、K2、K3反应放出的热量,ρ为反应器液体密度,Cρ为反应器液体热容,Kw为冷却套的传热系数,AR冷却套传热面积,Ei为第i个反应的反应激活能量。

  该CSTR模型的常微分方程组由3个微分方程组成,即将CA、CB、T作为系统3个状态变量建立微分方程,取冷却剂温度Tk为控制系统的操作变量,反应器中物质浓度CA作为被控变量,反应器入口温度T0和浓度CA0、物料进料体积流量V是波动的,可以作为外部的扰动。各参数的值如表1所示。表1CSTR模型常微分方程组参数表变量名变量符号参数值单位物质A进料体积流量V14.19L/h反应器入口温度T079.7℃物料A进料初始浓度CA05.1mol/L反应放出的热量ΔHRAB-4.2KJ/mol反应放出的热量ΔHRBC11KJ/mol反应放出的热量ΔHRAD41.85KJ/mol反应器中液体的密度ρ0.934 2kg/L反应器液体热容Cρ3.01KJ/(kg·K)冷却套的传热系数kw4 032KJ/(h·m2·K)冷却套传热面积AR0.215m2反应器体积VR10L反应速率系数k101.287×1012h-1反应速率系数k201.287×1012h-1反应速率系数k309.043 2×109h-1反应的反应激活能量E1-9 758.3K反应的反应激活能量E2-9 758.3K反应的反应激活能量E3-8 560K

  根据CSTR模型的微分方程,以反应器温度T、反应器中物质A的浓度CA、反应器中物质B的浓度CB三者为状态变量,以冷却剂Tk为控制变量,建立关于微分方程的M文件。

  2CSTR过程仿真控制研究

  2.1控制算法

  本实验采用基于四阶五级RungeKutta的PID控制算法。四阶五级RungeKutta算法是一种求解微分方程近似解的数值方法,实际上是间接使用泰勒级数法的一种计算方法。该算法精度高,能对误差进行抑制。在区间[k,k+d]上用y(k)的值来估算或预测y(k+d)的值,得到预测值(k+d),将此预测值作为反馈信号与期望设定值进行比较得出偏差,作为PID控制的输入,依照PID控制律来设定控制器的输出,完成对被控对象的控制。

  MATLAB中的ode函数专门用于求解微分方程,而ode45表示采用四阶五级Runge-Kutta算法,它用四阶方法提供候选解,五阶方法控制误差,是一种自适应步长(变步长)的常微分方程数值解法,本实验M文件中就采用ode45求解微分方程。

  首先确定仿真的采样时间、起止时间以及每一步模型仿真的时间区间,并为龙格库塔算法设定初始值。

  然后初始化PID控制器,并设定PID参数和设定值。经过PID参数的调整,得到Kc=0.03;Ti=4;Td=0.05。

  最后运用循环语句的形式编写基于四阶五级RungeKutta法的PID控制算法。

  2.2添加扰动

  为模拟真实现场控制系统环境,这里需要添加3个扰动:物质A进料流量(QIn)扰动,反应器入口温度扰动(To),物质A进料浓度扰动(Ca0)。

  2.3控制结果

  将初始值设为y0=[0;1;80.7],在控制系统的作用下最终达到稳定,如图2所示。

  上排从左至右分别表示反应器物质A浓度(被控变量),反应器物质B的浓度和反应器温度,即3个状态变量。下排从左至右分别表示冷却剂的温度(操纵变量)和控制误差(即控制器的输入)。可发现由于加了积分作用,控制系统的余差为0,并且控制效果较好。

  3GUI界面的设计

  为了能够实时改变控制系统模型的参数,在程序运行过程中添加扰动,并使被控对象及状态变量的控制结果动态显示,需要添加一个GUI界面来实现这些功能。

  首先在命令窗口中键入guide,GUIDE实际上是一套MATLAB工具箱[3]。启动GUIDE后,会出现GUIDE Quick Start 对话框,选择新建一个GUI界面,这里选用GUI with Axes and Menu模板,点击OK后进入版面设计窗口。利用窗口左侧的工具箱可以选择需要添加的组件,用来输入扰动和采样时间,同时显示控制系统变量的数值和曲线。完成版面设计后,以CSTR_GUI文件名保存,此时设计内容会保存在两个文件中,一个是FIG文件,一个是M文件。GUI界面设计如图3所示。

  然后打开M文件CSTR_GUI.M,对GUI进行编程。设置并启动0.2 s定时器,

  另外,需要编写定时器中断响应函数TimerCallback,里面包含对界面扰动参数和采样时间的读取,同时实现对GUI人机界面中变量值的更新。

  4CSTR模型控制结果

  设定初始值:反应器中物质A的浓度CA=2.14 mol/L,反应器中B的浓度CB=1.05 mol/L,反应器温度T=80.7 ℃,即y0=[2.14;1.05;80.7]。

  在无扰动的情况下,控制结果如图4。

  此时PID参数为Kc=0.03,Ti=4 s,Td=0.05 s,从反应器中物料A浓度曲线来看,控制作用响应速度快,超调小,且没有稳态误差,控制效果较好。由于CSTR模型控制系统是完全用MATLAB进行仿真的,曲线平稳之后没有出现任何波动,这与实际现场控制状况不太相同。

  在控制系统达到稳态之后,可以在GUI界面上对某一参加数上一扰动(例如增加5%的物质A进料流量扰动),观察控制系统的调节能力,如图5所示。

  从图5可以看出在加入进料流量扰动后,反应器中A浓度曲线出现一定程度的波动,这一波动的大小受扰动值大小的影响,但物质A的浓度很快又恢复到设定值,说明控制系统有较强的调节能力。另外还能发现当进料流量增加5%后,反应器中物质B的浓度在重新达到稳态后并没有回到原值(从1.046 8 mol/L增加到1.047 8 mol/L),同时反应器的温度也有所增加,因此在通过改变进料流量改变物质B浓度的时候,需要注意反应器的温度,要避免触碰到反应器温度的上下限。

  5结论

  本文基于MATLAB建立CSTR对象模型,依据现实的生产环境以及各种干扰因素,通过整定PID参数完成对被控变量的控制,取得良好的控制效果。同时,为了达到实时修改模型参数、动态显示控制曲线和变量数值的目的,引入了GUI人机界面,使得用户可以方便、实时地对CSTR控制系统进行监控。此控制系统用于工业现场,对提升工作效率具有一定的实际意义。

推荐阅读