《计算机仿真技术》报告.doc

计算机仿真技术实验报告实验一 数字仿真方法验证一、实验目的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;ode45求解t-u图：（二）试用四阶RK法编程求解下列微分方程初值问题。仿真时间2s，取步长h=0.1。四阶RK法程序：12cleart=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; endytplot(t,y,r);结果：t=2,y=2.6109：（3） 试求示例3分别在周期为5s的方波信号和脉冲信号下的响应，仿真时间20s，采样周期Ts=0.1。已知系统的状态空间模型为程序：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)% square response of the systemu,t = gensig(square,5,20,0.1);subplot(212)lsim(sys,u,t) （四）某系统框图如图所示，试用SIMULINK进行仿真，并比较在无饱和非线性环节下系统仿真结果。模型：实验二 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% Create closed-loop feedback system modelsysclose=feedback(sys,1);% Quick access to zero-pole-gain data.z,p,k=zpkdata(sysclose,v);P结果：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);end% pidmargin.mfunction 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结果：ControllerKpTiTdP1.7849/PI1.42791.0882/PID2.14190.68010.16323计算比较原系统与P控制系统、PI控制系统、PID控制系统的瞬态性能指标。程序：sy23.m% Design PID controler% Create System Modelclfclearsys=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);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,SettlingTime 结果：SystemFinalvaluePeakTime,OvershootPercen