计算机仿真课程方案最近版独活草著

课程设计任务书
2018～2018学年第 2学期
学生姓名：独活草专业班级：08自动化一班
指导教师：自学成才工作部门：zhbit
作者QQ号：747238952
一、课程设计题目
《控制系统建模、分析、设计和仿真》
二、课程设计内容
<一）《控制系统建模、分析、设计和仿真》课题设计内容
最少拍有波纹控制系统
[0号题] 控制系统建模、分析、设计和仿真
设连续被控对象的实测传递函数为：
用零阶保持器离散化，采样周期取0.1秒，分别设计一单位加速度信号输入时的最少拍有波纹控制器Dy(z>和一单位速度信号输入时的最少拍无波纹控制器Dw(z>。

具体要求见<二）。

<二）《控制系统建模、分析、设计和仿真》课题设计要求及评分标准【共100分】
1、求被控对象传递函数G(s>的MATLAB描述。

<2分）
Matlab输入程序：
num=[888 6216 8880]。

den=[1 11 31 21 0 0]。

Gs=tf(num,den>
运行结果为：
Transfer function:
888 s^2 + 6216 s + 8880
------------------------------
s^5 + 11 s^4 + 31 s^3 + 21 s^2
2、求被控对象脉冲传递函数G(z>。

<4分）
Matlab输入程序：
Gz=c2d(Gs,0.1，‘zoh’>
运行结果为：
Transfer function:
0.1347 z^4 + 0.2989 z^3 - 0.5225 z^2 + 0.08645 z + 0.05483
----------------------------------------------------------
z^5 - 4.142 z^4 + 6.772 z^3 - 5.45 z^2 + 2.153 z - 0.3329 sampling time: 0.1
3、转换G(z>为零极点增益模型并按形式排列。

<2分）
[ab k]=zpkdata(Gz>。

Gz=zpk(a,b,k,0.1,'variable','z^-1'>
zero/pole/gain:
0.13473 z^-1 (1+3.403z^-1> (1-0.8187z^-1> (1-0.6065z^-1> (1+0.2408z^-1>
-----------------------------------------------------------------------
(1-z^-1>^2 (1-0.9048z^-1> (1-0.7408z^-1> (1-0.4966z^-1>
sampling time: 0.1
4、确定误差脉冲传递函数Ge(z>形式,满足单位加速度信号输入时闭环
稳态误差为零和实际闭环系统稳定的要求。

<6分）
5、确定闭环脉冲传递函数Gc(z>形式，满足控制器Dy(z>可实现、最少
拍和实际闭环系统稳定的要求。

<8分）
6、根据4、5、列写方程组，求解Gc(z>和Ge(z>中的待定系数并最终
求解Gc(z>和Ge(z> 。

<12分）
7、求针对单位加速度信号输入的最少拍有波纹控制器Dy(z>并说明
Dy(z>的可实现性。

<3分）
分析：<1）选择误差脉冲传递函数Ge(z>，按单位加速度信号输入，Ge(z>中应含有(1-z^-1>^3因子；所以设Gez= (1-z^-1>^3*(Z>
<2）选择闭环函数Gc(z>，因为Gz分子有z^-1的零极点，又有不稳定的零点因子(1+3.403*z^-1>；所以设
Gcz= z^-1 *(1+3.403*z^-1>*<z）
根据Gez中(1-z^-1>^3因子为3阶，可知<z）中含有三个未知数，因此设<z）=h0+h1*z^-1+h2*z^-2,故
Gcz=z^-1*(1+3.403*z^-1>*<h0+h1*z^-1+h2*z^-2）<3）根据Gez与Gcz的最高阶数应该保持一致，故可设
(Z>=<e0+e1*z^-1）,则有Gez= (1-z^-1>^3*<e0+e1*z^-1）
<4）且有Gez=1-Gcz
故只要求解出未知数h0 h1 h2 e0 e1，即可得出误差脉冲传递函数Gez形式与闭环函数Gcz的表达式。

列方程思路：
1.对于Gez，当z=1时，Gez=0，Gcz=1-Gez=1，所以当
z=1,Gcz-1=0。

2.对于Gez的一阶导数diff(Gez,1>，当z=1
时,(diff(Gez,1>,z,1>=0,(diff(Gcz,1>,z,1>=(diff(Gez,1
>,z,1>=0。

3.对于Gez的二阶导数diff(Gez,2>，当z=1
时,(diff(Gez,2>,z,1>=0,(diff(Gcz,2>,z,1>=(diff(Gez,2>,z
,1>=0。

4.对于Gcz,当z=-3.403时,Gcz=0，所以当z=-3.403，Gez-
1=0。

5.对于Gcz,当z=inf时,Gcz=0，所以当z=inf,Gez-1=0。

建立test1.m文件：
syms z e0 e1 h0 h1 h2
Gz=0.13473*z^-1*(1+3.403*z^-1>*(1-0.8187*z^-1>*(1-0.6065*z^-1>*(1+0.2408*z^-1>/(1-z^-1>^2/(1-0.9048*z^-1>/(1-0.7408*z^-1>/(1-0.4966*z^-1>
Gz =
13473/100000/z*(1+3403/1000/z>*(1-8187/10000/z>*(1-
1213/2000/z>*(1+301/1250/z>/(1-1/z>^2/(1-1131/1250/z>/(1-
463/625/z>/(1-2483/5000/z>
Gcz=z^-1*(1+3.403*z^-1>*(h0+h1*z^-1+h2*z^-2>
f1=subs(Gcz,z,1>-1
f2=subs(diff(Gcz,1>,z,1>
f3=subs(diff(Gcz,2>,z,1>
[h0j,h1j,h2j]=solve(f1,f2,f3>
A=double([h0j h1j h2j]>
Gcz=subs(Gcz,[h0 h1 h2],A>
Gez=(1-z^-1>^3*(e0+e1*z^-1>
f4=subs(Gez,z,-3.403>-1
f5=subs(Gez,z,inf>-1
[e0j,e1j]=solve(f4,f5>
B=double([e0 e1]>
Gez=subs(Gez,[e0 e1],B>
运行结果：
Gcz =1/z*(1+3403/1000/z>*(h0+h1/z+h2/z^2>
f1 =4403/1000*h0+4403/1000*h1+4403/1000*h2-1
f2 =-3903/500*h0-12209/1000*h1-4153/250*h2
f3 =11209/500*h0+11709/250*h1+4003/50*h2
h0j =99706454000/853********
h1j =-126270272000/853********
h2j =45950227000/853********
A =1.1681 -1.4793 0.5383
Gcz=1/z*(1+3403/1000/z>*(2630310348345151/2251799813685248-3331078278342513/2251799813685248/z+2424383833506079/450359962 7370496/z^2>
Gez =(1-1/z>^3*(e0+e1/z>
f4=853********/39408131827*e0-85358358827000/134105872607281*e1-1
f5 =e0-1
e0j =1
e1j =156368622481/853********
B = 1.0000 1.8319
Gez=(1-1/z>^3*(1+4125089092710593/2251799813685248/z>
Dyz=Gcz/Gez/Gz
Dyz =
100000/13473*(2630310348345151/2251799813685248-3331078278342513/2251799813685248/z+2424383833506079/450359962 7370496/z^2>/(1-
1/z>/(1+4125089092710593/2251799813685248/z>/(1-
8187/10000/z>/(1-1213/2000/z>/(1+301/1250/z>*(1-
1131/1250/z>*(1-463/625/z>*(1-2483/5000/z>
8、用程序仿真方法分析加速度信号输入时闭环系统动态性能和稳态性能。

<7分）
闭环传递函数：Hz=Dyz*Gz/(1+Dyz*Gz>
Hz =
1/z*(1+3403/1000/z>*(2630310348345151/2251799813685248-3331078278342513/2251799813685248/z+2424383833506079/45035 99627370496/z^2>/(1-
1/z>^3/(1+4125089092710593/2251799813685248/z>/(1+1/z*(1+3 403/1000/z>*(2630310348345151/2251799813685248-
3331078278342513/2251799813685248/z+2424383833506079/45035
99627370496/z^2>/(1-
1/z>^3/(1+4125089092710593/2251799813685248/z>>
simplify(Hz>。

[N,D]=numden(Hz>。

num1=sym2poly(N>。

den1=sym2poly(D>。

t=0:0.1:2。

u=0.5*(t.^2>。

dlsim(num1,den1,u>
软件仿真波形为：
9、用图形仿真方法(Simulink>分析单位加速度信号输入时闭环系统动
态性能和稳态性能。

<8分）
simple(Dyz>
combine(trig>:
(6575775870862877500000000000000*z^5-
22414322766418738680500000000000*z^4+306514098645577368355 020*********z^3-
21067987329182255412020959360000*z^2+727972881814105348618 6719680000*z-
1008722418934002773470566720000>/(758462472244533657600000 000000*z^5-267352290718713834436440000000*z^4-
2020441993717059232870790736000*z^3+1833093331500834295591 733264574*z^2-137630882650940498926607155515*z-
166130636658654386957895373059>
10、确定误差脉冲传递函数Ge(z>形式,满足单位速度信号输入时闭环
稳态误差为零和实际闭环系统稳定的要求。

<6分）
11、确定闭环脉冲传递函数Gc(z>形式，满足控制器Dw(z>可实现、无
波纹、最少拍和实际闭环系统稳定的要求。

<8分）
12、根据10、11、列写方程组，求解Gc(z>和Ge(z>中的待定系数并最
终求解Gc(z>和Ge(z> 。

<12分）
13、求针对单位速度信号输入的最少拍无波纹控制器Dw(z>并说明
Dw(z>的可实现性。

<3分）
Gz=
0.13473 z^-1 (1+3.403z^-1> (1-0.8187z^-1> (1-0.6065z^-1> (1+0.2408z^-1>
-----------------------------------------------------------------------
(1-z^-1>^2 (1-0.9048z^-1> (1-0.7408z^-1> (1-0.4966z^-1>
分析：<1）选择误差脉冲传递函数Ge(z>，按单位速度信号输入，故Ge(z>中应含有(1-z^-1>^2因子，Gz无不稳定的极点，所以
设Gez= (1-z^-1>^2*(Z>
<2）选择闭环函数Gc(z>，因为Gz有z^-1的因子和零点z=-
3.403,z=0.8187,z=0.6065,z=-0.2408，所以取Gcz=a*z^-
1*(1+3.403*z^-1>*(1-0.8187*z^-1>*(1-0.6065*z^-
1>*(1+0.2408*z^-1>
<3）根据Gez与Gcz的最高阶数应该保持一致，故可设
(Z>=<1+t1*z^-1+t2*z^-2+t3*z^-3）,
则有Gez= (1-z^-1>^2*<1+t1*z^-1+t2*z^-2+t3*z^-3）
<4）且有Gez=1-Gcz
故只要求解出未知数 a t1 t2 t3 ,即可得出误差脉冲传递函数Gez形式与闭环函数Gcz的形式。

列方程思路：
1,对于Gez的一阶导数diff(Gez,1>，当z=1时,(diff(Gez,1>,z,1>=0,(diff(Gcz,1>,z,1>=(diff(Gez,1>,z,1>=0。

2.对于Gcz,当z=-
3.403时,Gcz=0,所以当z=-3.403,Gez-1=0。

3.对于Gcz,当z=0.8187时,Gcz=0,所以当z=0.8187,Gez-1=0。

4.对于Gcz,当z=0.6065时,Gcz=0,所以当z=0.6065,Gez-1=0。

建立test2.m文件：
syms z at1 t2t3
Gz=0.13473*z^-1*(1+3.403*z^-1>*(1-0.8187*z^-1>*(1-0.6065*z^-1>*(1+0.2408*z^-1>/(1-z^-1>^2/(1-0.9048*z^-1>/(1-0.7408*z^-1>/(1-0.4966*z^-1>
Gz =
13473/100000/z*(1+3403/1000/z>*(1-8187/10000/z>*(1-
1213/2000/z>*(1+301/1250/z>/(1-1/z>^2/(1-1131/1250/z>/(1-
463/625/z>/(1-2483/5000/z>
Gcz=a*z^-1*(1+3.403*z^-1>*(1-0.8187*z^-1>*(1-0.6065*z^-
1>*(1+0.2408*z^-1>
Gcz=a/z*(1+3403/1000/z>*(1-8187/10000/z>*(1-
1213/2000/z>*(1+301/1250/z>
f1=subs(Gcz,z,1>-1。

[aj]=solve(f1>
aj =
25000000000000/9743904521043
Gcz=subs(Gcz,a,aj>
Gcz=25000000000000/9743904521043/z*(1+3403/1000/z>*(1-
8187/10000/z>*(1-1213/2000/z>*(1+301/1250/z>
Gez=(1-z^-1>^2*<1+t1*z^-1+t2*z^-2+t3*z^-3）。

f2=subs(Gez,z,-3.403>-1。

f3=subs(Gez,z,0.8187>-1。

f4=subs(Gez,z,0.6065>-1。

[t1j,t2j,t3j]=solve(f2,f3,f4>
t1j =1897977623066786822/152385161533393411
t2j =673755091791454319/21769308790484773
t3j =-3511818139920991356/152385161533393411
Gez=subs(Gez,[t1,t2,t3],[t1j,t2j,t3j]>
Gez=(1-
1/z>^2*(1+1897977623066786822/152385161533393411/z+6737550 91791454319/21769308790484773/z^2-
3511818139920991356/152385161533393411/z^3>
Dwz=Gcz/Gez/Gz
Dwz=1504546253029590625/79006507525472256/(1+3505808583160501/ 281474976710656/z+4355792841102475/140737488355328/z^2-
3243389708025163/140737488355328/z^3>*(1-1131/1250/z>*(1-
463/625/z>*(1-2483/5000/z>
>>
14、用程序仿真方法分析单位速度信号输入时闭环系统动态性能和稳态
性能。

<7分）
Tz=Dwz*Gz/(1+Dwz*Gz>
Tz =
14443644029084070363/5629499534213120000/(1+3505808583160501/2 81474976710656/z+4355792841102475/140737488355328/z^2-
3243389708025163/140737488355328/z^3>/z*(1+3403/1000/z>*(1 -8187/10000/z>*(1-1213/2000/z>*(1+301/1250/z>/(1-
1/z>^2/(1+14443644029084070363/5629499534213120000/(1+3505 808583160501/281474976710656/z+4355792841102475/1407374883 55328/z^2-
3243389708025163/140737488355328/z^3>/z*(1+3403/1000/z>*(1 -8187/10000/z>*(1-1213/2000/z>*(1+301/1250/z>/(1-1/z>^2>
求出Gz的零极点形式的表达式：
sysy=zpk([-3.4033 0.81873 0.60653 -0.24078],[1 1 0.90484
0.74082 0.49659],[0.13473],0.1>
Zero/pole/gain:
0.13473 (z+3.403> (z-0.8187> (z-0.6065> (z+0.2408>
--------------------------------------------------
(z-1>^2 (z-0.9048> (z-0.7408> (z-0.4966>
Sampling time: 0.1
将离散的时间模型转换为连续的时间模型：
>> sysyc=d2c(sysy,'tustin'>
Zero/pole/gain:
0.036184 (s-36.64> (s-20> (s+32.69> (s+4.898> (s+1.993>
----------
----------
----------
----------
----------
-----
s^2
(s+0.9991>
(s+2.978>
(s+6.727>
hold on
td=0:0.1:2；
ud=td；
lsim(sysy,ud,td>；
t=0:0.008:2；
u=t
lsim(sysyc,u,t>
软件仿真波形为：
15、用图形仿真方法(Simulink>分析单位速度信号输入时闭环系统动态
性能和稳态性能。

<8分）
simple(Dwz>
combine(trig>:
=(4187670488280980468750000*z^3-
8970827719995516360156250*z^2+6229079448044583452440625*z-1393903736283674639102469>/(219902325555200000000000*z^3+2 738912955594141406250000*z^2+6805926314222617187500000*z-5067796418789317187500000>
Dwz=(4187670488280980468750000-8970827719995516360156250*z^-1+6229079448044583452440625*z^-2-
1393903736283674639102469>*z^-
3/(219902325555200000000000+2738912955594141406250000*z^-1+6805926314222617187500000*z^-2-
5067796418789317187500000>*z^-3
16、根据8、9、14、15、的分析，说明有波纹和无波纹的差别和物理
意义。

<4分）
三、进度安排
6月13至6月14：下达课程设计任务书；复习控制理论和计算机仿真知
识，收集资料、熟悉仿真工具；确定设计方案和步骤。

6月14至6月16：编程练习，程序设计；仿真调试，图形仿真参数整定；
总结整理设计、仿真结果，撰写课程设计说明书。

6月16至6月17：完成程序仿真调试和图形仿真调试；完成课程设计说明
书；课程设计答辩总结。

四、基本要求
1．学生应按照课程设计任务书的要求独立分析、解决问题，按计划完成课程设计任务；
2．不得抄袭或找人代做，否则按考试作弊处理；
3. 学生在完成课程设计时须提交不少于3000字课程设计说明书；说明书
结构为：
<1）封面，<2）任务书，<3）摘要，<4）关键词，<5）目录，<6）正
文，<7）参考文献；
教研室主任签名：
年月日。