上海电力学院-面向结构图的仿真-实验报告

学生姓名： 学号： 实验题目： 面向结构图的数字仿真 课程名称： 控制系统仿真

实验目的：

? 理解并掌握典型环节的差分方程的计算方法 ? 理解并掌握非线性环节的特性及其仿真

实验内容：

? 应用MATLAB编写一个可以实现面向结构图的仿真程序

? 在以上程序的基础上，实现包含非线性环节的面向结构图的仿真

报告内容：

（1） 给出m文件的程序框图及P118的例题3-2和P141作业3-8的运行结果，并

将结果与simulink 的仿真结果进行比较，并说明编程过程中出现的问题与错误。

（2） 尽量做到界面友好，可以选择不同的非线性环节进行仿真。

此次实验约占整个科目成绩的25%，（其中程序部分占15%，报告部分占10%）提交日期：2013-6-5

实际提交日期：

声明：（包括报告内容和实验程序：此次实验完全自己完成；若得到了同学的帮助，请注明就哪一部分内容请教了同学，说明理由。） 签名： 注：本表打印出来作为封面，空格部分打印出来，手写填写姓名，学号和签名部分，实际提交日期由老师填写。 教师评语及成绩

主程序（非线性）：

clear all clc

P=input('各环节参数P='); WIJ=input('连接阵非零元素WIJ='); n=input('环节个数（系统阶次）n='); Y0=input('阶跃输入幅值Y0='); Tf=input('终止时间Tf=');

L1=input('打印间隔点数（每隔l1输出一次）L1='); T0=input('起始时间T0='); h=input('计算步长h='); X0=input('各环节初始值X0=');

Z=input('非线性标志向量Z=');%有几阶向量列数就为几 S=input('参数向量S=');%与Z同列数 nout=n;%输出环节的编号 Yt0=X0;

%另外加了两个非线性环节

%z=7 线性环节前有继电非线性环节

%z=8 线性环节前有带死区三位继电非线性环节 %z=9 线性环节后有继电非线性环节

%z=10 线性环节后有带死区三位继电非线性环节

%形成闭环各系数阵

% 用“A=diag(P(:,1));B=diag(P(:,2));C=diag(P(:,3));D=diag(P(:,4));”这个语句 %调用时说“Warning: Divide by zero.”得不到理想曲线

%网上说是分母太小，而上诉ABCD中都有0出现，故需要手动给这些点赋值，故有以下语句 A=zeros(1,n);B=zeros(1,n);C=zeros(1,n);D=zeros(1,n);%代替书上语句给ABCD赋值 for i=1:n A(1,i)=P(i,1); B(1,i)=P(i,2); C(1,i)=P(i,3); D(1,i)=P(i,4); end

m=length(WIJ(:,1)); 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

%按环节离散化数字仿真 for i=1:n if(A(i)==0); FI(i)=1;

FIM(i)=h*C(1,i)/B(1,i); FIJ(i)=h*h*C(1,i)/B(1,i)/2; FIC(i)=1; FID(i)=0; if(D(i)~=0);

FID(i)=D(1,i)/B(1,i); else end else

FI(i)=exp( -h*A(1,i)/B(1,i)); FIM(i)=(1-FI(i))*C(i)/A(i); FIJ(i)=h*C(i)/A(i)-FIM(i)*B(i)/A(i); FIC(i)=1;FID(i)=0; if(D(i)~=0);

FIM(i)=(1-FI(i))*D(i)/A(i); FIJ(i)=h*D(i)/A(i)-FIM(i)*B(i)/A(i); FIC(i)=C(i)/D(i)-A(i)/B(i); FID(i)=D(i)/B(i); else end end end

%非线性系统仿真 主程序（满足条件就调用各子函数） Y=zeros(n,1);X=Y;y=0;Uk=zeros(n,1);Ubb=Uk; t=T0:h*L1:Tf; N=length(t); for k=1:N-1 for j=1:L1 Ub=Uk;

Uk=W*Y+W0*Y0; for i=1:n if(Z(i)~=0) if(Z(i)==1)

Uk(i)=satu(Uk(i),S(i)); end if(Z(i)==2)

Uk(i)=dead(Uk(i),S(i)); end if(Z(i)==3)

[Uk(i),Ubb(i)]=backlash(Ubb(i),Uk(i),Ub(i),S(i));

end if(Z(i)==7)

Uk(i)=jidian(Uk(i),S(i)); end if(Z(i)==8)

Uk(i)=deadjidian(Uk(i),S(i)); end end end

Udot=(Uk-Ub)/h; Uf=2*Uk-Ub;

X=FI'.*X+FIM'.*Uk+FIJ'.*Udot; Yb=Y;

Y=FIC'.*X+FID'.*Uf; for i=1:n if(Z(i)~=0) if(Z(i)==4)

Y(i)=satu(Y(i),S(i)); end if(Z(i)==5)

Y(i)=dead(Y(i),S(i)); end if(Z(i)==6)

[Y(i),Ubb(i)]=backlash(Ubb(i),Y(i),Yb(i),S(i)); end if(Z(i)==9)

Y(i)=jidian(Y(i),S(i)); end

if(Z(i)==10)

Y(i)=deadjidian(Y(i),S(i)); end end end end

y=[y,Y(nout)]; end plot(t,y,'g'); grid on hold on

面向结构图（不含非线性）：

clc

%面向结构图仿真（不包含非线性，用四阶龙格库塔） P=input('各环节参数输入P=');

WIJ=input('连接阵非零元素输入WIJ='); n=input('环节个数n='); Yt0=input('各环节初值Yt0='); y0=input('阶跃输入幅值y0='); L1=input('打印间隔点数L1='); T0=input('起始时间T0='); Tf=input('终止时间Tf='); h=input('步长h=') nout=n;%输出环节的编号

A1=diag(P(:,1));B1=diag(P(:,2));C1=diag(P(:,3));D1=diag(P(:,4)); m=length(WIJ(:,1)); 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

Q=B1-D1*W;Qn=inv(Q);R=C1*W-A1;V1=C1*W0;Ab=Qn*R;b1=Qn*V1; Y=Yt0';y=Y(nout);t=T0; N=round((Tf-T0)/(h*L1)); for i=1:N;

for j=1:L1 %四阶龙格库塔 K1=Ab*Y+b1*y0;

K2=Ab*(Y+h*K1/2)+b1*y0; K3=Ab*(Y+h*K2/2)+b1*y0; K4=Ab*(Y+h*K3)+b1*y0; Y=Y+h*(K1+2*K2+2*K3+K4)/6; end

y=[y,Y(nout)]; t=[t,t(i)+h*L1]; end [t',y'] plot(t,y,'g') hold on

子程序就不另加写出来，在电子版中有m文件;