实验二:结构图法数字仿真
一. 实验目的
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');
四、实验结果