系统仿真上机作业.docx

1、系统仿真上机作业系统仿真上机作业学号：04103166姓名：丛旭亚系统仿真上机作业一、 计算机辅助系统分析：系统如下图所示其中 r 是单位阶跃， GN 是非线性器件，G0s=K(1+s)s(10s+1)(0.625s+1)(0.025s+1)1.当 GN=1 、K = 40时，用 MATLAB 画出开环 Bode 图，求出 C、B。 由其估计出tr 、ts 、 %。解：使用程序如下：num=40 40;den=conv(conv(10 1 0,0.625 1),0.025 1);bode(num,den);grid onG=tf(num,den);Gm Pm Wcg Wcp=margin(G)

2、运行程序得到如下结果：Gm =4.3241Pm =10.0298Wcg =5.1672Wcp =2.3985其中C= Wcp =2.3985rad/s，B= Pm =10.0298绘制出Bode图如下：通过系统的开环Bode图，可将系统的开环传递函数简化成G0s=40s(10s+1)，其相应的闭环传递函数为GB(s)=4010s2+s+40。可求得=0.025，n=2。所以可估计出 %=e-/1-2*100%=92.44%，ts=3.5/(n)=70s,tr=-n1-2=0.80s(=arctan1-2)。 2. 当 GN=1 、K = 40时，用 MATLAB 画出根轨迹图，并求出K = 4

3、0时的闭环极点； 由其估计出tr 、ts 、 %。解：使用程序如下：num=40 40;den=conv(conv(10 1 0,0.625 1),0.025 1);rlocus(num,den)G=tf(num,den);G1=feedback(G,1)p,z=pzmap(G1)运行程序得到如下结果：p = -40.1616 -0.2274 + 2.4146i -0.2274 - 2.4146i -1.0837 z =-1其中得到的p的值就是K=40时的闭环极点。绘出根轨迹图如下：由求得的闭环传递函数GB(s)=40（s+1）s+40.1616 s+0.2274- 2.4146is+0.22

4、74+ 2.4146i(s+1.0837 ) ,将其可近似为一个二阶系统GBs=40s+0.2274- 2.4146is+0.2274+ 2.4146i5.882 s2 + 0.4548 s + 5.882。可求得=0.0938，n=2.4253。所以可估计出 %=e-/1-2*100%=74.38%，ts=3.5/(n)=15.39s,tr=-n1-2=0.69s(=arctan1-2)。3. 当GN=1 ，K = 40时，仿真之，并由仿真结果求出tr 、ts 、 %。用自适应变步长法仿真结果如下：由仿真结果可知 tr=0.44s，ts=17.2s(=2%)， %=81%用定步长 RK-2

5、法，并观察步长与仿真结果的收敛于发散关系，以及仿真精度。当步长h=0.0498时，仿真结果收敛，仿真图如下：由仿真结果可知 tr=0.42s，ts=18.32s(=2%)， %=81%当步长h=0.0499时,仿真结果发散，仿真图如下：4. 令图中的K = 40。 GN分别为：分别仿真之，并由仿真结果求出 tr、ts 、 %。当GN为饱和特性时，仿真结果如下： tr=1.37s、ts= 17.6s、 %=56%。当GN为死区特性时，仿真结果如下： tr=0.44s，ts= 32.5s、 %=100%。GN=1，在 G0(s) 之后，反馈点之前加上。仿真之，并由仿真结果求出 tr、ts 、 %。

6、 tr=0.4s，ts=18.25s、 %=100%。对 3 和4 中、的 tr、ts 、 %。比较，并解释差异的原因。1.通过比较第三问中自适应变步长法和定步长RK-2法的仿真结果发现，后者的ts比前者的ts要长，而前后仿真的 tr、 %相差较小。其原因是自适应变步长法中步长hk的大小与y的变化快慢有关，当系统趋于稳定时，y变化很慢，hk变大，而定步长时hk一直为一个较小值，所以定步长RK-2法的ts要比自适应变步长法的ts长。而在t较小的一段时间内,自适应变步长法和定步长RK-2法的步长hk都比较小，所以它们的 tr、 %相差不大。2.4中的非线性环节加在G0(s)之前，相当于死区特性作用

