《计算机仿真技术》报告.pdf
计算机仿真技术实验报告2 实验一数字仿真方法验证一、实验目的1掌握基于数值积分法的系统仿真、了解各仿真参数的影响；2掌握基于离散相似法的系统仿真、了解各仿真参数的影响；3掌握 SIMULINK 动态仿真；4熟悉 MATLAB 语言及应用环境。二、实验环境网络计算机系统，MATLAB 语言环境三、实验内容、要求（一）试将示例 1 的问题改为调用 ode45函数求解，并比较结果。示例 1：设方程如下，取步长h=0.1。上机用如下程序可求出数值解。调用 ode45函数求解：1）建立一阶微分方程组du=u-2*t/u2）建立描述微分方程组的函数m文件function du=sy11vdp(t,u)du=u-2*t/u3）调用解题器指令 ode45求解yt,u=ode45(sy11vdp,0 1,1)plot(t,u,r-);xlabel(t);ylabel(u);结果对比：euler法：t=1,u=1.7848;RK法：t=1,u=1.7321;ode45求解:t=1,u=1.7321;1,01)0(2tuutudtdu3 ode45求解t-u图：00.10.20.30.40.50.60.70.80.9111.11.21.31.41.51.61.71.8tu（二）试用四阶RK 法编程求解下列微分方程初值问题。仿真时间2s，取步长h=0.1。1)0(2ytydtdy四阶 RK 法程序：cleart=2;h=0.1;n=t/h;t0=0;y0=1;y(1)=y0;t(1)=t0;for i=0:n-1 k1=y0-t02;k2=(y0+h*k1/2)-(t0+h/2)2;k3=(y0+h*k2/2)-(t0+h/2)2 k4=(y0+h*k3)-(t0+h)2;y1=y0+h*(k1+2*k2+2*k3+k4)/6;t1=t0+h;y0=y1;t0=t1;y(i+2)=y1;t(i+2)=t1;4 endytplot(t,y,r);结果：t=2,y=2.6109 00.511.522.511.21.41.61.822.22.42.62.83：（三）试求示例 3 分别在周期为 5s的方波信号和脉冲信号下的响应，仿真时间20s，采样周期 Ts=0.1。已知系统的状态空间模型为uXXXX0107814.07814.05572.02121214493.69691.1XXY程序：clear%Create system modelA=-0.5572-0.7814;0.7814 0;B=1;0;C=1.9691 6.4493;D=0;sys=ss(A,B,C,D);%impulse response of the systemsubplot(211)impulse(sys,0:0.1:20)5%square response of the systemu,t=gensig(square,5,20,0.1);subplot(212)lsim(sys,u,t)02468101214161820-20246Impulse ResponseTime(sec)Amplitude024681012141618200510Linear Simulation ResultsTime(sec)Amplitude（四）某系统框图如图所示，试用SIMULINK进行仿真，并比较在无饱和非线性环节下系统仿真结果。模型：6 2Out11Out220s +12s +20s32Transfer Fcn3s+0.5s+0.1Transfer Fcn220s +12s +20s32Transfer Fcn1s+0.5s+0.1Transfer FcnStepScope3Scope2Scope1Scope1.5Gain11.5Gain0.501234567891000.10.20.30.40.50.60.70.80.91tyy1(saturation)y27 8 实验二PID 控制器设计一、实验目的1了解 PID 控制原理，掌握相应PID 控制器设计仿真程序的应用；2掌握计算机辅助系统瞬态性能指标的计算；3掌握计算机辅助系统频率性能分析；二、实验环境网络计算机系统（采矿楼四楼测试实验室），MATLAB 语言环境三、实验内容、要求已知如图所示单位反馈系统要求：1绘制系统的开环Nyquist 图和 Bode 图，并判断该闭环系统是否稳定。程序：Sy21.m%Create system modelsys=tf(500 5000,1 33 337 1775 4950 5000);%Draw the Nyquist plotfigure(1)nyquist(sys)grid onRe,Im,w1=nyquist(sys);%Draw the Bode plotfigure(2)bode(sys)mag,phase,w2=bode(sys);grid on 5000495017753373350005002345ssssss9%Create closed-loop feedback system modelsysclose=feedback(sys,1);%Quick access to zero-pole-gain data.z,p,k=zpkdata(sysclose,v);P 结果：-1-0.8-0.6-0.4-0.200.20.40.60.81-1-0.8-0.6-0.4-0.200.20.40.60.810 dB-20 dB-10 dB-6 dB-4 dB-2 dB20 dB10 dB6 dB4 dB2 dBNyquist DiagramReal AxisImaginaryAxis-200-150-100-500Magnitude(dB)10-210-1100101102103-360-270-180-900Phase(deg)Bode DiagramFrequency (rad/sec)p=-19.9388-5.1469+2.5108i-5.1469-2.5108i-1.3837+3.6577i-1.3837-3.6577i 无闭环右极点，系统稳定。2应用 ZieglerNichols 方法设计 P 控制器、PI 控制器和 PID 控制器。程序：sy22.m%Design PID controler%Create System Modelclf,clcsys=tf(500 5000,1 33 337 1775 4950 5000);sysgroup=feedback(sys,1);%Design P-PI-PIDControlerfor i=1:3 type=i;sysc,Kp,Ti,Td=pidmargin(sys,type);sysopen=sysc*sys;sysclose=feedback(sysopen,1);sysgroup=append(sysgroup,sysclose);10 end%pidmargin.m function sysc,Kp,Ti,Td=pidmargin(sys,type)margin(sys)Gm,Pm,Wg,Wc=margin(sys);Kcr=Gm;Wcr=Wg;Tcr=2*pi/Wcr;switch typecase 1 disp(P Controler)Kp=0.5*Kcr Ti=No Design Td=No Design sysc=Kp;case 2 disp(PI Controler)Kp=0.4*Kcr Ti=0.8*Tcr Td=No Design sysc=Kp*(1+tf(1,Ti,0);case 3 disp(PID Controler)Kp=0.6*Kcr Ti=0.5*Tcr Td=0.12*Tcr sysc=Kp*(1+tf(1,Ti,0)+tf(Td,0,1);end结果：Controller Kp Ti Td P 1.7849/PI 1.4279 1.0882/PID 2.1419 0.6801 0.1632 3计算比较原系统与P 控制系统、PI 控制系统、PID 控制系统的瞬态性能指标。程序：sy23.m%Design PID controler%Create System Modelclfclearsys=tf(500 5000,1 33 337 1775 4950 5000);11 sysgroup=feedback(sys,1);%Design P-PI-PIDControlerfor i=1:3 type=i;sysc,Kp,Ti,Td=pidmargin(sys,type);sysopen=sysc*sys;sysclose=feedback(sysopen,1);sysgroup=append(sysgroup,sysclose);endclffor i=1:4 subplot(2,2,i)step(sysgroup(i,i)pidsys(i)=sysgroup(i,i);%This program performs an analysis of step response of Design P-PI-PIDControlernum,den=tfdata(pidsys(i),v);%Compute steady valueFinalvalue(i)=polyval(num,0)/polyval(den,0);%Compute overshooty,t=step(pidsys(i);%Ymax,k=max(y);PeakTime(i)=t(k);%OvershootPercent(i)=100*(Ymax-Finalvalue)/Finalvalue;%Compute rise timen=1;while y(n)0.1*Finalvalue,n=n+1;endm=1;while y(m)0.98*Finalvalue&y(r)1.02*Finalvalue)r=r-1;endSettlingTime(i)=t(r);%endFinalvalue,PeakTime,OvershootPercent,RiseTime,SettlingTime12 结果：0123400.20.40.60.8sysTime(sec)AmplitudeP-sysTime(sec)Amplitude0246800.511.5 PI-sysTime(sec)AmplitudePID-sysTime(sec)Amplitude0246800.5101234500.511.5System Finalvalue PeakTime,OvershootPercent RiseTime SettlingTime sys 0.5000 1.1572 22.9897 0.4788 2.9528 P_sys 0.6409 1.0266 59.9731 0.2933 7.9680 PI_sys 1.0000 1.1604 45.2370 0.3094 7.9678 PID_sys 1.0000 0.8954 45.7340 0.2132 4.9886
