matlab中ode45函數的用法「matlab求解微分方程組例題」

摘要：連續攪拌反應釜（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控制系統進行監控。此控制系統用於工業現場，對提升工作效率具有一定的實際意義。