自控原理实验

实验一面向微分方程的数值积分法仿真一、实验目的1、掌握数值积分法的基本概念、原理及应用；2、用龙格-库塔法解算微分方程，增加编写仿真程序的能力；3、分析数值积分算法的计算步长与计算精度、速度、稳定性的关系；4、对数值算法中的“病态问题”进行研究。二、实验原理1、数值积分法的基本概念及定义用被积函数的有限个抽样值的离散和或加权平均值近似地代替定积分的值。在求函数/(X)的定积分- -■- ■-■■■■■■ ■ 时，常常无法用初等函数表示原函数，因此能按牛顿-莱布尼茨公式匸冷血”⑻一能)计算积分值的定积分是不多的。另外，当/(x)是列表函数时，也不能使用式计算它的积分值。上述事实说明，必须研究近似估算积分的数值积分方法。2、龙格-库塔法：龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施误差进行抑制，所以其实现原理也较复杂。该算法是构建数学支持的基础之上的。对于一阶精度的欧拉公示有：y(i+1)=y(i)+h*K1，K1=f(xi，yi)，当用点xi处的斜率近似值K1与右端点xi+1处的斜率K2的算术平均值作为平均斜率K*的近似值，那么就会得到二阶精度的改进欧拉公式y(i+1)=(y(i)+h)*(K1+K2)/2，K1=f(xi，yi)，K2=f(x(i)+h),y(i+h)*K1)依次类推，如果在区间［xi,xi+l］内多预估几个点上的斜率值KI、K2、……Km,并用他们的加权平均数作为平均斜率K*的近似值，显然能构造出具有很高精度的高阶计算公式。经数学推导、求解，可以得出四阶龙格－库塔公式，也就是在工程中应用广泛的经典龙格－库塔算法：y(i+l)=y(i)+h*(K1+2*K2+2*K3+K4)/6,KI二f(x(i),y(i)),K2二f(x(i)+h/2,y(i)+h*K1/2)K3=f(x(i)+h/2,y(i)+h*K2/2)K4=f(x(i)+h,y(i)+h*K3)通常所说的龙格-库塔法是指四阶而言的，我们可以仿二阶、三阶的情形推导出常用的标准四阶龙格-库塔公式。三、实验内容1、 已知系统微分方程及初值条件》=£+”丁(0)=1取步长h=0.1，试分别用欧拉方程法和RK4法求t=2h时的y值，并将求得的值与解析解列°二牯T-1比较(将三个解绘于同一坐标中，且用数值进行比较)，说明造成差异的原因。(①编程完成；②选用MATLABode函数完成。)③、= 4C■&2、已知系统的传递函数为 在单位阶跃输入下，系统响应的解析解为皿=1刚-必站陶-卩严’-山严’试分别用欧拉方程法和RK4法对系统进行仿真(编程完成)：比较两种数值积分解与解析解的逼近程度；(绘图)改变步长，分析步长对数值解精度的影响；不断加大步长，分析计算稳定性的变化；3、 求下图所示系统的阶跃响应y(t)的数值解。(v=1,k=1,开始时间10=0,结束时间tf=10,h=0.25)分析k、v对系统响应的影响。(①编程用RK4求解；②Simulink)(修正：G(s)表达式的分母中，0.25s+1改为(0.25s+1)的二次方)

