天津大学控制系统设计与仿真上机题概要Word下载.docx
《天津大学控制系统设计与仿真上机题概要Word下载.docx》由会员分享，可在线阅读，更多相关《天津大学控制系统设计与仿真上机题概要Word下载.docx（32页珍藏版）》请在冰豆网上搜索。
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）
r=1;
%输入单位阶跃信号
k1=Ab*x+B*r;
%采用四阶龙格库塔法
%输出值
den=conv（[1,10,100],[5,1]）;
0;
%绘制在单位阶跃输入下响应的曲线图
4、自学OED45指令用法，并求解题2中二阶系统的单位阶跃响应。
ODE45是Matlab中龙格库塔专门的指令，可大大简化程序。
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
%绘制单位阶跃响应曲线
二、第二次上机任务
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取值对响应曲线的影响。
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
Systemisstable'
pzmap（num,den）;
%绘制零极点图，更直观的显示
显示结果：
Systemisunstable
0.1263+0.5642i
0.1263-0.5642i
法2：
根据特征值判断系统的稳定性
P=poly（den）;
%将分母转换为多项式的形式，即系统的特征方程
r=roots（P）;
%系统特征方程的跟就是闭环极点
ii=find（real（r）>
disp（r（ii））;
%显示不稳定的同时输出不稳定的根
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）;
%将传递函数模型转换为部分分式模型
n=length（ii）;
if（n>
SystemisUnstable'
%显示不稳定的同时输出不稳定的点
SystemisStable'
SystemisUnstable
利用李雅普诺夫第二法判断系统的稳定性
A=[1,3;
5,2];
Q=eye（size（A））;
P=lyap（A,Q）;
%求解矩阵P
ii=find（P（1,1）>
%判断一阶行列式的值是否大于0
jj=find（det（P）>
%判断二阶行列式的值是否大于0
w1=length（ii）;
w2=length（jj）;
if（w1>
0&
&
w2>
0）%w1、w2均大于0说明稳定
thesystemisstable'
thesystemisunstable'
thesystemisunstable
解状态方程的特征值并根据特征值判断系统稳定性A=[1,3;
r=poly（A）;
p=roots（r）;
%利用多项式求解特征值
theunnstablepolesare:
'
%显示不稳定点
else
disp（'
systermisstable'
end
5.4051
解状态方程的特征值并根据特征值判断系统稳定性，使用另一种方法求解特征值
A=[1,3;
p=eig（A）;
%求矩阵A的特征值，利用MATLAB中自带求解方式eig（）
ii=find（real（p）>
if（n1>
else
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];
[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>
n=1;
whiley（n）<
c
n=n+1;
y3=t（n）;
tr=[tr,y3]%有超调时，上升时间=第一次到达稳态值所需的时间
n=1;
whiley（n）<
0.1*c
n=n+1;
m=1;
whiley（m）<
0.9*c
m=m+1;
y3=t（m）-t（n）;
tr=[tr,y3]%无超调时，上升时间=稳态值的10%到稳态值的90%
e=length（t）;
whiley（e）>
0.98*c&
y（e）<
1.02*c
e=e-1;
y4=t（e）;
ts=[ts,y4]%以正负2%规定误差带，并从时间序列的倒序开始
m=1;
f=1;
whilem<
3&
f<
length（y）
if（y（f）>
y（f+1）&
y（f-1）<
y（f））
pm（m）=y（f）;
m=m+1;
f=f+1;
ifpm
（1）>
c&
pm
（2）>
y5=（pm
（1）-c）/（pm
（2）-c）;
%衰减比公式
sj=[sj,y5]
zeta=zeta+0.2;
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）;
den=conv（[2,1],conv（[1,1],[0.1,1]））;
rlocus（num,den）;
[K,poles]=rlocfind（num,den）;
%绘制根轨迹图
K='
disp（K）;
Selectapointinthegraphicswindow
selected_point=
0.0143+3.9706i
K=
35.3516
2）
figure
（2）;
num=10;
%K=10
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'
）%绘制原系统的阶跃响应，判断是否满足指标
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（'
校正前'
校正后'
r=
-4.4000+8.5887i
-4.4000-8.5887i
Transferfunction:
0.4756s+1
-------------
0.01039s+1
升级会员 