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

二、實驗內容

　　1、實驗要求

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

　　模擬步長：0.001，模擬時間：1.5s。

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

　　（1）環節連線關係：

　　即：

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

　　（2）輸出方程：

　　其中： 、

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

　　3、系統引數輸入與閉環引數的生成：
　　（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演算法求解所得的閉環狀態方程。

三、實驗完整程式碼

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');

四、實驗結果

&n