4、已知系统状态状态方程为-2.1 19i=19 -2140'-40-20-2.1 19i=19 -2140'-40药0扎忑初二［10-if40（修正：x的系数矩阵中最后一个元素为-40,而非40）采用RK4法，步长分别取h=0.01、0.04、0.06,求系统的零输入响应，并绘图分析各状态变量的响应状态及产生的原因。（提示：病态系统）四、实验步骤1、打开电脑选择桌面上matlab。若桌面上没有则单机开始，选择程序，单击matlab程序，打开matlab程序。2、 新建一个M文件，并命名为“微分方程”。3、 按照前面的分析结果，分析系统整体的性能，进行编程，将编程的结果输入用户界面。4、 点击运行。5、 对程序代码中出现的错误进行改正，并重新运行。6、 观察结果，记录数据，进行参数变化，记录变化的结果。7、 将实验结果记录至之前的表格中。8、截图，按照要求完成实验报告。五、实验代码及结果1.t3=0t0=0;tf=2;h=0.1;y1=1;y2=1;y3=1;t1=0;t2=0;n=round(tf-t0)/h;t3=0fori=1:ny1(i+1)=y1(i)+h*(2*h+y1(i));t1=[t1,t1(i)+h];endfori=1:nk1=y2(i)+t2(i);k2=y2(i)+h*k1/2+t2(i)+h/2;k3=y2(i)+h*k2/2+t2(i)+h/2;k4=y2(i)+h*k3+t2(i)+h;y2(i+1)=y2(i)+h*(k1+2*k2+2*k3+k4)/6;t2=[t2,t2(i)+h];endfori=1:ny3(i+1)=2*exp(t3(i))-t3(i)-1;t3=[t3,t3(i)+h];endplot(t1,y1,'r',t2,y2,'g',t3,y3,'b')2.a=[010;001;-22.06-27-10];b=[0;0;1];c=[40.600];X1=[0;0;0];Y1=0;u=1;Y2=0;Y3=0;X2=[0;0;0];h=0.1;t0=0;tf=2;t1=0;t2=0;t3=0;N=round(tf-t0)/h;fori=1:Nk1=a*X1+b;k2=b+a*(h*k1/2+X1);k3=b+a*(h*k2/2+X1);k4=b+a*(h*k3+X1);X1=X1+h*(k1+2*k2+2*k3+k4)/6;Y1=[Y1,c*X1];t1=[t1,t1(i)+h];endfori=1:Nx=X2(:,i)+h*(a*X2(:,i)+b*u);y=c*x;X2=[X2,x];Y2=[Y2,y];t2=[t2,t2(i)+h];endfori=1:NY3(i+1)=1.84-4.95*t3(i)*exp(-1.88*t3(i))-1.5*exp(-1.88*t3(i))-0.34*exp(-6.24*t3(i))t3=[t3,t3(i)+h];endplot(t1,Y1,'r',t2,Y2,'g',t3,Y3,'b')3.k=1;a=conv([100],conv([0.251],[0.251]))b=[2*kk]X0=[0000]v=1;n0=4;tf=10;h0=0.25;r=1;V=v;n=n0;T0=0;Tf=tf;h=h0;R=r;b=b/a(1);a=a/a(1);A=a(2:n+1);A=[rot90(rot90(eye(n-1,n)));-fliplr(A)];B=[zeros(1,n-1),1]';m1=length(b);C=[fliplr(b),zeros(1,n-m1)];Ab=A-B*C*V;X=X0';y=0;t=T0;N=round(Tf-T0)/h;fori=1:Nk1=Ab*X+B*R;k2=Ab*(X+h*k1/2)+B*R;k3=Ab*(X+h*k2/2)+B*R;k4=Ab*(X+h*k3)+B*R;X=X+h*(k1+2*k2+2*k3+k4)/6;y=[y,C*X];t=[t,t(i)+h];end[t',y']plot(t,y)4.A=[-2119-20;19-2120;40-40-40];x=[1;0;-1];X=x;t=0;t0=0;tf=2;h=0.01;n=round(tf-t0)/h;fori=1:nk1=A*x;k2=A*(x+h*k1/2);k3=A*(x+h*k2/2);k4=A*(x+h*k3);x=x+h*(k1+2*k2+2*k3+k4)/6;X=[X,x];t=[t,t(i)+h];endy1=X(1,:);y2=X(2,:);y3=X(3,:);plot(t,y1,'r',t,y2,'g',t,y3,'b')六、预习要求及思考题要求实验前，必须预习实验知识点，按实验内容的要求的确定仿真方案，完成程序设计，以便在实验时进行调试。分析思考题：1、在进行仿真计算时，是否选用的数值积分法的阶次越高越好？答：不是，在进行仿真计算时。若选取的数值积分阶数过高，则会影响系统的稳定性，在选取的时候要综合各个方面进行考虑。2、选用数值积分法进行仿真的原则答：1、精度仿真结果的精度主要受三项误差影响，即截断误差、舍入误差、累计误差。2、计算速度计算速度取决于所用数值方法和计算步长。3、稳定性数值稳定性主要与计算步长h有关，而且与仿真对象的时间常数有关。七实验小节本次试验是一次综合性的实验，首先巩固了Matlab/Simuli仿真环境，然后掌握Simulin图形化建模方法，最后要完成的是“双闭环直流电动机调速系统”的性能分析，前面过程还比较顺利，到了K值的确定这个点就出了问题，最后在老师的帮助下才理解了K是一个定值，在前面的说明中提到过。最后我还是磕磕绊绊的做出了结果，可惜的是后面的“电流环跟随性能仿真实验”，以及“转速环抗扰性能仿真”实验由于时间关系没能来及得及完成，但是听老师说这两个实验还是比较简单的，按步就班的方法来操作就行了。这门课程的实验给我的感觉就是实用性很强，在实际工程中遇到的各种问题怎样通过数学方法来建模？以及建模后用什么工具来实现仿真等等，这些专业性知识都是我们学自动化必须要掌握的。控制系统数字仿真与CAD实验二第一部分面向结构图的数值积分法仿真一、 实验目的1、 加深理解连续系统面向结构图仿真的原理及特点。2、 进一步掌握数字积分法解算微分方程的方法。3、 增加编写仿真程序的能力。二、 实验原理1、 连续系统：时间和各个组成部分的变量都具有连续变化形式的系统。连续系统的运动过程常可用一个或一组反映其变量间的相互联系和因果关系的微分方程来描述，通常称这一方程为它的数学模型。当微分方程的系数为常数时称为定常系统，当系数随时间而变化时则称为时变系统。如果连续系统的数学模型具有常微分方程的形式，则该系统属于集中参数系统;而当具有偏微分方程的形式时,相应的连续系统属于分布参数系统。对自动控制系统，只有当受控过程和控制方式同为连续时的系统才称为连续控制系统。。2、 数字积分法:数字积分法又称微分分析法(DigitaldifferentialAnalyzer),是在数字积分器的基础上建立起来的一种插补算法。数字积分法的有点是，易于诗选多坐标联动，叫容易地实现二次曲线、高次曲线的插补，并具有运算速度快,应用广泛的特点。三、 实验内容1、用面向方框图的数字仿真方法对下列系统进行仿真。2、求解下图所示系统在f=-1(t)阶跃扰动作用下第④、第⑤环节的动态过程。分别用面向框图的数值积分法(RK4法)、MATLAB中有关系统建模的命令和Simulink三种方法求解。四、 实验步骤1、 打开电脑选择桌面上matlab。若桌面上没有则单机开始，选择程序，单击matlab程序，打开matlab程序。2、 新建一个M文件，并命名为“微分方程”。3、 按照前面的分析结果，分析系统整体的性能，进行编程，将编程的结果输入用户界面。4、 点击运行。5、 对程序代码中出现的错误进行改正，并重新运行。6、 观察结果，记录数据，进行参数变化，记录变化的结果。7、 将实验结果记录至之前的表格中。8、 截图，按照要求完成实验报告。五、 实验代码及结果P=[0,0.07,1,0.14;1,0.012,1,0;0,0.05,1,0.15;10,1,1,0;1,0.01,0,0.0008];WIJ=[1,0,1;1,4,-1;2,1,1;3,2,1;3,5,-1;4,3,1;5,4,1]n=5;Y0=1;Yt0=[00000];h=0.01;t=0;L1=5;T0=0;Tf=2;nout=4;A=diag(P(:,1));B=diag(P(:,2));C=diag(P(:,3));D=diag(P(:,4));m=length(WIJ(:,1));W0=zeros(n,1);W=zeros(n,n);fork=1:mif(WIJ(k,2)==0)W0(WIJ(k,1))=WIJ(k,3);elseW(WIJ(k,1),WIJ(k,2))=WIJ(k,3);end;end;Q=B-D*W;Qn=inv(Q);R=C*W-A;V1=C*W0;Ab=Qn*R;b1=Qn*V1;Y=Yt0';y=Y(nout);t=T0;N=round((Tf-T0)/(h*L1));fori=1:N;forj=1:L1;k1=Ab*Y+b1*1;k2=Ab*(Y+h*k1/2)+b1*1;k3=Ab*(Y+h*k2/2)+b1*1;k4=Ab*(Y+h*k3)+b1*1;Y=Y+h*(k1+2*k2+2*k3+k4)/6;end;y=[y,Y(nout)];t=[t,t(i)+L1*h];end;plot(t,y)2.P=[10.014910;10.0025426.66670;10.39100.4199;00.24810;12.8810];WIJ=[151;211;341;421;43-1;501;54-1;];n=5;Y0=-1;Yt0=[00000];h=0.005;L1=10;T0=0;Tf=10;nout=5;A=diag(P(:,1));B=diag(P(:,2));C=diag(P(:,3));D=diag(P(:,4));m=length(WIJ(:,1));W0=zeros(n,1);W=zeros(n,n);fork=1:mif(WIJ(k,2)==0);W0(WIJ(k,1))=WIJ(k,3);elseW(WIJ(k,1),WIJ(k,2))=WIJ(k,3);end;end;Q=B-D*W;Qn=inv(Q);R=C*W-A;V1=C*W0;Ab=Qn*R;b1=Qn*V1;Y=Yt0';y=Y(nout);t=T0;N=round((Tf-T0)/(h*L1));fori=1:N;forj=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)gridon;xlabel('t/s');ylabel('Amplitude');axis([010-0.080.02]);六、实验小节在做向结构图的数值积分法仿真的实验前，我以为不会难做，就像以前做实验一样做完实验，然后两下子就将实验报告做完。直到做完测试实验时，我才知道其实并不容易做，但学到的知识与难度成正比,使我受益匪浅.。在做实验前一定要将课本上的知识吃透，因为这是做实验的基础，否则，在老师讲解时就会听不懂，这将使你在做实验时的难度加大，浪费做实验的宝贵时间•比如做应变片的实验，你要清楚电桥的各种接法，如果你不清楚，在做实验时才去摸索，这将使你极大地浪费时间,使你事倍功半•做实验时，一定要亲力亲为，务必要将每个步骤，每个细节弄清楚，弄明白,实验后，还要复习,思考这样，你的印象才深刻，记得才牢固。否则，过后不久你就会忘得一干二净，这还不如不做•做实验时，老师还会根据自己的亲身体会，将一些课本没有的知识教给我们拓宽我们的眼界，使我们认识到这门课程在生活中的应用是那么的广泛•。通过这次向结构图的数值积分法仿真的实验，使我学到了不少实用的知识,更重要的是，做实验的过程，思考问题的方法，这与做其他的实验是通用的真正使我们受益匪浅。第二部分面向结构图的离散相似法仿真一、 实验目的1、 掌握离散相似法的基本概念和原理；2、 典型环节离散系数的求取及差分方程表示；3、 掌握非线性系统的数字仿真方法二、 实验原理离散相似法：离散相似法是指将一个连续系统进行离散化处理，然后求得与它等价的离散模型（差分方程）的方法。主要应用于连续系统建模与仿真领域中。从连续系统离散化的角度出发，建立连续系统模型的等价离散化模型，并用采样系统的理论和方法介绍另一种常用的仿真算法。这种算法使得连续系统在进行（虚拟的）离散化处理后仍保持与原系统“相似”，故称之为离散相似算法。三、 实验内容1、已知控制系统结构图如图所示，设输入阶跃函数幅值Y0=10,滞环非线性参数s=1（滞环宽度），请用离散相似法编程和Simulink法对系统进行如下分析：1） 不考虑非线性环节影响时，求解y（t）的阶跃响应；2） 考虑非线性环节影响，其余参数不变，求解y（t）并与线性情况所得结果进行比较；3）改变的滞环非线性参数s,分析该非线性对系统的影响。2、系统结构图如图所示，先理论分析该系统是否会产生自激振荡，若会求出振荡的振幅和频率，并用离散相似法编程和Simulink法验证分析的结果。四、实验步骤1、 打开电脑选择桌面上matlab。若桌面上没有则单机开始，选择程序，单击matlab程序，打开matlab程序。2、 新建一个M文件，并命名为“微分方程”。3、 按照前面的分析结果，分析系统整体的性能，进行编程，将编程的结果输入用户界面。4、 点击运行。5、 对程序代码中出现的错误进行改正，并重新运行。6、 观察结果，记录数据，进行参数变化，记录变化的结果。7、将实验结果记录至之前的表格中。8、截图，按照要求完成实验报告。五、实验代码及结果

