控制系统仿真实验报告.docx
- 文档编号:4271719
- 上传时间:2022-11-28
- 格式:DOCX
- 页数:27
- 大小:987.18KB
控制系统仿真实验报告.docx
《控制系统仿真实验报告.docx》由会员分享,可在线阅读,更多相关《控制系统仿真实验报告.docx(27页珍藏版)》请在冰豆网上搜索。
控制系统仿真实验报告
控制系统仿真实验报告
专业自动化
班级0807班
学号--------------
姓名------
日期:
2011年12月25日
实验一、MATLAB语言编程
一、实验目的
熟悉MATLAB语言及其环境,掌握编程方法。
要求认真听取实验指导老师讲解与演示。
二、实验内容
1、运行交互式学习软件,学习MATLAB语言。
2、在MATLAB的命令窗口下键入如下命令:
INTRO
(注意:
intro为一个用MATLAB语言编写的幻灯片程序,主要演示常用的MATLAB语句运行的结果。
)
然后,根据显示出来的幻灯片右面按钮进行操作,可按START——NEXT——NEXT按钮一步步运行,观察。
3、自编程序并完成上机编辑,调试,运行,存盘:
(1)用MATLAB命令完成矩阵的各种运算。
A=
11121314
21222324
31323334
41424344
求出下列运算结果
(1)A(:
,1)
(2)A(2:
)(3)A(1:
2,2:
3)(4)A(2:
3,2:
3)(5)A(:
,1:
2)(6)A(2:
3)(7)A(:
)(8)A(:
,:
)(9)ones(2,2)(10)eye
(2)
程序代码:
A=[11:
14;21:
24;31:
34;41:
44];
disp('A=');disp(A);
disp('
(1)A(:
1):
');disp(A(:
1));
disp('
(2)A(2,:
):
');disp(A(2,:
));
disp('(3)A(1:
2,2:
3):
');disp(A(1:
2,2:
3));
disp('(4)A(2:
3,2:
3):
');disp(A(2:
3,2:
3));
disp('(5)A(:
1:
2):
');disp(A(:
1:
2));
disp('(6)A(2:
3):
');disp(A(2:
3));
disp('(7)A(:
):
');disp(A(:
));
disp('(8)A(:
:
):
');disp(A(:
:
));
disp('(9)ones(2,2):
');disp(ones(2,2));
disp('(10)eye
(2):
');disp(eye
(2));
运行结果:
(1)A(:
1):
11
21
31
41
(2)A(2,:
):
21222324
(3)A(1:
2,2:
3):
1213
2223
(4)A(2:
3,2:
3):
2223
3233
(5)A(:
1:
2):
1112
2122
3132
4142
(6)A(2:
3):
2131
(7)A(:
):
11
21
31
41
12
22
32
42
13
23
33
43
14
24
34
44
(8)A(:
:
):
11121314
21222324
31323334
41424344
(9)ones(2,2):
11
11
(10)eye
(2):
10
01
(2)绘制数学函数与矩阵运算功能:
(可自选)
例如:
y(t)=1-2e-tsin(t)(0 理解数组运算与矩阵运算功能。 代码如下: t=0: 0.1: 8; y=1-2*exp (1)-t.*sin(t); figure('name','huizhishuxuemoxing'); plot(t,y); title('y(t)=1-2e-tsin(t)'); xlabel('t'); ylabel('y'); 运行结果: 4、理解函数文件与命令文件的区别,并自编编函数文件并调用。 自定义函数: functions=sqrt(r) s=pi*r^2+r^3 调用函数主代码: r=2; s=sqrt(r) 运行结果: s= 1.4142 5、学会通过Help,熟悉MATLAB中为用户提供的功能各异的函数与命令。 实验二、数值积分算法练习与函数调用 一、实验目的 理解数值积分法,熟练掌握MATLAB的函数调用。 二、实验示例介绍: 1、用Euler法求初值问题的数值解: 设方程如下: du/dt=u-2t/u(2-1)式 u(0)=1t=[0,1] 取步长h=0.1,名为FZSYZ1.M: t0=0; tf=1; x0=1; h1=0.1; t=[t0: h1: tf]; n=length(t) u=x0; uu (1)=u; fori=2: n du=2*t(i-1)/u; u=du*h1+u; uu(i)=u; end figure (1) plot(t,uu) 运行结果: 2、在MATLAB中提供了现成的数值积分法的函数,如ode1,ode23,ode45求解微分方程组,下面介绍ode23它为二阶/三阶的RKF方法在MATLAB的ToolBox文件夹中MATLAB/funfun下的M文件,在此介绍其调用方法与应用例子如下: [t,x]=ode23(‘系统函数名’,to,tf,x0,to1,trace) 其中,系统函数名为描述系统状态方程的M函数的名称,该函数名在调用时应该用引号括起来(文件字符串) t0,tf为起始与终止时间 x0: 系统初始状态变量的值(列向量) to1: 控制解的精度。 (缺省值为t01=10在ode23中) trace: 输出形式控制变量,非零则程序运行的每步都显示出来。 Trace为缺省值——不显示中间结果。 t: 输出参数返回积分的时间离散值(列向量) x: 输出参量,返回每个时间点值的解的列向量 注意: 系统函数的编写格式为固定的 Functionxdot=fun21(t,x) Xdot=x-t^2 并以fun21.m存盘。 第二步: 编写如下程序并以fzsy22.m存盘。 t0=0;tf=3;t01=e-6; x0=1;trace=1; [t,x]=ode23(‘fun21’,t0,tf,x0,t01,trace) Plot(t,x) 第三步: 在命令窗口下运行fzsy21即可求出x的解,并画出曲线。 3、实验具体内容、步骤、要求: (1)运行交互式软件中函数调用,学习程序: (2)试将(2-2)方程改为用Euler编程求解试比较用ode23求解结果。 先编写如下程序: functiondx=fun22(t,x) dx=x-t.^2; 代码如下: t0=0; tf=3; options=odeset('RelTol',1e-6,'AbsTol',1e-6); [T,X]=ode23('fun22',[t0,tf],1,options); figure,plot(T,X); title('ode23'); t0=0; tf=3; x0=1; h1=0.1; t=[t0: h1: tf]; n=length(t) x=x0; xx (1)=x; fori=2: n dx=x-t(i-1)^2; x=dx*h1+x; xx(i)=x; end figure (2) plot(t,xx) 运行结果: (3)试将(2-1)方程改为用ode23算法调用函数求解,并试比较结果。 自定义函数fu233.m functionudot=fun233(t,u) udot=u-2*t/u; end 代码如下: t0=0; tf=1;tol=1e-6; u0=1;trace=1; [t,u]=ode23('fun233',t0,tf,u0,tol,trace); subplot(211) title('ode23') holdon plot(t,u) t0=0; tf=1; x0=1; h1=0.1; t=[t0: h1: tf]; n=length(t) u=x0; uu (1)=u; fori=2: n du=u-2*t(i-1)/u; u=du*h1+u; uu(i)=u; end subplot(212) title('Euler') holdon plot(t,uu) 运行结果: (4)利用ode23或ode45求解线性时不变系统微分方程 自定义函数fun34.m functiondy=fun34(t,y) A=[-0.5,1;-1,-0.5]; dy=A*y; 代码如下: t0=0;tf=4; y0=[01]; options=odeset('RelTol',1e-6,'AbsTol',[1e-61e-6]); [T,Y1]=ode23('fun34',[t0,tf],y0,options); figure('name','name'); plot(T,Y1(: 1),'g*',T,Y1(: 2),'m-'); title('ode23'); xlabel('t'); ylabel('y'); legend('y1','y2'); 运行结果: (5)求出G1(S)=2/(S^2+2S+1)与G2(S)=1/(2S^3+3S^2+1)的单位阶跃响应,并分别求出状态空间模型。 先作题一代码如下: n1=[2]; d1=[121]; step(n1,d1); [A1,B1,C1,D1]=tf2ss(n1,d1) n2=1;d2=[2331]; figure (2) step(n2,d2); [A2,B2,C2,D2]=tf2ss(n2,d2) n=conv(200,[12]); d=conv([11],[11042]); [Z,P,K]=tf2zp(n,d) figure(3) step(n,d); 运行结果: (6)设方程式为y=-40y,y(0)=2用各种数值积分方法与不同步长求方程式的数值解,并比较之。 代码如下: figure('name','Eulerdy=-40y'); t0=0; tf=1; h=0.01; t=t0: h: tf; n=length(t); y0=2; y_Euler (1)=y0; fori=2: n Q(i-1)=-40*y_Euler(i-1)*h; y_Euler(i)=y_Euler(i-1)+Q(i-1); end [T,Y]=ode23('fun36',[t0,tf],y0); plot(T,Y,'m*'); holdon plot(t,y_Euler); title(['h=',num2str(h)]); xlabel('t'); ylabel('y'); 运行结果: 实验三、控制工具箱与SIMULINK软件应用 一、实验目的: 熟悉工具箱及其使用,进行系统仿真分析,通过仿真对系统进行校正校验。 二、实验预习要求: 必须先复习教材及上课介绍的有关控制工具箱命令与SIMULINK仿真工具的使用,并对实验题目作好准备。 三、学会调出、运行已由SIMULINK建立的仿真模型。 在打开MATLAB的窗口下, (1)、输入SIMULINK或双击命令窗口的工具栏左起第二个按钮(NEWSIMULINKModel)将会打开,LibrarySIMULINK然后指向FILE菜单下拉菜单open调出fzsy31.mdl文件。 然后在fzy31.mdl文件的菜单观察并记录有关设置参数,然后指向Start下拉菜单单击一次则观察输出图形。 该仿真模型(SIMULINK系统结构图)如下: (2)按照下图修改上述仿真模型,观察分析系统输出。 四、实验设计题目与要求 解: 对于前3问的实验程序代码为: clc;clearall; figure('name','根轨迹') num1=[1]; den1=[100]; figure (1),rlocus(num1,den1); title('增加比例补偿器'); num2=[11]; den2=conv([100],[15]); figure (2),rlocus(num2,den2); title('增加超前补偿器'); num3=[15]; den3=conv([100],[11]); figure(3),rlocus(num3,den3); title('增加滞后补偿器'); 运行结果: 分析: 通过以上的图像可知, (1)增加比例补偿器后,系统仍处在临界稳定状态,属于不稳定系统,增加Kc的作用仅是增加闭环纯虚数极点的大小,不能影响系统的稳定。 (2)增加超前补偿器后系统的根轨迹都在左半平面内,当Kc≠0时,系统处稳定,在一定的范围内增加Kc的作用是增强系统的稳定性。 (3)增加滞后补偿器后,系统的三条根轨迹有两条都在右半平面,系统处在不稳定,增加Kc就进一步增强系统的不稳定性。 对于第4问,仿真系统图如下: 模块I的仿真详细图 系统总图 分析,模块I、II、III调节器分别对应前三问中的三种情况,并且在每个模块中分别设置Kc=0.1、0.5、1,通过示波器观察Kc值不同时对同一个系统的影响。 仿真图如下: 由上图可知,随着Kc的增大系统的快速性逐渐的增大,但是系统的震荡性越大。 解: 实验才程序如下,图像如下: clc;clearall; z1=[];p1=[0-2];k1=4; G=zpk(z1,p1,k1); %PI补偿器K1(s)=(s+1)/s,K2(s)=(s+0.5)/s K1=tf([11],[10]); K2=tf([10.5],[10]); A=K1*G; sysA=feedback(A,-1) sysB=feedback(G,K1) C=K2*G; sysC=feedback(G,-1) subplot(231),step(sysA)%;title('系统A阶跃响应'); subplot(232),step(sysB);title('系统B阶跃响应'); subplot(233),step(sysC);title('系统C阶跃响应'); t=0: 0.001: 10; u=t;%斜坡输入信号 subplot(234),lsim(sysA,u,t),title('系统A斜坡响应'); subplot(235),lsim(sysB,u,t);title('系统B斜坡响应'); subplot(236),lsim(sysC,u,t);title('系统C斜坡响应'); 运行结果: 分析: 只有系统B是稳定的,只有针对系统B而言才有稳态误差;阶跃响应终值为0,稳态误差为1;斜坡响应终值约为1,稳态误差为无穷大。 取Kp=9,Ki=7,凑试Kd如下: Kp=9;Ki=7; Kd=[5: 13]; figure(3) fori=1: length(Kd) G3=tf([Kd(i),Kp,Ki],[1,0]); sysG3=G1*G3; sysG33=feedback(sysG3,1); step(sysG33); axis([0,13,0.2,1.6]) holdon; 解: 实验程序如下: 取Kp=9,凑试Ki如下: Kp=9; Ki=[1: 2: 10]; figure (2) fori=1: length(Ki) G2=tf([Kp,Ki(i)],[1,0]); sysG2=G2*G1; sysG22=feedback(sysG2,1); step(sysG22); holdon; 首先凑试Kp方法如下: num=[1]; den=conv([1,-1],[1,-1]); G1=tf(num,den); Kp=[1: 3: 15]; fori=1: length(Kp) G=Kp(i)*G1; sysG=tf(G); sysG11=feedback(sysG,1); step(sysG11); holdon 凑试Kp如上图凑试Ki如上图凑试Kd如上图 综合上面的PID参数分别凑试可得出当Kp=9,Ki=7,Kd=12~13之间,校正后的系统阶跃输出如下: 实验四、数字控制系统仿真与综合应用 解,将上面的代码输入后图像如下: 解: (1)建立了仿真系统如下: 仿真系统的仿真如上 仿真示波器的输出如下: (2)先设计连续系统PID调节器,然后将其离散化,建立simulink仿真系统,观察加入零阶保持器后系统的输出,根据指标再将微调一下离散PID参数,使其达到理想指标。 连续的PID凑试法程序如下: Kd凑试如下: Kp=0.4;Ki=0.001; Kd=[0: 0.02: 0.1]; figure(3) fori=1: length(Kd) G3=tf([Kd(i),Kp,Ki],[1,0]); sysG3=G1*G3; sysG33=feedback(sysG3,1); step(sysG33); holdon axis([0,10,0.2,1.6]) end Ki凑试如下: p=0.4; Ki=[0.001: 0.1: 1]; figure (2) fori=1: length(Ki) G2=tf([Kp,Ki(i)],[1,0]); sysG2=G2*G1; sysG22=feedback(sysG2,1); step(sysG22); holdon axis([0,13,0,2]) end Kp凑试如下: clc;closeall;clearall; num=[2000]; den=conv(conv([1,0],[1,10]),[1,20]); G1=tf(num,den); Kp=[0.1: 0.1: 1]; fori=1: length(Kp) G=Kp(i)*G1; sysG=tf(G); sysG11=feedback(sysG,1); step(sysG11); holdon axis([0,13,0.2,1.6]) end 有Kp凑试图像可知,Kp从0.1~1逐渐变大过程中,系统反应速度逐渐增大,但考虑Kp越大导致的超调也会越大,所以取Kp=0.4,超调为20%; 由图像可知随着Ki的增大系统的超调增大速度增快,但系统震荡越大,取Ki=0.001,超调为小于20%; 由Kd凑试图像可知,当Kd=0.02时,系统效果较理想。 建立simulink仿真系统如下: 仿真系统总图 调节器subsystem详细图 综合simulink仿真中加入零阶保持器的仿真结果与该PID凑试m文件的仿真结果可知当存在D时系统经过长时间的有文波调节过程后文波越来越大变得不稳定, 所以只用PI调节,取Kp=0.3,Ki=0.001,Kd=0。 最终的系统仿真图如下: 解: 控制系统仿真图如下: 仿真输出图: 分析: 离散环节的采样周期T1根据以下因素确定: (1)系统的频带宽度; (2)实际采样开关硬件的性能;(3)实现数字控制器的计算程序的执行时间的长短。 对于种种原因采样周期T1比较大,但对系统中德连续部分按照采样周期选择仿真步长T2,将会出现较大的误差,因此T1>T2。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 控制系统 仿真 实验 报告