7、在输入e(s)=r(s)-y(s)上，所以对原系统的动态性能影响较大， tr、ts 、 %与原系统相差较多；而中的非线性环节加在G0(s)之后，则死区特性并没有直接作用于输入e(s)上，而是作用于e(s) G0(s)上，所以对原系统的动态特性影响较小， tr、ts 、 %与原系统相差较小。二、病态系统（stiff）仿真（simulink） R y r :单位阶跃；G(s)=1(10-41s+1)(1042s+1)1.用自适应变步长法（RK45）仿真之解：当1=103，2=10-3时。仿真结果如下：2.用定步长四阶龙格库塔法仿真，并试着搜索收敛的步长h 的范围；若找不到h，将1 增大，2减小，用

8、定步长四阶龙格库塔法仿真，寻找h。解：当1=103，2=10-3时。收敛的步长h的范围是：h50，该系统是一个病态系统，所以用病态仿真算法仿真最精确，而RK45和RK4则都会产生一定的误差。三计算机辅助控制器设计： A R（s） + Y(s) - B 其中 Gs=10s(s+1) 要求：开环B45，c4.2，且tr0.4s，ts1.5s， %25% 1. 开关处于 A 时，系统性能满足上述要求否？解：编写如下程序求解系统的各项指标如下：num=10;den=1 1 0;G=tf(num,den);Gm Pm Wcg Wcp=margin(G)num=10;den=1 1 0;sys=tf(nu

9、m,den);sys=feedback(sys,1);y,t=step(sys);ytr=find(y=1);risetime=t(ytr(1)ymax,tp=max(y);maxovershoot=ymax-1s=length(t);while y(s)0.98&y(s)1.02s=s-1;endsettlingtime=t(s+1)程序运行结果为：Gm =InfPm =17.9642Wcg =InfWcp = 3.0842risetime =0.6447maxovershoot = 0.6045settlingtime=7.2762由此可见B=Pm=17.964245，c= 3.08420

10、.4s，ts= 7.2762s1.5s， %=60.45%25%，均不符合要求。2.开关处于 B 时，计算机辅助设计GC(s)，使系统性能满足上述要求。由第一问可知，要使系统满足上述要求必须串联超前校正环节。设计频率法超前校正的MATLAB子函数如下：function Gc=plsj(G,kc,yPm)G=tf(G);mag,pha,w=bode(G*kc);Mag=20*log10(mag);Gm,Pm.Wcg,Wcp=margin(G*kc);phi=(yPm-getfield(Pm,Wcg)*pi/180;alpha=(1+sin(phi)/(1-sin(phi);Mn=-10*log1

11、0(alpha);Wcgn=spline(Mag,w,Mn);T=1/(Wcgn*sqrt(alpha);Tz=alpha*T;Gc=tf(Tz 1,T 1);End主程序如下：num=10;den=1 1 0;G=tf(num,den);kc=2;the greater, the better!yPm=18+40;Gc=plsj(G,kc,yPm)这一步得到的串联超前校正环节传递函数GCs= 0.3503 s + 10.05938 s + 1，校正环节的增益为2。Gc1=feedback(G*Gc*kc,1);Gm,Pm,Wcg,Wcp=margin(G*Gc*kc)y,t=step(Gc1

12、);ytr=find(y=1);risetime=t(ytr(1)ymax,tp=max(y);maxovershoot=ymax-1s=length(t);while y(s)0.98&y(s)45，c= 6.93374.2，tr=0.2554s0.4s，ts=0.9398s1.5s， %=20.21%25%，均符合要求。校正后仿真结果如下图所示：四、一级倒立摆仿真收敛性研究仿真在位移x处加阶跃干扰的倒立摆非线性模型（已建好，但仿真不成功）， 除了倒立摆系统和控制器不能改动，其他地方可以酌情修改，如更改算法等等。要求： 1.仿真收敛； 2.说明修改过程并解释物理意义； 3.说明如何验证仿真结果（不）稳定是系统本身（不）稳定还是仿真原因造成的？解：在位移X处加如下阶跃干扰，如下图所示： 其物理意义是：在时间t=1s时,使小车的位移发生由0到0.15的突变。在t=0时刻，在输入端加一个单位阶跃信号作为输入的外力F。采用自适应变步长ode45算法，得到的仿真结果如下：位移x:摆杆与竖直方向的夹角:有仿真结果可知仿真收敛。通过理论计算可以得到系统本身的稳定性，若该结果与仿真结果一致，则可以验证仿真结果的（不）稳定是系统（不）稳定造成的；若该结果图仿真结果不一致，则可以验证仿真结果的（不）稳定是仿真原因造成的。