【控制系統數字仿真與CAD】實驗二：結構圖法數字仿真

一． 實驗目的
　　1. 掌握結構圖法仿真複雜控制系統的方法；
　　2. 掌握複雜系統聯接矩陣W和輸入聯接矩陣W0的求解過程；
　　3. 掌握複雜系統的環節鏈接,矩陣A、 B、 C、D的求解過程；
　　4. 掌握MATLAB編制結構圖數字仿真的程序，並運用MATLAB程序對仿真結構進行處理。
　　

　　  

2、實驗內容

　　一、實驗要求

　　上圖系統中，各子系統的初始值均爲零，輸入Un*爲單位階躍。試編寫仿真程序，求各環節輸出響應（求解閉環狀態方程時，指定採用RK4方法）。記錄仿真結果及程序中計算獲得的Ab和Bb。

　　仿真步長：0.001，仿真時間：1.5s。

　　二、根據結構圖寫出聯接矩陣

　　（1）環節鏈接關係：

　　

　　即：

　　其中爲聯接矩陣， 爲輸入聯接矩陣。爲各環節輸入向量，爲各環節輸出向量 。

　　（2）輸出方程：

　　

　　其中：、

　　（3）環節輸入輸出關係：(A+Bs)y(s)=(C+Ds)u(s)

　　三、系統參數輸入與閉環參數的生成：
　　（1）輸入輸出關係參數輸入矩陣P：

P=[ 0, tn,    kn,  kn*tn ;
    0, ti,    ki,  ki*ti ;
    1, ts,    ks,  0     ;
    1, tl,    1/R, 0     ;
    0, ce*tm, R,   0     ];

　　（2）閉環係數陣A， B， C， D的生成：

A=diag( P(:,1) );
B=diag( P(:,2) );
C=diag( P(:,3) );
D=diag( P(:,4) );

　　（3）系統參數輸入與閉環參數的生成：

　　鏈接矩陣輸入方法：創建非零元素矩陣WIJ（m*3型），將非零元素按照i， j， wij次序輸入。其中i爲被做用環節號； j爲做用環節號； wij爲做用關係值； m爲非零元素個數（包括W0陣）。

　　 

WIJ=[1, 0, 1     ;
     1, 5, -Alpha;
     2, 1, 1     ;
     2, 4, -Beta ;
     3, 2, 1     ;
     4, 3, 1     ;
     4, 5, -ce   ;
     5, 4, 1     ];

　　（4）閉環係數陣W， W0的生成：

m=length(WIJ(:,3));
n=5;
W0=zeros(n,1);
W =zeros(n,n);
for k=1:m
    if(WIJ(k,2)==0)
        W0(WIJ(k,1))=WIJ(k,3);
    else
        W( WIJ(k,1),WIJ(k,2) ) = WIJ(k,3);
    end
end

　　（5）系統參數輸入與閉環參數的生成：
　　閉環係數陣Ab， Bb的生成：

Q=B-D*W;    Qn=inv(Q);
R=C*W-A;    V1=C*W0;
Ab=Qn*R;    Bb=Qn*V1;

　　（6）數值積分求解

　　採用RK4算法求解所得的閉環狀態方程。

3、實驗完整代碼

clear;
clc;
%*******  初始化參數  ********%
kn=26.7;        tn=0.03;	ki=0.269;	ti=0.067;
ks=76;          ts=0.00167;
R=6.58;         tl=0.018;   tm=0.25;	ce=0.131;
Alpha=0.00337;	Beta=0.4;	Idl=0;
%*******  輸入矩陣  ********%
P=[ 0, tn,    kn,  kn*tn ;
    0, ti,    ki,  ki*ti ;
    1, ts,    ks,  0     ;
    1, tl,    1/R, 0     ;
    0, ce*tm, R,   0     ];
%*******  閉環係數陣  ********%
A=diag( P(:,1) );
B=diag( P(:,2) );
C=diag( P(:,3) );
D=diag( P(:,4) );
WIJ=[1, 0, 1     ;
     1, 5, -Alpha;
     2, 1, 1     ;
     2, 4, -Beta ;
     3, 2, 1     ;
     4, 3, 1     ;
     4, 5, -ce   ;
     5, 4, 1     ];
%*******  計算W、W0  ********%
m=length(WIJ(:,3));
n=5;
W0=zeros(n,1);
W =zeros(n,n);
for k=1:m
    if(WIJ(k,2)==0)
        W0(WIJ(k,1))=WIJ(k,3);
    else
        W( WIJ(k,1),WIJ(k,2) ) = WIJ(k,3);
    end
end
%*******  閉環系統陣Ab、Bb生成  ********%
Q=B-D*W;    Qn=inv(Q);
R=C*W-A;    V1=C*W0;
Ab=Qn*R;    Bb=Qn*V1;
h  = 0.001;
x  = [];
y  = [0;0;0;0;0];
y1 = [];
y2 = [];
y3 = [];
y4 = [];
y5 = [];
% ********  4階龍格-庫塔法  **********%
for i = 0:1:1000
    t = i*h;     % simulation time: 1s
    x(i+1) = t;
    k1 = Ab *  y         + Bb;
    k2 = Ab * (y+h*k1/2) + Bb;
    k3 = Ab * (y+h*k2/2) + Bb;
    k4 = Ab * (y+h*k3)   + Bb;
    y = y + (k1 + 2*k2 +2*k3 + k4) * h/6;
    y1(i+1) = y(1,1);
    y2(i+1) = y(2,1);
    y3(i+1) = y(3,1);
    y4(i+1) = y(4,1);
    y5(i+1) = y(5,1);
end
%*******  繪圖  ********%
plot(x, y1, x, y2,'--', x, y3,'*', x, y4,'.-', x, y5,'x');
xlabel(' time  /s ');ylabel(' Output ');title('系統輸出狀況 ');
legend('Y1','Y2','Y3','Y4','Y5');

4、實驗結果