天津大学控制系统设计与仿真上机题

1、控制系统设计与仿真上机实验报告学院：自动化学院班级：自动化姓名：学号：一、 第一次上机任务1、熟悉matlab软件的运行环境，包括命令窗体，workspace等，熟悉绘图命令。2、采用四阶龙格库塔法求如下二阶系统的在幅值为1脉宽为1刺激下响应的数值解。 ，解：解题思路：二阶系统输入的是一个脉宽为1的信号，并非阶跃信号，故输入较特殊。题中给出的是传递函数形式，将其转换为状态空间的形式更容易加入输入。M文件：function t,y = ode4( g,x0,h,v,t0,tf ) g=ss(g); Ab=g.A-g.B*v*g.C;B=g.B;C=g.C;x=x0;y=0;t=t0; %将传递函

2、数转换成状态空间形式 N=round(tf-t0)/h); for i=1:N if(i<=20) r=1; else r=0; end %幅值脉宽为1的输入 k1=Ab * x+B*r; k2=Ab * (x+h*k1/2)+B*r; k3=Ab * (x+h*k2/2)+B*r; k4=Ab * (x+h*k3)+B*r; x=x+h*(k1+2*k2+2*k3+k4)/6; %采用四阶龙格库塔法 y=y,C*x; %输出值 t=t,t(i)+h; end 程序： num=100;den=1,10,100;g=tf(num,den); %传递函数 tf=10; t0=0; h=0.0

3、5; x=0;0; %设置参数 t,y=ode4(g,x,h,0,t0,tf); %调用ode4()函数 plot(t,y); %绘制在幅值脉宽为1的输入下响应 的曲线图响应曲线：3、采用四阶龙格库塔法求高阶系统阶单位跃响应曲线的数值解。 ，解：解题思路：该题与题2类似，区别在于该题分母与系统输入不同（分母成为三阶的系统）。故可沿用上面的方法解答。M文件: function t,y = ode4( g,x0,h,v,t0,tf ) g=ss(g); Ab=g.A-g.B*v*g.C;B=g.B;C=g.C;x=x0;y=0;t=t0; %将传递函数转换成状态空间形式 N=round(tf-t0

4、)/h); r=1; %输入单位阶跃信号 for i=1:N k1=Ab * x+B*r; k2=Ab * (x+h*k1/2)+B*r; k3=Ab * (x+h*k2/2)+B*r; k4=Ab * (x+h*k3)+B*r; x=x+h*(k1+2*k2+2*k3+k4)/6; %采用四阶龙格库塔法 y=y,C*x; %输出值 t=t,t(i)+h; end 程序： num=100;den=conv(1,10,100,5,1);g=tf(num,den); %传递函数 tf=10; t0=0; h=0.05; x=0;0;0; %设置参数 t,y=ode4(g,x,h,0,t0,tf);

5、 %调用ode4()函数 plot(t,y); %绘制在单位阶跃输入下响应的曲线图 响应曲线：4、 自学OED45指令用法，并求解题2中二阶系统的单位阶跃响应。解：解题思路：ODE45是Matlab中龙格库塔专门的指令，可大大简化程序。M文件： function g=ab(t,x) g=x(2);-10*x(2)-100*x(1)+100; %转换为微分方程程序： t,y=ode45('ab',0,10,0 0); %调用ode45（）函数t变化范围0-10，初值为0 plot(t,y); %绘制单位阶跃响应曲线响应曲线：二、第二次上机任务1、试用simulink方法解微分方程

6、，并封装模块，输出为。得到各状态变量的时间序列，以及相平面上的吸引子。参数入口为的值以及的初值。（其中，以及初值分别为） 提示：模块输入是输出量的微分。解：用simulink进行仿真,模块封装图如下：封装内部结构图：各状态变量的时间序列图：（图中振荡的绿、红、蓝色线为三个变量的时间序列图）相平面的吸引子（以x2、x3为例）：2、用simulink搭建PI控制器的控制回路，被控对象传递函数：，分别分析 （1）、比例系数由小到大以及积分时间由小到大对阶跃响应曲线的影响。 （2）、控制器输出有饱和以及反馈有时滞情况下，阶跃响应曲线的变化。 （3）、主控制回路传递函数为：，副回路为：，主回路采用PI控