3.P=[110510;10.513.P=[110510;10.51WIJ=[101;211;32X0=[0000];Z=[0006];S=[000h=0.01;L1=25;n=4;0;10.110;0110];1;431;14-1];1];T0=0;Tf=10;nout=4;Y0=10;A=(P(:,1));B=(P(:,2));C=(P(:,3));D=(P(:,4));m=length(WIJ(:,1));W0=zeros(n,1);W=zeros(n,n);fork=1:mif(WIJ(k,2)==0);W0(WIJ(k,1))=WIJ(k,3);elseW(WIJ(k,1),WIJ(k,2))=WIJ(k,3);endendfori=1:nif(A(i)==0);I(i)=1;FIM(i)=h*C(i)/B(i);FIJ(i)=h*h*C(i)/B(i)/2;FIC(i)=1;FID(i)=0;if(D(i)~=0);FID(i)=D(i)/B(i);elseendelseFI(i)=exp(-h*A(i)/B(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=(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);elseendendendY=zeros(n,1);X=Y;y=0;Uk=zeros(n,1);Ubb=Uk;t=T0:h*L1:Tf;N=length(t);fork=1:N-1forl=1:L1Ub=Uk;Uk=W*Y+W0*Y0;fori=1:nif(Z(i)~=0)if(Z(i)==1)Uk(i)=satu(Uk(i),S(i));endif(Z(i)==2)Uk(i)=dead(Uk(i),S(i))endif(Z(i)==3)[Uk(i),Ubb(i)]=backlash(Ubb(i),Uk(i),Ub(i),S(i));EndendendUdot=(Uk-Ub)/h;Uf=2*Uk-Ub;X=FI'.*X+FIM'.*Uk+FIJ'.*Udot;Yb=Y;Y=FIC'.*X+FID'.*Uf;fori=1:nif(Z(i)~=0)if(Z(i)==4)Y(i)=satu(Y(i),S(i));endif(Z(i)==5)Y(i)=dead(Y(i),S(i));endif(Z(i)==6)[Y(i),Ubb(i)]=backlash(Ubb(i),Y(i),Yb(i),S(i));endendendendy=[y,Y(nout)];endplot(t,y)4.P=[0,1,1,0;1,1,1,0;2,1,1,0];WIJ=[101;13-1;211;321];n=3;Y0=10;X0=[0000];Yt0=[0000];h=0.01;L1=25;t0=0;tf=10;nout=3;Z=[006];S=[001];A=(P(:,1));B=(P(:,2));C=(P(:,3));D=(P(:,4));m=length(WIJ(:,1));W0=zeros(n,1);W=zeros(n,n);fork=1:mif(WIJ(k,2)==0);W0(WIJ(k,1))=WIJ(k,3);elseW(WIJ(k,1),WIJ(k,2))=WIJ(k,3);end;end;sp3_3sp3_4plot(t,y)六、实验小节本次实验主要是通过改变离散步长，对结构图的各个环节逐个刷新进行仿真，得到仿真曲线与simulink的相对比，分析系统的稳定性，了解离散相似法的原理，熟悉仿真程序的结构。实验过程中，遇到的最大困难是对matlab的使用问题，自己上网搜索解决了全部问题，如多项式相乘函数，示波器背景颜色变为白色，将曲线导出到文档中等。仿真程序由两部分构成，分别是仿真调用函数m和逐个环节刷新的仿仿真程序w_DisNodesSimi。编程序的关键是理解离散相似法的原理，首先将各个矩阵求出来，然后将各个变量矩阵放进仿真调用函数即可。控制系统数字仿真与CAD实验三采样控制系统的数字仿真一、实验目的1、 掌握采样控制系统数字仿真的一般方法；2、 掌握采样周期与计算步长的关系与确定方法。二、 实验原理1、 采样控制系统：通常把系统中的离散信号是脉冲序列形成的离散系统,称为采样控制系统或脉冲控制系统；而把数字序列形成的离散系统,称为数字控制系统或计算机控制系统.2、 采样周期：在周期性测量过程变量（如温度、流量……）信号的系统中，相邻两次实测之间的时间间隔。离散控制系统（包括计算机数字控制系统）都采用周期性测量方式，采样间隔之内的变量值是不测量的。如采样周期过长，将引起有用信号的严重丢失，使系统品质变差。反之，如采样周期过短，则两次实测值的变化量太小，亦不相宜。采样周期的选择甚为重要，一般取为回复时间（即大体上达到稳态所需时间）的十分之一左右。三、 实验内容1、已知采样系统结构如图所示，分别用程序设计方法和Simulink法求系统的输出响应。（取仿真时间：Tf=10；采样周期:T=0.2；计算步长：h=0.01）（注：第一个框图中，分子分母中z的指数均为T；第二个框图中e的指数为-Ts）2、设某数字控制系统如图所示，采样周期为T=0.1s，初始状态可〔°）二花©=0（注：第二个框图中e的指数为-Ts）数字控制器为PI调节器，其中^=15j^=0-2D(z)= 七兰+瓦(1-z-1)数字控制器为PID调节器 ，其中^=15,^=0.2,^=40要求：1)用程序设计法、MATLAB控制工具箱时域响应分析函数、Simulink法求系统在单位阶跃输入信号=临)作用下的输出响应；2)比较两种控制器作用下系统的响应，由此得出微分调节的作用。四、实验步骤1、 打开电脑选择桌面上matlab。若桌面上没有则单机开始，选择程序，单击matlab程序，打开matlab程序。2、 新建一个M文件，并命名为“微分方程”。3、 按照前面的分析结果，分析系统整体的性能，进行编程，将编程的结果输入用户界面。4、 点击运行。5、 对程序代码中出现的错误进行改正，并重新运行。6、 观察结果，记录数据，进行参数变化，记录变化的结果。7、 将实验结果记录至之前的表格中。8、截图，按照要求完成实验报告。五、实验代码及结果clc;clear;G=[2.72-1];F=[0.717];P=[0110;1110];WIJ=[101;211];n=2;Y0=1;X0=[00];h=0.01;T0=0;Tf=10;Ts=0.2;nout=2;A=P(:,1);B=P(:,2);C=P(:,3);D=P(:,4);m=length(WIJ(:,1));W0=zeros(n,1);W=zeros(n,n);fork=1:mif(WIJ(k,2)==0);W0(WIJ(k,1))=WIJ(k,3);elseW(WIJ(k,1),WIJ(k,2))=WIJ(k,3);endendfori=1:nif(A(i)==0)FI(i)=1;FIM(i)=h*(C(i)/B(i));FIJ(i)=h*h*C(i)/B(i)/2;FIC(i)=1;FID(i)=0;if(D(i)~=0)FID(i)=D(i)/B(i);elseendelseFI(i)=exp((-h)*A(i)/B(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);elseendendendY=zeros(n,1);X=Y;y=0;Uk=zeros(n,1);Ub=Uk;U=0;uk=0;ek=0;E=zeros(n,1);x2=0;t=T0:h:Tf;t0=T0:h:Ts;N=length(t0);t1=T0:Ts:Tf;M=length(t1);fork=1:M-1;U=uk;ek=Y0-x2;E=[ek;E(1)];uk=-F*U+G*E;forl=1:N-1Ub=Uk;Uk=W*Y+W0*uk;Udot=(Uk-Ub)/h;Uf=2*Uk-Ub;X=FI'.*X+FIM'.*Uk+FIJ'.*Udot;Yb=Y;Y=FIC'.*X+FID'.*Uf;y(((N-1)*(k-1)+l),1)=Y(nout);endx2=Y(nout);endplot(1:1000,y);2.(1)G=[15.2-15];F=[-1];P=[1110;1410];WIJ=[101;211];n=2;Y0=1;X0=[00];h=0.01;L1=25;T0=0;Tf=10;Ts=0.1;nout=2;A=P(:,1);B=P(:,2);C=P(:,3);D=P(:,4);m=length(WIJ(:,1));W0=zeros(n,1);W=zeros(n,n);fork=1:mif(WIJ(k,2)==0);W0(WIJ(k,1))=WIJ(k,3);elseW(WIJ(k,1),WIJ(k,2))=WIJ(k,3);endendfori=1:nif(A(i)==0)FI(i)=1;FIM(i)=h*(C(i)/B(i));FIJ(i)=h*h*C(i)/B(i)/2;FIC(i)=1;FID(i)=0;if(D(i)~=0)FID(i)=D(i)/B(i);elseendelseFI(i)=exp((-h)*A(i)/B(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);elseendendendY=zeros(n,1);X=Y;y=0;Uk=zeros(n,1);Ub=Uk;U=0;uk=0;ek=0;E=zeros(2,1);x2=0;t=T0:h:Tf;t0=T0:h:Ts;N=length(t0);t1=T0:Ts:Tf;M=length(t1);fork=1:M-1;U=uk;ek=Y0-x2;E=[ek;E(1)];uk=-F*U+G*E;forl=1:N-1Ub=Uk;Uk=W*Y+W0*uk;Udot=(Uk-Ub)/h;Uf=2*Uk-Ub;X=FI'.*X+FIM'.*Uk+FIJ'.*Udot;Yb=Y;Y=FIC'.*X+FID'.*Uf;y(((N-1)*(k-1)+l),1)=Y(nout);endx2=Y(nout);endt=0:0.01:9.99;plot(t,y,'--');gridon;xlabel('t/s');ylabel('Amplitude');holdon(2)G=[55.2-9540];F=[-1];P=[1110;1410];WIJ=[101;211];n=2;Y0=1;X0=[00];h=0.01;L1=25;T0=0;Tf=10;Ts=0.1;nout=2;A=P(:,1);B=P(:,2);C=P(:,3);D=P(:,4);m=length(WIJ(:,1));W0=zeros(n,1);W=zeros(n,n);fork=1:mif(WIJ(k,2)==0);W0(WIJ(k,1))=WIJ(k,3);elseW(WIJ(k,1),WIJ(k,2))=WIJ(k,3);endendfori=1:nif(A(i)==0)FI(i)=1;FIM(i)=h*(C(i)/B(i));FIJ(i)=h*h*C(i)/B(i)/2;FIC(i)=1;FID(i)=0;if(D(i)~=0)FID(i)=D(i)/B(i);elseendelseFI(i)=exp((-h)*A(i)/B(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);elseendendendY=zeros(n,1);X=Y;y=0;Uk=zeros(n,1);Ub=Uk;U=0;uk=0;ek=0;E=zeros(3,1);x2=0;t=T0:h:Tf;t0=T0:h:Ts;N=length(t0);t1=T0:Ts:Tf;M=length(t1);fork=1:M-1;U=uk;ek=Y0-x2;E=[ek;E(1:2)];uk=-F*U+G*E;forl=1:N-1Ub=Uk;Uk=W*Y+W0*uk;Udot=(Uk-Ub)/h;Uf=2*Uk-Ub;X=FI'.*X+FIM'.*Uk+FIJ'.*Udot;Yb=Y;Y=FIC'.*X+FID'.*Uf;y(((N-1)*(k-1)+l),1)=Y(nout);endx2=Y(nout);endt=0:0.01:9.99;plot(t,y);gridon;xlabel('t/s');ylabel('Amplitude');六、实验小节通过本次实验我掌握了采样控制系统数字仿真的一般方法，掌握采样周期与计算步长的关系与确定方法。实验中遇到很多困难，但是在老师耐心的知道下我们还是完成了实验。本次实验要感谢老师耐心的指导和同学间的帮助，没有这些是不能完成的。通过此次实验我明白了控制系统数字仿真与CAD在实际中的应用是相当广泛的，通过在Matlab上的仿真实验，我们明白了各个参数变化对系统的影响，进一步体会到了采样周期与计算步长的关系。此次实验收益很多。

控制系统数字仿真与CAD】实验数字仿真技术的综合应用一、 实验目的1、 使学生在系统地学习了计算机仿真的基本概念和方法后。2、 对计算机仿真的全过程有一个清楚的概念和认识。3、 能够将所学的知识应用到实际系统中去。二、 实验原理一、一阶倒立摆动力学建模1.1直线一阶倒立摆数学模型的推导对系统建立数学模型是系统分析、设计的前提，而一个准确又简练的数学模型将大大简化后期的工作。为了简化系统分析，在实际的模型建立过程中，要忽略空气流动阻力，以及各种次要的摩擦阻力。这样，可将倒立摆系统抽象成小车本系统内部各相关参数定义如下：和匀质刚性杆组成的系统，如下图所示。本系统内部各相关参数定义如下：MmMmbl杆质心的长度IF摆杆质量小车摩擦系数摆杆转动轴心到摆杆惯量加在小车上的力小车位置0 摆杆与垂直向上方向的夹角图1-1小车及摆杆受力分析9 摆杆与垂直向下方向的夹角（考虑到摆杆初始位置为竖直向下）下图是系统中小车和摆杆的受力分析图。其中，N和P为小车与摆杆相互作用力的水平和垂直方向的分量。注意：在实际倒立摆系统中检测和执行装置的正负方向已经完全确定，因而矢量方向定义如图，图示方向为矢量正方向。应用Newton方法来建立系统的动力学方程过程如下：分析小车水平方向所受的合力，可以得到如下方程：TOC\o"1-5"\h\zMX=F-bx-N (1-1)由摆杆水平方向的受力进行分析可以得到下面等式：d2\o"CurrentDocument"N=m (x+lsin) (1-2)dt2\o"CurrentDocument"即 N=mx+ml^co@-ml°2si® (1-3)把这个等式代入上式中，就得到系统的第一个运动方程：\o"CurrentDocument"(M+m)x+bx+ml®cos0一ml®2sin0=F (1-4)为了推出系统的第二个运动方程，我们对摆杆垂直方向上的合力进行分析，可以得到下面方程：d2\o"CurrentDocument"P-mg=-m (lcos0) (1-5)dt2\o"CurrentDocument"即： P一mg=ml®si®+ml®2co0 (1-6)力矩平衡方程如下：\o"CurrentDocument"-Plsin®-Nlcos®=I® (1-7)注意：此方程中力矩的方向，由8=兀+申,cos©=-cos®,sin=-sin®，故等式前面有负号。合并这两个方程，约去P和N,得到第二个运动方程：\o"CurrentDocument"(I+ml2)®+mglsin®=-mlxcos® (1-8)1．微分方程模型设®=兀+申，当摆杆与垂直向上方向的夹角Q与1(单位是弧度)相比很小,即©vvl时，则可以进行近似处理：cose=-l,sin0=-^,( )2=0。为了与控制dt理论的表达习惯相统一，即u一般表示控制量，用u来代表被控对象的输入力F，线性化后得到该系统数学模型的微分方程表达式：I(I+ml2)(f)-mgl©=mix

[(M+m)x+bx-ml©=u2.传递函数模型对方程组(2-9)进行拉普拉斯变换，得到(l-9)(I+ml2)①(s)s2-mgl①(s)=mlX(s)s2(M+m)X(s)s2+bX(s)s-ml①(s)s2=U(s)注意：推导传递函数时假设初始条件为0。(l-l0)由于输出角度为©，求解方程组(2-10)的第一个方程，可以得到X(s)「(八ml2)-邑[①(s)mls2(1-11)①(s) mls2X(s) (I+ml2)s2一mgl(1-12)如果令v=x，则有:①(s) mlV(s) (I+ml2)s2一mgl(1-13)(I+ml2) g①(s)s2+b(I+ml2)*gml s2ml s2把上式代入方程组(1-10)的第二个方程，得到(M+m)①(s)s一ml①(s)s2=U(s)(1-14)整理后得到以输入力u为输入量，以摆杆摆角©为输出量的传递函数：mls2.(1-15)①(s)_ q.(1-15)U(s) b(I+ml2) (M+m)mgl bmglS4+ S3 S2一 Sqq其中q_[M+m)(I+ml2)-(mlJ实际参数如下：M 小车质量 0.5Kgm 摆杆质量 0.2Kgb 小车摩擦系数 0.1N/m/secl 摆杆转动轴心到杆质心的长度 0.3m0.006Kg*m*mI摆杆惯量0.006Kg*m*m把上述参数代入，可以得到系统的实际模型。摆杆角度和小车位移的传递函数为：1-16)1-17)1-18)①(s)_ 0.06s21-16)1-17)1-18)X(s)_0.024s2-0.588摆杆角度和小车加速度之间的传递函数为：①(s)_ 0.06V(sy_0.024s2-0.588摆杆角度和小车所受外界作用力的传递函数：①(s)_ 4.5455sU(s)_s3+0.