《计算机仿真技术》实验报告实验一 数字仿真方法验证一、实验目的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求解y[t,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(2∈⎪⎩⎪⎨⎧=-=t u u t u dt duode45求解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(2y t y dt dy 四阶RK 法程序：clear t=2; h=0.1; n=t/h; t0=0; y0=1;y(1)=y0; t(1)=t0;for i=0:n-1 k1=y0-t0^2;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;end y tplot(t,y,'r'); 结果：t=2,y=2.61090.511.522.511.21.41.61.822.22.42.62.83：（三）试求示例3分别在周期为5s 的方波信号和脉冲信号下的响应，仿真时间20s ，采样周期Ts=0.1。
已知系统的状态空间模型为u X X X X ⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡--=⎥⎦⎤⎢⎣⎡0107814.07814.05572.02121[]⎥⎦⎤⎢⎣⎡=214493.69691.1X X Y程序：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 system subplot(211)impulse(sys,[0:0.1:20])% square response of the system[u,t] = gensig('square',5,20,0.1); subplot(212) lsim(sys,u,t)02468101214161820-20246Impulse ResponseTime (sec)A m p l i t u d e02468101214161820510Linear Simulation ResultsTime (sec)A m p l i t u d e（四）某系统框图如图所示，试用SIMULINK 进行仿真，并比较在无饱和非线性环节下系统仿真结果。
模型：2Out11Out220s +12s +20s 32Transfer Fcn3s+0.5s+0.1Transfer Fcn220s +12s +20s 32Transfer Fcn1s+0.5s+0.1Transfer FcnStepScope3Scope2Scope1Scope1.5Gain11.5Gain0.51234567891000.10.20.30.40.50.60.70.80.91tyy1(saturation)y2实验二 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 plot figure(1)nyquist(sys) grid on[Re,Im,w1]=nyquist(sys); % Draw the Bode plot figure(2) bode(sys)[mag,phase,w2]=bode(sys); grid on5000495017753373350005002345++++++s s s s s s% Create closed-loop feedback system model sysclose=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 dB 10 dB 6 dB 4 dB2 dBNyquist DiagramReal AxisI m a g i n a r y A x is-200-150-100-50M a g n i t u d e (d B )10-210-110101102103-360-270-180-900P h a s e (d e g )Bode DiagramFrequency (rad/sec)p = -19.9388 -5.1469 + 2.5108i -5.1469 - 2.5108i -1.3837 + 3.6577i -1.3837 - 3.6577i无闭环右极点，系统稳定。
2．应用Ziegler —Nichols 方法设计P 控制器、PI 控制器和PID 控制器。
程序：sy22.m% Design PID controler % Create System Model clf,clcsys=tf([500 5000],[1 33 337 1775 4950 5000]); sysgroup=feedback(sys,1); %Design P-PI-PIDControler for 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 1disp('P Controler')Kp=0.5*KcrTi='No Design'Td='No Design'sysc=Kp;case 2disp('PI Controler')Kp=0.4*KcrTi=0.8*TcrTd='No Design'sysc=Kp*(1+tf(1,[Ti,0]));case 3disp('PID Controler')Kp=0.6*KcrTi=0.5*TcrTd=0.12*Tcrsysc=Kp*(1+tf(1,[Ti,0])+tf([Td,0],1));end结果：Controller\ Kp Ti TdP 1.7849 / /PI 1.4279 1.0882 /PID 2.1419 0.6801 0.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:3type=i;[sysc,Kp,Ti,Td]=pidmargin(sys,type);sysopen=sysc*sys;sysclose=feedback(sysopen,1);sysgroup=append(sysgroup,sysclose);endclffor i=1:4subplot(2,2,i)step(sysgroup(i,i))pidsys(i)=sysgroup(i,i); %%% This program performs an analysis of step response of Design P-PI -PIDControler[num,den]=tfdata(pidsys(i),'v'); %%% Compute steady valueFinalvalue(i)=polyval(num,0)/polyval(den,0); %%% Compute overshoot[y,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.9*Finalvalue,m=m+1;endRiseTime(i)=t(m)-t(n); %%% Compute settling timer=length(t);while (y(r)>0.98*Finalvalue & y(r)<1.02*Finalvalue)r=r-1;endSettlingTime(i)=t(r) ; %%endFinalvalue,PeakTime,OvershootPercent,RiseTime,SettlingTim e结果：012340.20.40.60.8sysTime (sec)A m p l i t u d eP -sysTime (sec)A m p l i t u d e024680.511.5 P I-sysTime (sec)A m p l i t u d eP ID-sysTime (sec)A m p l i t u d e024680.510123450.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_sys1.00000.895445.73400.21324.9886。