基於Simulink的簡單微分方程組仿真

　　到目前為止，我的所有仿真都是自己敲代碼，一般利用四階龍格庫塔算法、歐拉算法、預估校正算法（分數階）等對系統進行仿真。最近我看了點Simulink的內容，發現很多情況下直接利用Simulink比敲代碼方便得多，但是對於里面很多模塊我不了解，現在對最簡單的微分方程組進行仿真做點筆記，這里所謂的最簡單就是沒有時滯，自治系統，不考慮脈沖、間歇、采樣等因素，就是單純的連續自治微分方程組。
　　當然，這一類系統自己寫代碼或者用MATLAB自帶的ode45函數都比較容易實現，但是，為了學習Simulink，我想選擇這些簡單的系統熟悉一些模塊還是比較好的。利用Simulink搭建微分方程組的方法有很多，我個人感覺，利用自帶的“Matlab Function”模塊是最方便的，只需要將方程輸進去，再加個積分器即可，我們考慮如下Rossler方程：

\(\begin{array}{l} \dot x = - y - z\\ \dot y = x + ay\\ \dot z = b + (x - c)y \end{array}\)
其中，參數\(a=b=0.2\)，\(c=5.7\)，初值可都取零，這個系統實現的Simulink框圖如下：

下面說一下每個模塊的含義（從左到右）：
1“Matlab Function”模塊：將這個模塊拖到空白的simulink搭建界面后，雙擊，即可輸入代碼，本例中，輸入如下代碼：

function y = fcn(u)
y=u;
a=0.2;
b=0.2;
c=5.7;
y(1)=-u(2)-u(3);
y(2)=u(1)+a*u(2);
y(3)=b+(u(1)-c)*u(3);

2 積分器模塊：作用就是積分一次，將這個模塊拖到空白的simulink搭建界面后，雙擊，在Initial condition框里輸入初值，本例[0;0;0]；
3 輸出模塊：作用就是將運行（積分）完的結果輸出到工作空間；
　　這樣一來，咱們的Simulink框圖搭建就結束了，我想這應該是比寫代碼方便很多的，下面的問題是，如何調出工作空間的函數。上面的框圖按照要求搭建完了之后，我們將它保存在你需要的文件夾里，比如命名為“test2.slxc”，后面的擴展名不用管，可能每個版本略有區別。接下來我們要做的就是新建一個m文件來作圖，需要將你新建的m文件和這個搭建的框圖放在一個文件夾，這個務必注意，新建的m文件里代碼如下：

[t,x]=sim('test2',[0,100]);
plot3(x(:,1),x(:,2),x(:,3));

　　這個代碼應該是非常簡單的，而且相信大家應該很容易理解，這里簡單說明一下。首先，第一行，[t,x]是用來接收返回值的，t接收的是時間變量，x則是咱們剛才保存到空間的系統狀態，等號右邊的“sim”函數不用咱們管，函數的第一個參數則是咱們剛保存的simulink框圖文件，第二個參數“[0,100]”的意思是時間范圍，即從0到100，這個可以自由調整。當然，做完圖，咱們可以正常的去加一些圖例、坐標軸的標題等，這和正常的編程是一樣的。咱們運行這里新建的m文件得到下圖：

　　大家是否感覺到這里的圖很不光滑？這與解方程組的算法有關，默認的算法是變步長的，實際上，咱們是可以調整的，一般來說，我還是建議用咱們之前的四階龍格庫塔方法，固定步長（0.01或者0.001等視具體情況而定）。調整算法和步長的步驟如下：首先在咱們的test2這個simulink的搭建界面的菜單欄上找到“simulation”按鈕，點擊，下拉框找到“Model Configuration Parameters”按鈕，點擊，則有如下界面：

　　看到右邊“Solver Selection”選項，點開Type下拉框，選擇“Fixed-Step”，選擇Solver下拉框，選擇“ode4（Runge-Kutta）”，再點開Solver details，設置你需要的固定步長，結果如下圖所示，點擊OK按鈕，然后回到Simulink框圖界面，點擊保存按鈕。

　　之后，再次運行剛才建立的m文件，得到如下結果：

　　不難發現，這次得到的相圖很光滑，到此為止，所有的簡單微分方程組都可以這樣解決了，大家可以自己嘗試其它的單個混沌系統仿真。下面我們考慮一下兩個系統的同步問題，為了方便起見，我們考慮兩個神經網絡系統同步，系統的模型如下：
\(\begin{array}{l} \dot x = Ax + Bf(x)\\ \dot y = Ay + Bf(y) + u\\ u = - D(y - x) \end{array}\)
其中，\(f(x) = (|x + 1| - |x - 1|)/2\)，\(A = - {I_3}\)，\(D=5{I_3}\)，\(B = \left[ {\begin{array}{*{20}{c}} {1.25}&{ - 3.2}&{ - 3.2}\\ { - 3.2}&{1.1}&{ - 4.4}\\ { - 3.2}&{4.4}&1 \end{array}} \right]\)，初值：\(x(0) = [0.1;0.1;0.1]\)，\(y(0) = [0.1;0.1;0.100001].\)如果不控制，即沒有反饋項，盡管兩個系統的初值很近，但是隨着時間的推移，狀態圖如下：

　　這是混沌系統的特性，對初值極為敏感。下面取控制參數如前文所示，得到如下結果：

　　在這里處理同步問題時，我們也是將其視為多個微分方程構成的一個大的微分方程組，大家應該能看到，每個程序的區別就是“Matlab Function”模塊、積分器中的初值和最后的m文件，這里一共有六個方程，總體思路是將前面三個看成第一個系統，后面三個看成相應系統，由此可得積分器中的初值輸入為“[0.1;0.1;0.1;0.1;0.1;0.100001]”， “Matlab Function”模塊的代碼如下，請大家好好理解一下：

function y = fcn(u)
y=u;
A=-eye(3);D=5*eye(3);
B=[1.25 -3.2 -3.2;-3.2 1.1 -4.4;-3.2 4.4 1];
x1=u(1:3);
y1=u(4:6);
y=[A*x1+0.5*B*(abs(x1+1)-abs(x1-1));A*y1+0.5*B*(abs(y1+1)-abs(y1-1))-D*(y1-x1)];

最后，咱們自己寫的m文件代碼如下：

[t,xx]=sim('test2',[0,150]);
x=xx(:,1:3);
y=xx(:,4:6);
subplot(1,3,1)
plot(t,x(:,1),'b',t,y(:,1),'r');
legend('x_{1}','y_{1}');
subplot(1,3,2)
plot(t,x(:,2),'b',t,y(:,2),'r');
legend('x_{2}','y_{2}');
subplot(1,3,3)
plot(t,x(:,3),'b',t,y(:,3),'r');
legend('x_{3}','y_{3}');

　　我不對這里的ｍ文件中的代碼做任何解釋，請大家自行理解。在這里，我們想做驅動系統或者相應系統的相圖應該也是輕而易舉的。另外，自適應控制方法也可以很容易實現，包括復雜網絡的同步等都可以實現。總的來說，只要不含時滯、非自治、脈沖、間歇等的簡單微分方程組，利用這個框圖都很容易實現，只是“Matlab Function”模塊、積分器中的初值和最后的m文件不同。這里省去了復雜的代碼，只需要將方程組右端的函數輸入即可，麻煩的地方就是最后的m文件，但是，無論用什么方法，都無法省略最后的m文件中的相關語句。簡單來說，利用simulink框圖可以得到系統各個時刻的狀態，最后的m文件就是將這些狀態提取出來，作圖顯示。

2019.10.19 王飛