控制系统数字仿真实验指导书

1、淮 阴 工 学 院控制系统数字仿真实验指导书编者：王建国适用学院：电子与电气工程学院电子与电气工程学院2015年 6 月 1 日实验一面向一阶微分方程组的连续系统数字仿真实验一、实验目的1掌握以一阶微分方程组形式描述的连续系统的数字仿真方法和步骤。2. 了解如何编程实现对一阶微分方程组模型的仿真研究。二、实验内容单位反馈控制系统的结构如图A.1所示。图A.1 系统结构1根据系统闭环传递函数，求出对应的状态空间模型，并利用MATLAB中的ODE解函数编程，对系统进行仿真研究，并绘制对应的单位阶跃响应曲线。2考虑微分方程模型试编程对其进行仿真（0t100），并用不同的色彩、线型绘制仿真结果曲线。三

2、、预习要求1. 复习ODE解函数使用方法和调用格式。2复习一阶微分方程组仿真程序的编制方法。四、实验报告要求1. 整理各种实验条件下的数据和曲线。2. 写出实验报告。五、教案1利用MATLAB编程求出图A.1所示的系统状态空间模型的程序（文件名为Appl1_1_1.m）如下：% 这是实验一中第一个系统求取状态空间模型的程序clear% 输入三个环节的系统传递函数sys1=tf(1,0.5,1);sys2=tf(1,0.2,1);sys3=tf(1,1,0);% 求系统的闭环传递函数sys=sys1*sys2*sys3;sys_cl_tf=feedback(sys, -1);% 求系统的状态空间

