控制系统的采样与仿真

1、 前言 连续系统数字仿真方法是指用数字计算机对连续系统进行仿真的方法。采用这 种方法时首先将连续系统的数学模型转变为适合在数字计算机上进行试验的仿真 模型，实现这种转变的计算方法主要有微分方程数值解法和离散相似法。 MATLAB 产品家族是美国 MathWorks 公司开发的用于概念设计、算法开发、 建模仿真、实时实现的理想的集成环境。是矩阵实验室（Matrix Laboratory）的 简称，是美国 MathWorks 公司出品的商业数学软件，用于算法开发、数据可视化、 数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括 MATLAB 和 SIMULINK 两大部分。 MATLAB

2、 由于其完整的专业体系和先进的设计开发思路，使得 MATLAB 在多种 领域都有广阔的应用空间，特别是在科学计算、建模仿真以及系统工程的设计开 发上已经成为行业内的首选设计工具，它将数值分析、矩阵计算、科学数据可视 化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗 环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供 了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如 C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。 MATLAB 软件工具在自动化专业、测控技术与仪器和电气工程及其自动化等专 业的本科生学习中，经

3、常用来计算、仿真和设计，尤其是 MATLAB 软件的仿真功能， 能使我们对所学知识有更加深入的理解和分析。 目录目录 一、系统分析：一、系统分析：.3 1. 绘制碾磨控制系统开环根轨迹图、绘制碾磨控制系统开环根轨迹图、BODE 图和奈奎斯特图，并判断稳定性；图和奈奎斯特图，并判断稳定性；.4 2.当控制器为当控制器为 () ( ) () c K sa G s sb ，试设计一个能满足要求的控制器（要求用根轨迹法，试设计一个能满足要求的控制器（要求用根轨迹法 和频率响应法进行设计）和频率响应法进行设计） ；.7 21 进行根轨迹校正： .8 2.2 频率校正： .12 3.将采样周期取为将采样周

4、期取为0.02Ts，试确定与，试确定与( ) c Gs对应的数字控制器对应的数字控制器( ) c Gz（要求用多（要求用多 种方法进行离散化，并进行性能比较）种方法进行离散化，并进行性能比较） ；.17 4、5、6：连续，离散单位阶跃输入响应比较：连续，离散单位阶跃输入响应比较.25 6.比较并讨论比较并讨论 4 和和 5 的仿真结果；的仿真结果；.27 7.讨论采样周期的不同选择对系统控制性能的影响；讨论采样周期的不同选择对系统控制性能的影响；.28 8.如控制器改为如控制器改为 PID 控制器，请确定满足性能指标的控制器，请确定满足性能指标的 PID 控制器参数。控制器参数。.33 8.1

5、 用最优 PID 控制法设计：.33 9.如希望尽可能的缩短系统的调节时间，请设计相对应的最小拍控制器，并画出校正如希望尽可能的缩短系统的调节时间，请设计相对应的最小拍控制器，并画出校正 后系统的阶跃响应曲线。后系统的阶跃响应曲线。.37 9.1 采用无纹波最少拍系统设计 .37 10.课程设计总结：课程设计总结：.39 11.参考资料：参考资料：40 一、系统分析：一、系统分析： 某工业碾磨系统的开环传递函数为 10 ( ) (5) G s s s 要求用数字控制器 D(z)来改善系统的性能，使得相角裕度大于45 o ，调节时间小 于 1s(2%准则) 1. 绘制碾磨控制系统开环根轨迹图、B

6、ode 图和奈奎斯特图，并判断稳定性； 2.当控制器为 () ( ) () c K sa G s sb ，试设计一个能满足要求的控制器（要求用根轨 迹法和频率响应法进行设计） ； 3.将采样周期取为 0.02Ts ，试确定与 ( ) c Gs 对应的数字控制器 ( ) c Gz （要求 用多种方法进行离散化，并进行性能比较） ； 4.仿真计算连续闭环系统对单位阶跃输入的响应； 5.仿真计算数据采样系统对单位阶跃输入的响应； 6.比较并讨论 4 和 5 的仿真结果； 7.讨论采样周期的不同选择对系统控制性能的影响； 8.如控制器改为 PID 控制器，请确定满足性能指标的 PID 控制器参数。 9

7、.如希望尽可能的缩短系统的调节时间，请设计相对应的最小拍控制器，并画出 校正后系统的阶跃响应曲线。 1. 绘制碾磨控制系统开环根轨迹图、绘制碾磨控制系统开环根轨迹图、Bode 图和奈奎斯图和奈奎斯 特图，并判断稳定性特图，并判断稳定性： G=zpk(,0 -5,10);sisotool(G);margin(G); 根轨迹图 Bode 图：截止频率为 1.88rad/s,相角裕度为 69 N=0;R=0;Z=P-R=0;该系统稳定。 2.当控制器为当控制器为，试设计一个能满足要求的控，试设计一个能满足要求的控 () ( ) () c K sa G s sb 制器（要求用根轨迹法和频率响应法进行设

