天津大学控制系统设计与仿真上机题概要.docx
- 文档编号:9735231
- 上传时间:2023-02-06
- 格式:DOCX
- 页数:32
- 大小:437.44KB
天津大学控制系统设计与仿真上机题概要.docx
《天津大学控制系统设计与仿真上机题概要.docx》由会员分享,可在线阅读,更多相关《天津大学控制系统设计与仿真上机题概要.docx(32页珍藏版)》请在冰豆网上搜索。
天津大学控制系统设计与仿真上机题概要
控制系统设计与仿真上机实验报告
学院:
自动化学院
班级:
自动化
姓名:
学号:
一、
第一次上机任务
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;%将传递函数转换成状态空间形式
N=round((tf-t0)/h);
fori=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.05;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)/h);
r=1;%输入单位阶跃信号
fori=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);%调用ode4()函数
plot(t,y);%绘制在单位阶跃输入下响应的曲线图
响应曲线:
4、自学OED45指令用法,并求解题2中二阶系统的单位阶跃响应。
解:
解题思路:
ODE45是Matlab中龙格库塔专门的指令,可大大简化程序。
M文件:
functiong=ab(t,x)
g=[x
(2);-10*x
(2)-100*x
(1)+100];%转换为微分方程
程序:
[t,y]=ode45('ab',[0,10],[00]);%调用ode45()函数t变化范围0-10,初值为0
plot(t,y);%绘制单位阶跃响应曲线
响应曲线:
二、第二次上机任务
1、试用simulink方法解微分方程,并封装模块,输出为
。
得到各状态变量的时间序列,以及相平面上的吸引子。
参数入口为
的值以及
的初值。
(其中
,以及初值分别为
)提示:
模块输入是输出量的微分。
解:
用simulink进行仿真,模块封装图如下:
封装内部结构图:
各状态变量的时间序列图:
(图中振荡的绿、红、蓝色线为三个变量的时间序列图)
相平面的吸引子(以x2、x3为例):
2、用simulink搭建PI控制器的控制回路,被控对象传递函数:
,分别分析
(1)、比例系数由小到大以及积分时间由小到大对阶跃响应曲线的影响。
(2)、控制器输出有饱和以及反馈有时滞情况下,阶跃响应曲线的变化。
(3)、主控制回路传递函数为:
,副回路为:
,主回路采用PI控制器,副回路采用P控制器,分析控制系统对主回路以及副回路的阶跃扰动的抑制。
注:
PI控制器表达式为
,串级控制如图所示。
解:
(1)用simulink进行仿真,仿真图如下:
比例系数由小到大:
积分时间(Ti)由小到大
:
由上面两组图可以总结:
比例系数由小到大,调节时间变短,响应变快,但超调量增大;积分时间(Ti)由小到大,超调量减小,但调节时间变长,响应变慢。
(2)用simulink进行仿真,仿真图如下:
结果:
(有输出饱和)(有延迟)
由上图可以总结:
有控制器输出有饱和的情况下,输出被控制在了饱和值;反馈有时滞情况下,输出出现了振荡。
(3)用simulink进行仿真,仿真图如下:
结果:
由上图可以总结:
对于副回路的阶跃扰动,主回路和副回路都能对其进行抑制;对于主回路的阶跃扰动,主回路能对其进行抑制,而副回路失去控制作用。
3、编写S函数模块,实现两路正弦信号的叠加,正弦信号相位差为60度。
解:
用simulink搭建仿真电路如图:
S函数的编写:
function[sys,x0,str,ts]=addsin(t,x,u,flag)
switchflag,
case0
[sys,x0,str,ts]=mdlInitializeSizes();%条用mdlInitializeSizes函数进行初始化%阶段,设置各参数值
case3
sys=mdlOutputs(u);%定义输出的计算算法
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;%输入直接传到输出口
sizes.NumSampleTimes=1;%有一个采样周期
sys=simsizes(sizes);
x0=[];
str=[];
ts=[-10];
functionsys=mdlOutputs(u)
sys=u
(1)+u
(2);%输出等于u
(1)+u
(2)
实验结果:
Scope+波形Scope波形
4、熟悉SimPowerSystem模块,试用ThyristorConverter模块以及Synchronized6-pulsegenerator模块实现三相电整流。
(自学内容)
解:
解题思路:
该题内容即三相交流电的整流,改变触发角改变整流波形。
题中要求采用一周期6波头的波形。
用simulink搭建仿真电路:
遇到问题:
试着改变simulation中configurationpatrameter中的选项,可以解决上图问题,但打开scope发现结果仍然不对
。
5、利用触发子系统构建以零阶保持器,实现对正弦信号的采样,并比较不同采用周期下的采样波形。
解:
用simulink搭建采集仿真电路如图:
改变采样周期值,T从小变大,采样结果如下组图:
由上组图可见采样周期越小,原正弦信号采集之后损失的越小,获得的正弦信号也就越接近模拟量。
6、若被控对象传递函数为
控制器为
,试用simulink搭建一单位反馈控制系统,分析采用周期T对系统单位阶跃响应的影响。
解:
解题思路:
控制器D(z)是离散量,故在搭建系统时要用到采样器和零阶保持器用于A/D的转换中。
用simulink搭建系统如图:
T分别取0.5,1,5结果显示:
7、设一单位反馈控制系统,控制器采用PI控制,Kp=200,Ki=10,控制器饱和非线性宽度为2,受控对象为时变模型,由微分方程给出,如下:
求系统单位阶跃响应,并分析不同Kp取值对响应曲线的影响。
解:
用simulink搭建仿真电路如图:
Kp从小到大变化,对应的响应曲线图:
由上组图可知,随着Kp的增大,系统响应变快,但受到输出饱和的影响,故Kp一定大时调节时间几乎不变。
三、第三次上机任务
1、熟悉控制系统各个模型表示方法的命令以及它们之间的相互转化。
(展开形式,零极点形式,状态空间形式以及部分分式形式。
)
解:
主要的相互转换的命令有:
residue:
传递函数模型与部分分式模型互换
ss2tf:
状态空间模型转换为传递函数模型
ss2zp:
状态空间模型转换为零极点增益模型
tf2ss:
传递函数模型转换为状态空间模型
tf2zp:
传递函数模型转换为零极点增益模型
zp2ss:
零极点增益模型转换为状态空间模型
zp2tf:
零极点增益模型转换为传递函数模型
2、试用至少三种方法,判断一下系统的稳定性:
解:
法1:
根据极点的位置判断系统稳定性
num=[1231];den=[15211];
[z,p,k]=tf2zp(num,den);%将传递函数模型转换为零极点增益模型
ii=find(real(p)>0);n1=length(ii);
if(n1>0)
disp('Systemisunstable');
disp(p(ii));%在显示不稳定的同时输出不稳定极点
else
disp('Systemisstable');
end
pzmap(num,den);%绘制零极点图,更直观的显示
显示结果:
Systemisunstable
0.1263+0.5642i
0.1263-0.5642i
法2:
根据特征值判断系统的稳定性
den=[15211];P=poly(den);%将分母转换为多项式的形式,即系统的特征方程
r=roots(P);%系统特征方程的跟就是闭环极点
ii=find(real(r)>0);n1=length(ii);
if(n1>0)
disp('Systemisunstable');
disp(r(ii));%显示不稳定的同时输出不稳定的根
else
disp('Systemisstable');
end
显示结果:
Systemisunstable
5.0000
2.0000
1.0000
1.0000+0.0000i
1.0000-0.0000i
法3:
根据部分分式中极点的位置判断系统的稳定性
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('SystemisUnstable');
disp(p(ii));%显示不稳定的同时输出不稳定的点
else
disp('SystemisStable');
end
显示结果:
SystemisUnstable
0.1263+0.5642i
0.1263-0.5642i
解:
法1:
利用李雅普诺夫第二法判断系统的稳定性
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('thesystemisstable');
else
disp('thesystemisunstable');
end
显示结果:
thesystemisunstable
法2:
解状态方程的特征值并根据特征值判断系统稳定性A=[1,3;5,2];
r=poly(A);p=roots(r);%利用多项式求解特征值
ii=find(real(p)>0);n1=length(ii);
if(n1>0)
disp('theunnstablepolesare:
');
disp(p(ii));%显示不稳定点
else
disp('systermisstable');
end
显示结果:
theunnstablepolesare:
5.4051
法3:
解状态方程的特征值并根据特征值判断系统稳定性,使用另一种方法求解特征值
A=[1,3;5,2];p=eig(A);%求矩阵A的特征值,利用MATLAB中自带求解方式eig()
ii=find(real(p)>0);n1=length(ii);
if(n1>0)
disp('theunnstablepolesare:
');
disp(p(ii));%显示不稳定点
else
disp('systermisstable');
end
显示结果:
theunnstablepolesare:
5.4051
3、试产生一周期为5秒,时长为30秒,最大值为1,最小值为0的三角波;得到如下一阶系统在三角波输入下的时间响应曲线。
解:
num=1;
den=[21];
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、对如下二阶系统做时域分析,得到阻尼比在0~1之间变化的时候,阶跃响应的上升时间,调节时间,峰值时间,超调量以及衰减比(第一个峰值与稳态值之差与第二个峰值与稳态值之差的比)其中
。
解:
解题思路:
在MATLAB中曲线图实际上是一连串的点绘制而成的。
根据这一性质,可以总结出个指标的求解方法
1上升时间:
有超调时,上升时间=第一次稳态值时间(即响应值=1),从序列开始依次检验是否已等于1;
无超调时,上升时间=从稳态值的10%到稳态值的90%所需要的时间,从序列开始分别找稳态值10%和90%的值,两者相减即可。
2调节时间:
时间序列到达稳定(以正负2%为规定误差)的时间,可以从序列的倒序开始依次进行判断,直到满足稳定
3峰值时间:
峰值时间即响应达到最大值的时间,故调用max()函数即可
4超调量:
在求峰值时间的同时记录下最大值,并用超调量公式求解即可
5衰减比:
第一组超调量由风值时间可以确定,第二组在第一组之后再找最大值,方法同上。
zeta=0.2;
tp=0;%峰值时间
mo=0;%超调量
tr=0;%上升时间
ts=0;%调节时间
sj=0;%衰减比
whilezeta<=1
num=25;%wn=5
den=[1,10*zeta,25];
g=tf(num,den);
[y,t]=step(g);
plot(t,y);%绘制图形
gridon;holdon
[ya,k]=max(y);%取最大值,并记录时间(即峰值时间)
y1=t(k);
tp=[tp,y1]
c=1;
y2=100*(ya-c)/c;%超调量公式
mo=[mo,y2]
ify2>0
n=1;
whiley(n) n=n+1; end y3=t(n); tr=[tr,y3]%有超调时,上升时间=第一次到达稳态值所需的时间 else n=1; whiley(n)<0.1*c n=n+1; end m=1; whiley(m)<0.9*c m=m+1; end y3=t(m)-t(n); tr=[tr,y3]%无超调时,上升时间=稳态值的10%到稳态值的90% end e=length(t); whiley(e)>0.98*c&&y(e)<1.02*c e=e-1; end y4=t(e); ts=[ts,y4]%以正负2%规定误差带,并从时间序列的倒序开始 m=1;f=1; whilem<3&&f if(y(f)>y(f+1)&&y(f-1) pm(m)=y(f); m=m+1; end f=f+1; end ifpm (1)>c&&pm (2)>c y5=(pm (1)-c)/(pm (2)-c);%衰减比公式 sj=[sj,y5] end zeta=zeta+0.2; end 显示结果: tp=00.62950.68200.78691.04921.5633 mo=052.570125.37869.47781.5164-0.3554 tr=00.36720.44590.55960.83930.6715 ts=03.88191.67871.17160.74751.1646 sj=03.612615.54155.80410.92860.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); 显示结果: Selectapointinthegraphicswindow selected_point= 0.0143+3.9706i K= 35.3516 2) figure (2); num=10;%K=10 den=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,再由特征方程求解极点);Ⅳ求校正器的传递函数;Ⅴ校验校正器。 kc=5;%由稳态响应误差为10%求得校正器增益kc=5 num=200; den=conv(conv([1,0],[12]),[110]); sys1=tf(num,den); sys=feedback(sys1,1); subplot(2,1,1);step(sys,'k')%绘制原系统的阶跃响应,判断是否满足指标 holdon sigma=0.2;%超调量指标20% zeta=((log(1/sigma))^2/((pi)^2+(log(1/sigma))^2))^(1/2);%由超调量公式算出zeta wn=4.4/zeta*1;%由调节时间公式确定振荡频率wn p=[12*zeta*wnwn^2]; r=roots(p)%求出期望极点位置 s_1=r(1,1);nk1=2;dk1=conv(conv([10],[0.51]),[0.11]); ngv=polyval(nk1,s_1);dgv=polyval(dk1,s_1);g=ngv/dgv; zetag=angle(g);mg=abs(g);ms=abs(s_1);zetas=angle(s_1); tz=(sin(zetas)-kc*mg*sin(zetag-zetas))/(kc*mg*ms*sin(zetag)); tp=-(kc*mg*sin(zetas)+sin(zetag+zetas))/(ms*sin(zetag)); nk=[tz,1];dk=[tp,1];Gc=tf(nk,dk)%求出校正装置的传递函数 n1=10;d1=conv(conv([10],[0.51]),[0.11]); s1=tf(n1,d1); sys3=feedback(s1*Gc,1);%校正之后的传递函数 subplot(2,1,2);step(sys3); gtext('校正前'); gtext('校正后'); 显示结果: r= -4.4000+8.5887i -4.4000-8.5887i Transferfunction: 0.4756s+1 ------------- 0.01039s+1
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 天津大学 控制系统 设计 仿真 上机 概要
![提示](https://static.bdocx.com/images/bang_tan.gif)