控制系统数字仿真作业.doc

控 制 系 统 数 字 仿 真 作 业学号：S309047020 姓名：张宇2010年1月5号2.2 将化成状态空间表达式，要求A阵为对角阵。解：=令： ， 则：将以上两式拉氏反变换： 所以可得： y（4 2）2.3 已知系统的状态空间描述为：A B C=求该系统的传递函数。解：=+=2.4 求图2-20所示系统的状态空间描述。图2-20解：所以可以写出系统的状态空间描述为：A B C2.5分别用欧拉法、二阶龙格库塔及四阶龙格库塔法计算系统：在阶跃函数下的过渡过程1. 选择相同的步距h0.005，试比较计算结果。2. 选择不同的步距：欧拉法h0.0001，二阶龙格库塔法h0.005，四阶龙格库塔法h0.05，试比较计算结果。解：= 求得：=可以写出系统的状态空间方程为： 即：=，=，=，利用MATLAB编写了程序，见附录1、附录2和附录3，1. 选择相同的步距h0.005，得到的三者的图形如下：图 1 欧拉法 图 2 2阶龙格库塔法 图 3 4阶龙格库塔法时间（s）欧拉法二阶四阶时间（s）欧拉法二阶四阶0.0500.00416660.00411781.050.976460.977970.978410.100.385340.413450.41231.100.96930.968150.967830.151.07141.09461.09311.150.983360.980810.980120.201.51191.51921.5191.201.00311.00111.00060.251.48231.47881.48081.251.01351.01341.01340.301.13751.13231.13511.301.01041.01181.01220.350.799320.796160.797551.350.999551.00121.00160.400.692320.690040.688951.400.990480.991180.991360.450.821650.818770.816131.450.988690.988150.988020.501.03261.03031.02811.500.993240.992090.99180.551.16031.16081.16071.550.999190.998350.998140.601.1431.1471.14871.601.00211.00211.00210.651.0341.03881.04091.651.00091.00151.00160.700.932680.93480.935721.700.997460.998140.998310.750.904630.902560.901881.750.994720.9950.995070.800.947530.942990.941471.800.994270.994040.993980.851.01221.00861.00751.850.99570.995230.995120.901.04891.04871.04871.900.997460.997130.997050.951.04111.0441.04491.950.998240.998230.998231.001.00681.01021.01122.000.997750.997990.99805从数据比较看出，三种方法均在0.2s达到峰值，其中欧拉法的峰值最小为1.5119，4阶龙格库塔法的最小为1.519，但4阶龙格库塔法的稳定时间最短，欧拉法的最长。2. 选择不同的步距：欧拉法h0.0001，二阶龙格库塔法h0.005，四阶龙格库塔法h0.05， 图 4 欧拉法 图 5 2阶龙格库塔法 图 6 4阶龙格库塔法时间（s）欧拉法二阶四阶时间（s）欧拉法二阶四阶0.0500.0041670.345281.050.980790.968150.968740.100.349990.413451.02111.100.967730.980810.976250.151.03031.09461.49271.150.978321.00110.996150.201.49741.51921.50951.200.99881.01341.01190.251.50281.47881.18881.251.01281.01181.01430.301.17371.13230.835261.301.01281.00121.00550.350.822580.796160.688971.351.00280.991180.994390.400.687730.690040.7871.400.99210.988150.988790.450.796990.818770.996321.450.9880.992090.990760.501.00881.03031.14871.500.991250.998350.996840.551.1541.16081.16071.550.997591.00211.00180.601.15541.1471.06391.601.00191.00151.00270.651.05281.03880.952261.651.00180.998141.00010.700.943440.93480.902791.700.998660.9950.996530.750.901520.902560.930221.750.99530.994040.994630.800.935590.942990.994881.800.993980.995230.99510.851.00151.00861.0441.850.994950.997130.996910.901.04661.0441.04981.900.996880.998230.998470.951.04691.01021.02061.950.998170.997990.998771.001.01490.977970.985392.000.998110.997930.99793结果：欧拉法的步距减小之后其超调和峰值都有都随之减小了，在各个时刻上对应的值都比h=0.005时要小，4阶龙格库塔法在步长调大之后响应过程变快了，峰值时间和幅度小了。3.1试求系统：的及。解：=+可得 即A= B= = = =3.2试求系统：的及（=c2 if y4=c2 f2c=y4-c2; else f2c=y4+c2; end else f2c=0; end u1=y0-f2c; if (abs(y1)=c1 if y1=c1 f1c=c1; else f1c=-c1 end else f1c=y1; end u2=f1c; x1=faiT1*x1+faiTM1*u1+faiTM_pie1*(u1-u10); y1=(0.5-0.1)*x1+u1; x2=faiT2*x2+faiTM2*u2+faiTM_pie2*(u2-u20); y2=x2; x3=faiT3*x3+faiTM3*u3+faiTM_pie3*(u3-u30); y3=x3; x4=faiT4*x4+faiTM4*u4+faiTM_pie4*(u4-u40); y4=x4; u10=u1; u20=u2; u30=u3; u40=u4; y=f2c; Y=Y,y;endplot(Y)grid on;xlabel(ms);附5 单纯形法程序附5.1 计算目标函数程序function Q=Q(k1,T1)k1=0.3;T1=95;num=10*k1*T1 0;den=2*T1 2+T1 1 0;a b c d=tf2ss(num,den);T=0.005;h=T;N=5000;x=0;0;0;R=1*ones(1,N);yy=1*ones(1,N);yy(1,1)=0;for k=2:N u=1-yy(1,k-1); k0=a*x+b*u; k2=a*(x+k0*h/2)+b*u; k3=a*(x+k2*h/2)+b*u; k4=a*(x+k3*h)+b*u; x=x+(k0+2*k2+2*k3+k4)*h/6; y=c*x+d*u; yy(1,k)=y;endQ=0;for k=1:N f1=yy(1,k); f2=abs(1-f1); Q=Q+k*T*T*f2;end附5.2 单纯形法主程序clear;clc;N=2;h=0.05;mu=2;lamda=0.4;k_max=20;%最大搜索次数k=0;%当前搜索次数disp(初始点);k1=0.3;T1=95;alpha=zeros(N+1,N+1);alpha(1,1)=k1;alpha(2,1)=T1;e=eye(N,N);alpha_h=zeros(2,1);%最坏点alpha_l=zeros(2,1);%最好点alpha_g=zeros(2,1);%次坏点alpha_r=zeros(2,1);%反射点alpha_e=zeros(2,1);%扩张点alpha_s=zeros(2,1);%压缩点%2.计算单纯形的初始点for i=1:N for j=2:N+1 alpha(i,j)=alpha(i,1)+h*e(i,j-1); endend%3.计算每个初始点对应的函数值for i=1:N+1 alpha(3,i)=Q(alpha(1,i),alpha(2,i);end alpha%4.找出最坏点、最好点和次坏点for i=1:N+1%函数值从大到小排列 for j=i+1:N+1 if alpha(3,i)e&k=k_max k=k+1 aH=aH,alpha_h; aL=aL,alpha_l; alpha_r=alpha_g+alpha_l-alpha_h;%9.计算反射点 c_r=Q(alpha_r(1,1),alpha_r(2,1); if c_rc_g u=(1-mu)*c_h+mu*c_r; if uc_l alpha_e=(1-mu)*alpha_h+mu*alpha_r; c_e=Q(alpha_e(1,1),alpha_e(2,1); if c_ec_r alpha_s=alpha_e; c_s=c_e; disp(扩张); else alpha_s=alpha_r; c_s=c_r; disp(反射点); end else alpha_s=alpha_r; c_s=c_r; disp(反射点); end alpha_h=alpha_s c_h=c_s if c_hc_g%重新找出最坏点、最好点和次坏点 l_Q=c_h; l_s=alpha_h; c_h=c_g; alpha_h=alpha_g; c_g=l_Q; alpha_g=l_s; end if c_gc_l l_Q=c_g; l_s=alpha_g; c_g=c_l; alpha_g=alpha_l; c_l=l_Q; alpha_l=l_s; end else if c_rc_h alpha_s=(1-lamda)*alpha_h+lamda*alpha_r; c_s=Q(alpha_s(1,1),alpha_s(2,1); else alpha_s=lamda*alpha_h+(1-lamda)*alpha_r; c_s=Q(alpha_s(1,1),alpha_s(2,1); end if c_sc_g disp(压缩); alpha_h=alpha_s c_h=c_s if c_hc_g%重新找出最坏点、最好点和次坏点 l_Q=c_h; l_s=alpha_h; c_h=c_g; alpha_h=alpha_g; c_g=l_Q; alpha_g=l_s; end if c_gc_l l_Q=c_g; l_s=alpha_g; c_g=c_l; alpha_g=alpha_l; c_l=l_Q; alpha_l=l_s; end else alpha_h=(alpha_l+alpha_h)/2; alpha_g=(alpha_l+alpha_g)/2; c_h=Q(alpha_h(1,1),alpha_h(2,1); c_g=Q(alpha_g(1,1),alpha_g(2,1); disp(收缩); alpha_h c_h alpha_g c_g if c_hc_g%重新找出最坏点、最好点和次坏点 l_Q=c_h; l_s=alpha_h; c_h=c_g; alpha_h=alpha_g; c_g=l_Q; alpha_g=l_s; end if c_gc_l l_Q=c_g; l_s=alpha_g; c_g=c_l; alpha_g=alpha_l; c_l=l_Q; alpha_l=l_s; end end end s=abs(c_h-c_l)/abs(c_l); end disp(最优点); alpha_l c_l Q1(k1,T1); Q2(alpha_l(1,1),alpha_l(2,1); ah_K=;ah_T=;al_K=;al_T=; for j=1:2:24 ah_K=ah_K,aH(j); ah_T=ah_T,aH(j+1); al_K=al_K,aL(j); al_T=al_T,aL(j+1); end figure(2) grid on subplot(2,2,1) plot(ah_K,-rs,LineWidth,2,. MarkerEdgeColor,k,. MarkerFaceColor,g,. MarkerSize,3) grid on; subplot(2,2,2) plot(ah_T,-rs,LineWidth,2,. MarkerEdgeColor,k,. MarkerFaceColor,g,. MarkerSize,3) grid on; subplot(2,2,3) plot(al_K,-rs,LineWidth,2,. MarkerEdgeColor,k,. MarkerFaceColor,g,. MarkerSize,3) grid on; subplot(2,2,4) plot(al_T,-rs,LineWidth,2,. MarkerEdgeColor,k,. MarkerFaceColor,g,. MarkerSize,3) grid on;附5.3 寻优前动态仿真程序function Q=Q1(k1,T1)k1=0.3;T1=95;num=10*k1*T1 0;den=2*T1 2+T1 1 0;a b c d=tf2ss(num,den);T=0.005;h=0.005;N=10000;x=0;0;0;R=1*ones(1,N);yy