8、计）制器（要求用根轨迹法和频率响应法进行设计）： 调节前 Gs=tf(10,1 5 0); Close_S=feedback(Gs,1); Step(Close_S,b); hold on 设计前调节时间为 1.18s 设计前截止频率为 1.88rad/s,相角裕度为 69（第一问中） 21 进行根轨迹校正： 2 2 1,2 =70 2 =0.8 12 4.4 2.55/.25/ 153.75 snn n nn arctg twrad swrad s w pwjwj 取度 由知取 ，求得5，取6 要求的主导极点为 要使得根轨迹向左转，要加入零点。考虑到校正装置的物理 可实现性，加入超前校正装置

9、。 1 1111 11 11 a ( ) b (a) ( ) (2)(b) , a2b180 40 50 c g o o oo o o c s G s s Ks G s s ss p pppp pp pp K（） （） 开环传递函数为 为了使得根轨迹通过根据相角条件 （）（0）（）（） 求得 （0）140，（2）90 （a）（b） 超前装置提供的超前相角为 1 1 111 111 53.75 1 6.512 1 0+511.499 0+511.499 100 6.512 =10 g gpj p K ppp ppp K p K 根据根轨迹的幅值条件 系统的开环增益为 3 3 3 6.512 (

10、) 11.499 6.512 6.499 c c c s G s s z p p zp 10（） 所以 （） 加校正装置后，除要求的主导极点，还有一个闭环零点和 一个非主导极点。 根据(-5+j 3. 75)+(-5-j 3. 75)+ 0+(-5)+(-11. 499)-第八法则 、对系统的影响，例如超调量可能会变大等，但闭环系统的性 能主要由复数极点确定。 ( )100(6.512) ( ) ( )(53.75)(53.75)(6.499) 1 ( )( ) ( )( ) C ss s R ssjsjs C ss R ss s 加校正装置后，系统的闭环传递函数为 系统的单位阶跃响应为 =

11、检验性能： Ds=tf(10*1 6.512,1 11.499); Gs=tf(10,1 5 0); Close_S=feedback(Ds*Gs,1); Step(Close_S,b); hold on 调节时间为 0.863s，符合要求。 G=zpk(-6.512,0,-5,-11.499,10); margin(G); G=zpk(-6.512,0,-5,-11.499,100); margin(G); 相角裕度为 48，符合要求。 2.2 频率校正： 详细设计 要求：静态速度误差为 20，相角裕度不小于 45，调节时间小于 1s（2%） 。 A. 根据静态误差指标确定开环增益 1 00

12、 110 lim( ) ( )lim220 1(5) 10 vc ss Ts KsG s G ss KK Tss s K B. 据确定的增益 K ，画出如下增益经调整后的未校正系统的 Bode 图 G=zpk(,0,-5,100); margin(G); 校正前的相角裕度为 28 C. 计算为满足设计要求所需增加的相位超前角度 从图可知为满足设计要求，还须 25 度左右的超前相角。即令 25 m D. 计算 1 sin 1 sin 2.4638 m m E. 选定最大超前角发生频率 因为校正环节在最大超前相角处有 10log a 的幅值提升，所以把 m m 10log10log(2.4638)

