一． 实验目的
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');
四、实验结果