7、制器，副回路采用P控制器，分析控制系统对主回路以及副回路的阶跃扰动的抑制。注：PI控制器表达式为，串级控制如图所示。解：（1）用simulink进行仿真,仿真图如下:比例系数由小到大：积分时间（Ti）由小到大： 由上面两组图可以总结：比例系数由小到大，调节时间变短，响应变快，但超调量增大；积分时间（Ti）由小到大，超调量减小，但调节时间变长，响应变慢。（2）用simulink进行仿真,仿真图如下:结果： （ 有输出饱和） （有延迟） 由上图可以总结：有控制器输出有饱和的情况下，输出被控制在了饱和值；反馈有时滞情况下，输出出现了振荡。（3）用simulink进行仿真,仿真图如下:结果： 由上图可

8、以总结：对于副回路的阶跃扰动，主回路和副回路都能对其进行抑制；对于主回路的阶跃扰动，主回路能对其进行抑制，而副回路失去控制作用。3、编写S函数模块，实现两路正弦信号的叠加，正弦信号相位差为60度。解：用simulink搭建仿真电路如图： S函数的编写： functionsys,x0,str,ts=addsin(t,x,u,flag) switch flag, case 0 sys,x0,str,ts=mdlInitializeSizes(); %条用mdlInitializeSizes函数进行初始化 %阶段，设置各参数值 case 3 sys=mdlOutputs(u); %定义输出的计算算法