3、模型sys_cl_ss=ss(sys_cl_tf)运行Appl1_1_1.m后，得到系统的状态空间模型为：据此，编制仿真程序（文件名为Appl1_1_1_fang.m）如下：% 这是实验一中第一个系统的仿真程序clearx0=0,0,0;tspan=0,10;t,x=ode45('Appl_1_1_func', tspan, x0);plot(t, 1.25*x(:,3);xlabel('t');ylabel('y(t)');grid;title('实验一 第一个系统的单位阶跃响应');其中，Appl_1_1_func.m程序如下

4、：% 这是实验一中第一个系统的对应的数学模型function xdot=Appl_1_1_func(t, x)xdot1(1)=-7*x(1)-2.5*x(2)-1.25*x(3)+1;xdot1(2)=4*x(1);xdot1(3)=2*x(2);xdot=xdot1'运行Appl1_1_1_fang.m后，得到系统的单位阶跃响应曲线如图A.1_1所示。图A.1_1 实验一第一个系统的单位阶跃响应2. 利用MATLAB编程编制求解系统的仿真程序（文件名为Appl1_1_2_fang.m）如下：% 这是实验一中第二个系统的仿真程序clearx0=10,10;tspan=0,100;t,

5、x=ode45('Appl_1_2_func', tspan, x0);plot(t, x(:,1), 'k', 'LineWidth',2);holdplot(t, x(:,2), 'r-.', 'LineWidth',2);legend('x(1)', 'x(2)');xlabel('t');ylabel('状态向量');grid;title('实验一 第二个系统的仿真输出');其中，Appl_1_2_func.m程序如下：% 这是

6、实验一中第二个系统的对应的数学模型function xdot=Appl_1_2_func(t, x)xdot1(1)=- x(1) *(x(1)*x(1)+x(2)*x(2)+x(2);xdot1(2)=-x(1)-x(2) * (x(1)*x(1)+x(2)*x(2);xdot=xdot1'运行Appl1_1_2_fang.m后，得到系统的仿真曲线如图A.1_2所示。图A.1_2 实验一第二个系统的仿真曲线实验2 面向系统结构图的连续系统数字仿真实验 一、实验目的（1）掌握以系统结构图形式描述的连续系统的数字仿真方法和步骤。（2）初步了解如何用仿真方法来分析系统的动态性能。（3）了解

7、不同的数值积分算法与仿真计算精度之间的关系。（4）学会一种初步寻求合理仿真步长的方法。（5）了解RK4法计算稳定性和步长的关系。二、实验内容单位反馈控制系统的结构如附录图1所示，其开环传递函数和闭环传递函数分别为在零初始条件下闭环系统单位阶跃响应的标准解为1 按实验目的、要求和已知条件，建立系统的Simulink模型。2. 按下面的经验公式（1）初选仿真步长 （为此，需要学习MATLAB中求取开环剪切频率 的函数margin的使用方法）。 （1） 为系统中反应最快的小闭环的开环剪切频率。 用margin函数编写如下程序求开环剪切频率：clearnum=1;den=0.1,0.7,1,0;G0=

8、tf(num,den);margin(G0)结果如图：由此我们得到。根据经验公式可得，仿真步长的取值范围为：。3. 选择RK4法，运行仿真模型，适当调整步长和仿真起止时间，以得到比较完整的过渡过程，观察纪录过渡过程的数据。根据实验内容2结果得到的仿真步长的取值范围：，结合实验内容1中建立的Simulink仿真模型，设定仿真参数时，我们选取定步长的ode4（即RK4）算法，仿真时间0.0到10.0，对比在仿真步长在之间变化时（选取h=0.06,0.10,0.14,0.18,0.22五组数据），步长对精度的影响。误差评价指标公式的程序：t(1)=0, yAE(1)=0;h=0.06(0.10,0.

9、14,0.18,0.22);M=10/h;For k=1:M t(k+1)= t(k)+h;yAE(k+1)=1-*exp(*t(k+1)- *exp(*t(k+1) .*sin( ;endyAE= yAEMAX_ERR=max(abs(yAE-y2);MEAN_ERR=sum(abs(yAE -y2)/size(yAE,1);ERR_NUM=MAX_ERR MEAN_ERR;求解得到的RK4算法步长与误差关系如下：仿真步长3.7394×10-63.4533×10-51.5475×10-45.2613×10-40.00142.4459×10-7

10、2.1378×10-69.3895×10-62.9336×10-57.4717×10-5上表中数据是我们按照RK4法步长选择经验公式，选取的仿真步长实验结果，为了比较在更小步长下RK4法精度变化，我们补做在=0.01，0.03,0.05时仿真，并比较精度，结果如下：仿真步长0.010.030.05680×10-52.0467×10-71.7180×10-6130×10-51.3853×10-81.1378×10-7通过以上两表比较，得出数值积分算法在舍入误差和截断误差二者的影响下，仿真总误差与步

11、长的关系并不是单调函数关系，而是一个具有极值的函数关系。同时可知，根据RK4法步长选取的经验公式选择步长，能让我们较为准确的找到最佳步长。4改变I环节的放大倍数（将1变为5），不断加大步长观测记录计算稳定性的变化。(1)h=0.01(2)h=0.1(3)h=0.25(4)h=0.4(5)h=0.5 经过以上五张图的对比得知，当仿真步长过大时，数值积分算法的解会发生错误。因此针对数值积分算法稳定性的要求，必须合理选择步长。5. 在相同的条件下，选择欧拉法和RK2法，再让仿真模型运行，观察纪录过渡过程的数据。我们选取=0.03,0.13,0.23三个步长，分别利用欧拉法、RK2法、RK4法进行仿真

12、，比较在相同步长情况下，各种方法的精度评价指标MAX_ERR、MEAN_ERR，实验结果如下：评价指标算法h=0.03h=0.13h=0.23欧拉法0.01540.07400.1473RK2法2.0535×10-40.00400.0138RK4法2.0467×10-71.0773×10-40.0017欧拉法0.00200.00970.0197RK2法2.7789×10-55.4921×10-40.0017RK4法6.9507×10-93.4069×10-64.6734×10-5由上表可以看出，在相同步长下仿真，RK

13、4法精度明显高于RK2法和欧拉法。三实验总结 本实验主要研究的是数值积分中RK4法的步长选取、精度问题和稳定性。从实验结果中我们得出，步长的选取对数值积分仿真的精度有着很大的影响。我们一般采用经验公式或选取仿真步长，以保证小误差仿真。如果步长选择不合理，会造成误差的积累越来越大，使计算出现不稳定情况，得出错误的结果。实验3 连续系统离散相似法的数字仿真实验一、实验目的（1）掌握以系统结构图形式描述的连续系时域离散相似法的数字仿真方法和步骤。（2）学会利用时域离散相似法分析线性和非线性控制系统的动态性能以及典型非线性环节对控制系统动态性能的影响。二、实验内容含死区环节的非线性控制系统的结构图如课

14、本所示。（1）具体见课本；（2）求出图A.2中传递函数对应的状态空间模型，并在该模型前加入虚拟的采样开关和零阶保持器，得到离散化状态空间模型。程序如下：clearsys1=tf(0.8 1,1 0)；sys2=tf(10,0.1 1)；sys3=tf(1,1 0)；sys= sys1* sys2*sys3;sysc=ss (sys);T=0.05sysd=c2d(sysc,T);（3）在c=0，0.1，0.5，1.0情况下，利用时域离散相似法编程完成对该系统的仿真研究。Ad=sysd.a;Bd=sysd.b;Cd=sysd.c;Dd=sysd.d;X=0;0;0;yt=0;tt=0;R=1;C

15、=0(0.1，0.5，1.0)M=30/T;for k=1:M yt1=0; if yt(end)>C; yt1=yt(end)-C; elseif yt(end)<-C; yt1=yt(end)+C; end U=R-yt1; X=Ad*X+Bd*U; Y=Cd*X+Dd*U; yt=yt,Y; tt=tt,k*T;endplot(tt,yt,'k'),grid三、预习要求 具体见课本。四、实验报告要求 具体见课本。实验四连续系统的快速数字仿真实验一、实验目的1熟悉快速数字仿真中双线性变换法、增广矩阵法和时域矩阵法的原理及步骤。2. 学会编程实现双线性变换法、增广

16、矩阵法和时域矩阵法，对控制系统进行仿真研究。二、实验内容1已知系统状态方程及初值条件 系统响应的解析解为试分别采用增广矩阵法和RK4法对系统进行仿真，并比较两种方法所得数值解的计算精度、所耗费的计算时间以及计算稳定性。（1）编写程序如下：clearA=-1000.25 999.75 0.5; 999.75 -1000.25 0.5; 0 0 0; t(1)=0; x1(1)=1; % 置解析解x1的初值x2(1)=-1;hh=0.05; M=round(10/hh); % 求解析解for n=1:M tt=t(n)+hh; t(n+1)=tt; x1(n+1)=-exp(-0.5*tt)+ex

17、p(-2000*tt)+1; x2(n+1)=-exp(-0.5*tt)-exp(-2000*tt)+1;end% 采用增广矩阵法求解T=0.003125; m=4; t_frog=zeros(1,M+5); x=zeros(3, M+5); t_frog(1)=0; x(:,1)=1; -1 ; 1; % 置状态向量的初始值star_time_frog=clock; F=expm(A*T); for k=0:m % 完成蛙跳式计算t_frog(k+2)= 2(k)*T;x(:, k+2)=F*x(:,1);if k<mF=F*F; end end qT=2(m)*T; for k=1:

18、M-1 t_frog(k+6)= t_frog(k+5)+ qT ;x(:, k+6)=F* x(:, k+5);endpass_time_frog=etime(clock, star_time_frog) % 采用RK4法求解h=0.001; % 置RK4法步长h=0.001st_rk4=zeros(1, 50*M+1); y=zeros(3, 50*M+1); t_rk4(1)=0; y(:,1)=1; -1 ; 1; % 置状态向量的初始值star_time_rk4=clock; for k=1:50*Mt_rk4(k+1)=k*h; % 自变量时间K1=A*y(:,k);K2=A*(y

19、(:,k)+h* K1/2);K3=A*(y(:,k)+h* K2/2);K4=A*(y(:,k)+h* K3);y(:, k+1)=y(:, k)+h*( K1+2* K2+2* K3+K4)/6; endpass_time_rk4=etime(clock, star_time_rk4) xlswrite('file.xlsx',x')（2）计算精度比较：下面三张图分别从excel表格中截取了解析解、RK4法和增广矩阵法最后一页的值进行比较。三者均是不断增大。由上图，增广矩阵法的计算精度比较高。（3）耗费计算时间比较>> y4_1pass_time_fro

20、g = 0.0010pass_time_rk4 =0.0870由以上的程序运行结果，增广矩阵法速度明显快于RK4法。（4）计算稳定性增广矩阵法是一种绝对稳定的算法，允许选择大采样周期T。 RK4法是条件稳定算法，步长h的大小除了与算法本身有关外，还与方程本身的性质有关。对于实际系统，由于其特征值不一定为实数。一般而言，步长h必须满足一定的条件，对于RK4法，步长的选择应满足：2已知控制系统的闭环传递函数为试采用双线性变换法（）对系统进行仿真，并与RK4法（）所得结果（作为标准解）进行计算精度比较。编写程序如下：clear;T=0.5;n1=8;d1=1 4 8;sys1=tf(n1,d1);n

21、2=8;d2=1 4 8;sys2=tf(n2,d2);n3=4;d3=1 4;sys3=tf(n3,d3);sysd=sys1*sys2*sys3;G=c2d(sysd,T,'tustin');A=G.den1;L=length(A);B=G.num1;B=B/A(1);A=A(2:L)/A(1);R=zeros(L,1);Y=zeros(L-1,1);M=10/T;yt=0; tt=0;for k=1:M R(k)=1; R=R(k);R(1:L-1); yk=-A*Y+B*R; Y=yk;Y(1:L-2); yt=yt,yk; tt=tt,k*T; endplot(tt,

22、yt,'k*');hold on;sim('y4');plot(tout,y4,'k')ttRK4法增广矩阵法00.0239843453605750.1200000000000000.5000000000000000.2641402369901490.41599999999999910.6736723136533440.7855999999999971.5000000000000000.9593505296240661.03471999999999421.044923198775741.0895999999999932.5000000000000

23、001.031322876999861.04129279999999331.007222706515090.9969894399999943.5000000000000000.9978893524483020.98741478399999440.9977727275212970.9946405887999944.5000000000000000.9993374528614051.00062679039999451.000068663152341.0016094904319945.5000000000000001.000137623044091.00058249256959461.0000494

24、84790660.9998793423257546.5000000000000000.9999998674351800.99980972487474670.9999922312849910.9999441673912267.5000000000000000.9999966842249181.00001928305245680.9999997736062331.0000210749898048.5000000000000001.000000407983231.00000473927007091.000000206933870.9999973080145069.5000000000000001.0

25、00000025992390.999997793096632根据上表，将增广矩阵法的结果与RK4法进行比较，发现t越大，计算的精度越高实验五采样控制系统的数字仿真实验一、实验目的1掌握采样控制系统数字仿真的特点。2. 了解数字控制器对系统动态性能的影响。3学会编制双重循环法的仿真程序。二、实验内容某单位反馈控制系统中，被控对象的传递函数为在单位速度信号输入下，选择零阶保持器在时设计的“最少拍”控制器为1. 按实验目的、要求和已知条件，建立系统的Simulink模型，并且编制双重循环法的仿真程序。 先用如下程序将连续系统转换为离散系统c2d函数，将系统离散化。clear;h=0.01; A=0

26、0;1 -1;B=10;0;G,H=c2d(A,B,h);运行程序后得到故连续系统被控对象的等价离散化状态方程为： 而题中控制系统的差分方程为： 上式中 双重循环法仿真程序：clearG=1 0;0.01 0.99;n=2;m=3;H=0.1;0.0005;c=0 1; h=0.01;T=1; N=round(T/h); TF=15; M=round(TF/T); x=0;0;u=zeros(n,1);e=zeros(m,1);r=zeros(M,1);A=0.282 0.718;B=0.543 -0.471 0.0999;t=0;xt=x;for a=1:M for b=1:N x=G*x+

27、H*u(1); xt=xt,x; t=t,b*h+(a-1)*T; end y=c*x; r(a)=a; ea=r(a)-y; e=ea;e(1:m-1); ua=A*u+B*e; u=ua;u(1:n-1);endx1t=xt(1,:);x2t=xt(2,:);plot(t,x1t,':k',t,x2t,'k');legend('x1', 'x2');grid;xlabel('tame(s)');2.分别运行双重循环法的仿真程序和Simulink模型，观察纪录系统动态性能的变化。（1）双重循环法仿真程序（2）simulink模型 从双重循环法和Simulink模型仿真的结果可以看出二者基本相同。系统的输出开始会有振荡，随时间的增加，振幅减小，最后趋于稳定。由此我们可以看出，双重循环法这种按连续系统离散相似算法建立连续部分各个环节仿真模型的方法，不仅仿真结果准确，而且在得到系统输出值的同时，也得到了被控对象内部状态变量的相应，实用性强，值得采用。3若数字控制器改为“无纹波”控制器试对两种情况下的输出进行比较。设定仿真时间为15S ，分别得出结果：由此可以看出，仿真中，在仿真快结束时系统会出现不稳定