13、3.916() =12rad / dB s 处选为： F.据式 T m 1 计算超前环节的时间常数因子 T 和校正环节的交接频率 1 122.463818.36 pm T 1 7.645 z T H. 对以上设计所得进行检验，看是否满足设计要求。)()(sGsGc /111010 ( ) ( ) 1(5)/1(5) z c p sTsK G s G sK Tss sss s (7.65) ( )24.638 (18.836) c s G s s I性能验证： Ds=tf(24.638*1 7.65,1 18.836); Gs=tf(10,1 5 0); Close_S=feedback(Ds*

14、Gs,1); Step(Close_S,b); hold on 调节时间为 0.573，满足设计要求。 G=zpk(-7.65,0,-5,-18.836,240.638); margin(G); 相角裕度为 48，满足设计要求。 3.将采样周期取为将采样周期取为，试确定与，试确定与对应的数字控对应的数字控0.02Ts( ) c Gs 制器制器（要求用多种方法进行离散化，并进行性能比（要求用多种方法进行离散化，并进行性能比( ) c Gz 较）较）： A选用根轨迹所得到的控制器函数： 10(6.512) 11.499 s Gs s B采用脉冲响应不变法，零阶保持器法，一阶保持器法，双线性变涣 法

15、，零极点匹配方法确定数字控制器 Gc（z） ； Gc=zpk(-6.512,-11.499,10); Gimp=c2d(Gc,0.02,imp)%脉冲响应不变法 Gzoh=c2d(Gc,0.02,zoh) %零阶保持器 Gfod=c2d(Gc,0.02,fod) %一阶保持器 Gtustin=c2d(Gc,0.02,tustin) %双线性 Gmatched=c2d(Gc,0.02,matched) %零极点匹配方 法 性能比较： G0=zpk(,0 -5,10);Gc=zpk(-6.512,-11.499,10); G=series(G0,Gc); G1=c2d(G,0.02,zoh);%零

16、阶保持器 G2=c2d(G,0.02,fod);%一阶保持器 G3=c2d(G,0.02,tustin); %双线性 G4=c2d(G,0.02,matched);%零极点匹配方法 G5=c2d(G,0.02,imp);%脉冲响应不变法 Gk1=feedback(G1,1);Gk2=feedback(G2,1); Gk3=feedback(G3,1);Gk4=feedback(G4,1); Gk5= feedback(G5,1); figure;margin(G1);grid figure;margin(G2);grid figure;margin(G3);grid figure;margin

17、(G4);grid figure;margin(G5);grid figure;step(Gk1,Gk2,Gk3,Gk4,Gk5);legend(zoh,fod,tustin, matched,imp);grid 零阶保持器 一阶保持器 双线性 零极点配置法 脉冲响应不变法 阶跃响应 各种离散化方法的动态性能比较： 离散化方 法/() /(rad/s)%tp/sts/s Zoh457.78250.3640.88 Fod497.76200.340.86 Tustin 487.79200.340.86 Matched438.15260.3750.86 imp 877.4830.060.82 c 4

18、、5、6：连续，离散单位阶跃输入响应比较：连续，离散单位阶跃输入响应比较 %连续系统的阶跃响应 Ds=tf(10*1 6.512,1 11.499); Ghs=tf(100,1 100);%保持器采用一节惯性环节 Gs=tf(10,1 5 0); Close_S=feedback(Ds*Ghs*Gs,1); Step(Close_S,b); hold on %离散系统的阶跃响应 Ts=0.02;i=100; Dz=c2d(Ds,Ts,tustin);%双线性变换 Gz=c2d(Gs,Ts,zoh);%零阶保持器 Close_Z=minreal(feedback(Dz*Gz,1); Y=dste

19、p(Close_Z.num1,Close_Z.den1,i); plot(Ts*(1:i),Y,-.r); hold off 蓝色为连续系统，红色为离散系统。 6.比较并讨论比较并讨论 4 和和 5 的仿真结果的仿真结果: 绘制连续闭环系统和数据采样系统对单位阶跃输入的响应的曲线图；根据曲 线图比较两系统的超调量、调节时间和峰值时间。 所得结果绘制成表格如下： %tp/sts/s 校正Gs250.360.893 校正Gz250.380.9 连续闭环系统与采样系统对单位阶跃输入响应性能几乎一样，差距不大，但得选 合适的离散化方法才行，上图是用双线性变换法。 7.讨论采样周期的不同选择对系统控制性

20、能的影响讨论采样周期的不同选择对系统控制性能的影响: Ds=tf(10*1 6.512,1 11.499); Gs=tf(10,1 5 0); %离散系统的阶跃响应 Ts=0.02;i=100;%把 Ts=0.02,0.04,0.06,0.08 分别代入 Dz=c2d(Ds,Ts,tustin);%双线性变换 Gz=c2d(Gs,Ts,zoh);%零阶保持器 Close_Z=minreal(feedback(Dz*Gz,1); Y=dstep(Close_Z.num1,Close_Z.den1,i); plot(Ts*(1:i),Y,-.r); hold off Ts=0.02s Ts=0.0

21、4s Ts=0.06s Ts=0.08s 将上述结果绘制成表格如下： 采样周期/s超调量%峰值时间tp/s调节时间ts/s Ts=0.02250.380.9 Ts=0.0430.40.40.92 Ts=0.0636.70.420.96 Ts=0.0842.70.450.98 根据表格可知，采样时间越短，性能越好。 8.如控制器改为如控制器改为 PID 控制器，请确定满足性能指标的控制器，请确定满足性能指标的 PID 控制器参数。控制器参数。 8.1 用最优 PID 控制法设计： 2 312 2 312 2 32 312 *3223 c(s) p(s)= cs = 101010( ) 0( )

22、( ) 1( ) 0( )(5 10)(10 10)10 ( )1.752.15 nnn PIDG K SK SK S K SK SKGc s Gs s Gc s GsSK SK SK ITAE D sSSS 设最优控制器及前置滤波器G p(s) 令G1，由于 G（） 闭环传递函数为 当采用指标时，查表得，最优控制多项式为： 为了满足调节时间的设 *32 2 123 2 2 32 4.4 1, =10 ( )17.52151000 ( ) 20.5,100,1.25 12.52051000 ( ) 17.52151000 ( ), ns n n t D sSSS s KKK SS s SSS

23、Gp s ITA 计要求。需要选取的合适值。由于 可选，则最优特征多项式为 将上式由于的分母多项式相比，应有 所以 为了进一步改善系统性能，在系统输入端接入前置滤波器从而 具有预期的最优 32 2 ( )( ) 0( )1000 ( ) 1( ) 0( )17.52151000 80 ( ) 16.480 E Gp s Gc s Gs s Gc s GsSSS Gp s SS 指标。 由此得 性能验证： K1=20.5;K2=100;K3=1.25; G0=zpk(,0 -5,10); Gc=tf(K3 K1 K2,1 0); Gp=tf(80,1 16.4 80); sys0=feedbac

24、k(Gc*G0,1);%无前置滤波器时的闭环系统 sysn=feedback(G0,Gc);%扰动端传递函数 sys=series(sys0,Gp);%加入前置滤波器之后的闭环系统 figure(1); subplot(2,1,1);step(sys0);grid;%绘制无前置滤波器的单位 阶跃输入响应曲线 subplot(2,1,2);step(sysn);%绘制无前置滤波器的单位阶跃 扰动响应曲线 axis(0,1.2,0,1.5);grid; figure(2); subplot(2,1,1);step(sys);grid;%绘制带前置滤波器的单位阶 跃输入响应曲线 subplot(2,

25、1,2);step(sysn);%绘制带前置滤波器的单位阶跃 扰动响应曲线 axis(0,0.8,0,1.5);grid; 无前置滤波单位阶跃响应：调节时间为 0.656s，满足设计指标。 有前置滤波单位阶跃响应：调节时间为 0.564s，满足设计指标。 9.如希望尽可能的缩短系统的调节时间，请设计相对应的如希望尽可能的缩短系统的调节时间，请设计相对应的 最小拍控制器，并画出校正后系统的阶跃响应曲线最小拍控制器，并画出校正后系统的阶跃响应曲线： 9.1 采用无纹波最少拍系统设计 概述：无纹波条件 )( )( )()(1 )( )( )( )( zG z zGzD zD zR zU zu )(

26、)( )( zG zG zG D N )( )()( )( )( )( zG zGz zG zU z N D u 的全部零点。应包括还必须 统的三个条件，计满足有纹波最少拍系无纹波最少拍系统的设 )()(zGz 详细设计： 10 ( )0.02 (5) G sTss S S Gc=zpk(,0 -5,10); Gzoh=c2d(Gc,0.02,zoh) 0. (z+0.9672) - (z-1) (z-0.9048) 11 11 (1) 0.0019351 0.9674 G( ) (1)(1 0.9048) zz z zz 广义对象的脉冲传递函数 11 1 1 2 2( ) ( )(1 0.9

27、674)( ) 1( )(1-)( ) z zzzF z zzF z （）闭环脉冲传递函数 -1 12 1 1 2 1 1 2 11 11 1 ( )( ) ( )1 ( ) .5 ( )1 0.491 ( )0.508(1 0.9674) 1( )(1)(1 0.491) 262.53 1 0.9048 1( ) ( ) ( ) 1( )1 0.4 F zF z F F zbZ z F F zZ zZZ zZZ z z D z G z z 和应取项数最少而能保证上式均成立的Z 的多项式 (z)=a 由可求得 (z)=0 08 由上式可得 （3）通过上式计算可得 1 91z 对设计结果验算)4

28、( 1 11 1 123 (1( )( ( )1 0.419 1 ( )( ) ( )0.508(1 0.9674) 1 0.508 zR zZ Y zz R zZZ Z ZZZ 检验误差序列 E(z)= 系统响应 由此可见，系统从第二拍开始跟踪上输入信号。 具体程序： Gc=zpk(,0 -5,10); Gzoh=c2d(Gc,0.02,zoh); Gs=zpk(,0 -5,10); Ts=0.02; Gz=c2d(Gs,Ts,zoh); Z=0.9048; P=-0.491; K=262.53; Dz=zpk(Z,P,K,Ts,0.02); Close_Z=feedback(Dz*Gz,1); step(Close_Z); axis(0,0.12,0,10);grid; 有图可知，用两拍就可实现误差为零，调节时