9、 case 1,2,4,9,sys=; otherwise error('ErrFLag=',num2str(flag); %其他情况显示出错 end function sys,x0,str,ts=mdlInitializeSizes() sizes=simsizes; sizes.NumContStates=0; %连续状态的个数 sizes.NumDiscStates=0; %离散状态的个数 sizes.NumOutputs=1; %输出变量的个数 sizes.NumInputs=2; %输入变量的个数 sizes.DirFeedthrough=1; %输入直接传到输出口

10、sizes.NumSampleTimes=1; %有一个采样周期 sys=simsizes(sizes); x0=; str= ; ts=-1 0; function sys=mdlOutputs(u) sys=u(1)+u(2); %输出等于u(1)+u(2)实验结果： Scope+波形 Scope波形4、熟悉SimPowerSystem模块，试用Thyristor Converter模块以及Synchronized 6-pulse generator模块实现三相电整流。（自学内容）解：解题思路：该题内容即三相交流电的整流，改变触发角改变整流波形。题中要求采用一周期6波头的波形。用simul

11、ink搭建仿真电路：遇到问题：试着改变simulation中configuration patrameter中的选项，可以解决上图问题，但打开scope发现结果仍然不对。5、利用触发子系统构建以零阶保持器，实现对正弦信号的采样，并比较不同采用周期下的采样波形。解：用simulink搭建采集仿真电路如图：改变采样周期值，T从小变大，采样结果如下组图： 由上组图可见采样周期越小，原正弦信号采集之后损失的越小，获得的正弦信号也就越接近模拟量。6、若被控对象传递函数为,控制器为，试用simulink搭建一单位反馈控制系统，分析采用周期T对系统单位阶跃响应的影响。解：解题思路：控制器D（z）是离散量，故

12、在搭建系统时要用到采样器和零阶保持器用于A/D的转换中。用simulink搭建系统如图： T分别取0.5,1,5结果显示：7、设一单位反馈控制系统，控制器采用PI控制，Kp=200,Ki=10, 控制器饱和非线性宽度为2，受控对象为时变模型，由微分方程给出，如下： 求系统单位阶跃响应，并分析不同Kp取值对响应曲线的影响。解：用simulink搭建仿真电路如图：Kp从小到大变化，对应的响应曲线图： 由上组图可知，随着Kp的增大，系统响应变快，但受到输出饱和的影响，故Kp一定大时调节时间几乎不变。三、第三次上机任务1、熟悉控制系统各个模型表示方法的命令以及它们之间的相互转化。（展开形式，零极点形式

13、，状态空间形式以及部分分式形式。）解：主要的相互转换的命令有：residue：传递函数模型与部分分式模型互换ss2tf：状态空间模型转换为传递函数模型ss2zp： 状态空间模型转换为零极点增益模型tf2ss： 传递函数模型转换为状态空间模型tf2zp： 传递函数模型转换为零极点增益模型zp2ss： 零极点增益模型转换为状态空间模型zp2tf： 零极点增益模型转换为传递函数模型2、试用至少三种方法，判断一下系统的稳定性: 解：法1：根据极点的位置判断系统稳定性 num=1 2 3 1; den=1 5 2 1 1; z,p,k=tf2zp(num,den); %将传递函数模型转换为零极点增益模型

14、 ii=find(real(p)>0); n1=length(ii); if(n1>0) disp('System is unstable'); disp(p(ii); %在显示不稳定的同时输出不稳定极点 else disp('System is stable'); end pzmap(num,den); %绘制零极点图，更直观的显示显示结果：System is unstable 0.1263 + 0.5642i 0.1263 - 0.5642i法2：根据特征值判断系统的稳定性 den=1 5 2 1 1; P=poly(den); %将分母转换为多

15、项式的形式，即系统的特征方程 r=roots(P); %系统特征方程的跟就是闭环极点 ii=find(real(r)>0); n1=length(ii); if(n1>0) disp('System is unstable'); disp(r(ii); %显示不稳定的同时输出不稳定的根 else disp('System is stable'); end显示结果： System is unstable 5.0000 2.0000 1.0000 1.0000 + 0.0000i 1.0000 - 0.0000i法3：根据部分分式中极点的位置判断系统的稳

16、定性 num=1,2,3,1; den=1,5,2,1,1; r,p,k=residue(num,den); %将传递函数模型转换为部分分式模型 ii=find(real(p)>0); n=length(ii); if(n>0) disp('System is Unstable'); disp(p(ii); %显示不稳定的同时输出不稳定的点 else disp('System is Stable'); end显示结果：System is Unstable 0.1263 + 0.5642i 0.1263 - 0.5642i解：法1：利用李雅普诺夫第二法

17、判断系统的稳定性 A=1,3;5,2; Q=eye(size(A); P=lyap(A,Q); %求解矩阵P ii=find(P(1,1)>0); %判断一阶行列式的值是否大于0 jj=find(det(P)>0); %判断二阶行列式的值是否大于0 w1=length(ii); w2=length(jj); if (w1>0&&w2>0) %w1、w2均大于0说明稳定 disp('the system is stable'); else disp('the system is unstable'); end 显示结果：th

18、e system is unstable法2：解状态方程的特征值并根据特征值判断系统稳定性A=1,3;5,2; r=poly(A); p=roots(r); %利用多项式求解特征值 ii=find(real(p)>0); n1=length(ii); if(n1>0) disp('the unnstable poles are:'); disp(p(ii); %显示不稳定点 else disp('systerm is stable');end显示结果：the unnstable poles are: 5.4051法3：解状态方程的特征值并根据特征值判

19、断系统稳定性，使用另一种方法求解特征值A=1,3;5,2; p=eig(A); %求矩阵A的特征值，利用MATLAB中自带求解方式eig()ii=find(real(p)>0); n1=length(ii);if(n1>0) disp('the unnstable poles are:'); disp(p(ii); %显示不稳定点else disp('systerm is stable');end显示结果：the unnstable poles are: 5.40513、试产生一周期为5秒，时长为30秒，最大值为1，最小值为0的三角波；得到如下一阶系

20、统在三角波输入下的时间响应曲线。 解：num=1;den=2 1;t=0:0.1:30;x=t/5*2*pi;u=0.5*(sawtooth(x,0.5)+1); %调用三角波函数sawtooth（）lsim(num,den,u,t); %绘制时间响应曲线 响应曲线如下图：5、对如下二阶系统做时域分析，得到阻尼比在01之间变化的时候，阶跃响应的上升时间，调节时间，峰值时间，超调量以及衰减比（第一个峰值与稳态值之差与第二个峰值与稳态值之差的比）其中。 解：解题思路：在MATLAB中曲线图实际上是一连串的点绘制而成的。根据这一性质，可以总结出个指标的求解方法 上升时间：有超调时，上升时间=第一次稳

21、态值时间（即响应值=1），从序列开始依次检验是否已等于1；无超调时，上升时间=从稳态值的10%到稳态值的90%所需要的时间，从序列开始分别找稳态值10%和90%的值，两者相减即可。 调节时间：时间序列到达稳定（以正负2%为规定误差）的时间，可以从序列的倒序开始依次进行判断，直到满足稳定 峰值时间：峰值时间即响应达到最大值的时间，故调用max（）函数即可 超调量：在求峰值时间的同时记录下最大值，并用超调量公式求解即可 衰减比：第一组超调量由风值时间可以确定，第二组在第一组之后再找最大值，方法同上。zeta=0.2; tp=0; %峰值时间mo=0; %超调量tr=0; %上升时间ts=0; %调

22、节时间 sj=0; %衰减比 while zeta<=1num=25; %wn=5den=1,10*zeta,25;g=tf(num,den);y,t=step(g);plot(t,y); %绘制图形 grid on ;hold onya,k=max(y); %取最大值，并记录时间（即峰值时间）y1=t(k);tp=tp,y1c=1;y2=100*(ya-c)/c; %超调量公式mo=mo,y2if y2>0n=1;while y(n)<c n=n+1;endy3=t(n);tr=tr,y3 %有超调时，上升时间=第一次到达稳态值所需的时间 else n=1; while y

23、(n)<0.1*c n=n+1; end m=1; while y(m)<0.9*c m=m+1; end y3=t(m)-t(n); tr=tr,y3 %无超调时，上升时间=稳态值的10%到稳态值的90%end e=length(t);while y(e)>0.98*c&&y(e)<1.02*c e=e-1;endy4=t(e);ts=ts,y4 %以正负2%规定误差带，并从时间序列的倒序开始m=1; f=1;while m<3&&f<length(y)if (y(f)>y(f+1)&&y(f-1)&l

24、t;y(f) pm(m)=y(f); m=m+1;endf=f+1;endif pm(1)>c&&pm(2)>cy5=(pm(1)-c)/(pm(2)-c); %衰减比公式sj=sj,y5 endzeta=zeta+0.2;end显示结果：tp = 0 0.6295 0.6820 0.7869 1.0492 1.5633mo = 0 52.5701 25.3786 9.4778 1.5164 -0.3554tr = 0 0.3672 0.4459 0.5596 0.8393 0.6715ts= 0 3.8819 1.6787 1.1716 0.7475 1.1646

25、sj = 0 3.6126 15.5415 5.8041 0.9286 0.9286 阻尼比从0-1变化时的阶跃响应曲线图 6、已知开环传递函数如下，1）试用根轨迹方法得到其临界稳定增益。2）若k=10,试用伯德图方法，判断其稳定性。解：1）figure(1);num=1;den=conv(2,1,conv(1,1,0.1,1);rlocus(num,den);K,poles=rlocfind(num,den); %绘制根轨迹图disp('K=');disp(K);显示结果：Select a point in the graphics windowselected_point

26、= 0.0143 + 3.9706iK= 35.3516 2）figure(2);num=10; %K=10den=conv(2,1,conv(1,1,0.1,1);bode(num,den); %绘制bode图grid; %绘制栅格线，便于判断伯德图： 由伯德图可看出，在对数幅频特性L(w)>0的区段内，对数相频特性曲线对-180线没有穿越，故判断K=10该系统稳定。7、已知系统开环传递函数如下试设计一超前校正环节，使得超调量为20%，调节时间为1s。系统单位斜坡稳态响应误差为10%。并作出校正前后后的系统单位阶跃响应时域曲线加以比较。解：解题思路：设计校正环节常用根轨迹法和频率特性法，分别用于指标给出时域形式和频域特性。该题中给出的控制指标是时域形式，故采用根轨迹法进行校正环节的设计。具体的流程可以归纳为:确定校正器增益 Kc ；校验原系统的阶跃响应超调量是否满足要求；确定期望极点位置（由超调量公式确定，由调节时间公式确定n，再由特征方程求解极点）；求校正器