系统建模与仿真作业程序

四川理工学院系统建模与仿真课程报告学生：张三学号：11000000000专业：自动化班级：自动化000班任课教师：李四四川理工学院自动化与电子信息学院二O一三年十一月1目录第一章………………………………………………………21.1题目……………………………………………………21.2解答……………………………………………………2第二章………………………………………………………42.1题目………………………………………………………42.2程序………………………………………………………52.3小结……………………………………………………10第三章………………………………………………………103.1题目………………………………………………………103.2解答………………………………………………………10第四章………………………………………………………154.1题目……………………………………………………154.2解答……………………………………………………15第五章………………………………………………………165.1题目……………………………………………………165.2解答……………………………………………………162第一章1.1题目、x＝1+/6+/7+/10，)1)(1)(1)(1(*2)(sxsxsxsxxsG，采用面积法和最小二乘法把系统G等效成sesTk1特性，求系统的kT等参数。写出等效过程及相关程序。最后把等效后的两个模型及原系统对阶跃输入的响应曲线绘制在一个图上进行比较，并分析优劣问题。解题ruxuenian=input('请输入入学年：');banji=input('请输入班级：');xuehao=input('请输入学号：');rxn=num2str(ruxuenian);xh=num2str(xuehao);a=ruxuenian;b=banji;3c=str2num(xh(length(xh)));d=str2num(xh(length(xh)-1));x=1+a/6+d/7+c/10;s=[1:0.1:10];gs=2.*x./((x.*s+1).^4);figure(1)plot(s,gs);a0=[000];a=lsqcurvefit(@(a,s)a(1)./(a(2).*s+1).*exp(-a(3).*s),a0,s,gs)holdonplot(s,a(1)./(a(2).*s+1).*exp(-a(3).*s),'m');sys1=zpk([],[-1/x-1/x-1/x-1/x],2*x)subplot(1,2,1)step(sys1);[numt,dent]=pade(a(3),5);syst=tf(numt,dent);num1=a(1);den1=[a(2)1];sysT=tf(num1,den1);sys2=series(sysT,syst);subplot(1,2,2)step(sys2);请输入入学年：2011请输入班级：2请输入学号：11071020101Optimizationterminated:first-orderoptimalitylessthanOPTIONS.TolFun,andnonegative/zerocurvaturedetectedintrustregionmodel.a=000Zero/pole/gain:672.5333--------------(s+0.002974)^44第二章题目编写微分方程dy/dx＝xy,当x＝0时y＝1+7/10+11/100=1.81,x属于0～3之间，编写积分程序，包括欧拉数值积分程序，预报校正数字积分程序、4阶龙格库塔积分程序，它们的积分步长分别取0.01，0.1,0.5，绘制积分结果曲线,比较在同一步长下不同算法的误差和同一算法在不同步长下的误差，得出结论说明（绿色线为欧拉法曲线，红色为预报校正法曲线，蓝色为4阶龙格库塔法曲线）解：(1)4阶龙格库塔5Funval.mfunctionfv=Funval(f,varvec,varval)var=findsym(f);varc=findsym(varvec);s1=length(var);s2=length(varc);m=floor((s1-1)/3+1);varv=zeros(1,m);ifs1~=s2fori=0:((s1-1)/3)k=findstr(varc,var(3*i+1));index=(k-1)/3;varv(i+1)=varval(index+1);endfv=subs(f,var,varv);elsefv=subs(f,varvec,varval);endLGKT4.mfunctiony=LGKT4(f,h,a,b,y0,varvec)formatlong;N=(b-a)/h;y=zeros(N+1,1);y(1)=y0;x=a:h:b;var=findsym(f);fori=2:N+1K1=Funval(f,varvec,[x(i-1)y(i-1)]);K2=Funval(f,varvec,[x(i-1)+h/2y(i-1)+K1*h/2]);K3=Funval(f,varvec,[x(i-1)+h/2y(i-1)+K2*h/2]);K4=Funval(f,varvec,[x(i-1)+hy(i-1)+h*K3]);y(i)=y(i-1)+h*(K1+2*K2+2*K3+K4)/6;endformatshort;plot(x,y,'r:')legend('h=0.01','h=0.1','h=0.5')gridonholdonsymsxy;z=x*y;LGKT4(z,0.01,0,3,1.81,[xy])symsxy;z=x*y;LGKT4(z,0.1,0,3,1.81,[xy])symsxy;z=x*y;LGKT4(z,0.5,0,3,1.81，[xy])600.511.522.53020406080100120140160180h=0.0100.511.522.53020406080100120140160180h=0.500.511.522.53020406080100120140160180小结：在用4阶龙哥库塔积分结果图像可以得出：在积分过程中，步长越短，其积分曲线越平滑即其近似程度越高，精度也就越高。（2）欧拉法Euler.m7functiony=Euler(f,h,a,b,y0,varvec)%±í´ïʽ£ºf%»ý·Ö²½³¤£ºh%×ÔÏÂÏÞ£ºa%×ÔÉÏÏÞ£ºb%³õÖµ£ºy0%±äÁ¿×飺varvecformatlong;N=(b-a)/h;y=zeros(N+1,1);y(1)=y0;x=a:h:b;fori=2:N+1y(i)=y(i-1)+h*Funval(f,varvec,[x(i-1),y(i-1)]);endformatshort;subplot(2,2,1)plot(x,y)legend('ou')gridonsymsxy;f=x*y;Euler(f,0.01,0,3,[xy])symsxy;f=x*y;Euler(f,0.1,0,3,[xy])symsxy;f=x*y;Euler(f,0.5,0,3,[xy])00.511.522.53020406080100120140160180800.511.522.5302040608010012000.511.522.53051015202530小结：在用欧拉积分法中，步长越短，精度越高！！预报校正functiony=yubaojiaozheng(f,h,a,b,y0,varvec,type,s)formatlong;N=(b-a)/h;y=zeros(N+1,1);y(1)=y0;x=a:h:b;y(2)=y0+h*Funval(f,varvec,[x(1)y(1)]);var=findsym(f);iftype==1fori=3:N+1v1=Funval(f,varvec,[x(i-1)y(i-1)]);t=y(i-2)+2*h*v1;9v2=Funval(f,varvec,[x(i)t]);y(i)=y(i-1)+h*(v1+v2)/2;endendformatshort;subplot(2,2,3)plot(x,y)legend('yu')gridonsymsxy;z=x*y;yubaojiaozheng(z,0.01,0,3,1.81,[xy],1,3)symsxy;z=x*y;yubaojiaozheng(z,0.1,0,3,1.81,[xy],1,3)symsxy;z=x*y;yubaojiaozheng(z,0.5,0,3,1.81,[xy],1,3)00.511.522.5302040608010012014016018000.511.522.530204060801001201401601801000.511.522.53050100150总结结论：在同一种算法下其步长越大误差越大；在同一步长下欧拉法误差相对较小。在数值积分中要根据不同的情况取不同的积分算法，同时在取步长的时候也要综合考虑，虽然步长取得越小越精确但是这样加大计算量，需要更长时间，所以应视情况综合考虑。第三章3.2解题（1）num=[2.1*1.2];den=conv([3.211],[133.21]);rlocus(tf(num,den))11（2）num=[2.1*1.2];den=conv([3.211],[133.21]);rlocus(tf(num,den))den=conv([3.211],[133.21]);rlocus(tf(num,den))num=[2.1*1.2];sys=tf(num,den,0.1)Transferfunction:2.52------------------------------------3.21z^3+10.63z^2+13.3z+3.21Samplingtime:0.1（3）step(num,den)[numcdenc]=cloop(num,den)12numc=0002.5200denc=3.210010.630013.30415.7300g=tf(numc,denc)Transferfunction:2.52------------------------------------3.21s^3+10.63s^2+13.3s+5.73margin(numc,denc)（4）s=tf('s');G1=2.1/(3.21*s+1);13G2=1.2/(s*s+s*3+3.21);G=G1*G2Transferfunction:2.52------------------------------------3.21s^3+10.63s^2+13.3s+3.21sisotool（5）num1=[2.1];den1=[3.211];num2=[1.2];den2=[133.21];[num3den3]=series(num1,den1,num2,den2);[num4den4]=series(num3,den3,num,den);num=[2.1*1.2];den=conv([3.211],[133.21]);[num4den4]=series(num3,den3,num,den);[numcdenc]=cloop(num4,den4);F=tf(numc,denc)Transferfunction:6.35--------------------------------------------------------------------------10.3s^6+68.24s^5+198.4s^4+303.5s^3+245.2s^2+85.41s+16.65[zpk]=tf2zp(numc,denc)14z