自控原理实验四：控制系统数字仿真一、实验目的通过本实验掌握利用四阶龙格-库塔法进行控制系统数字仿真的方法，并分析系统参数改变对系统性能的影响。
二、实验方法1、四阶龙格一一库塔法卄一阶微分方程如2 则在切L(如A巾｝处，F仇+J的近似值为:加1 =儿-：馆-2雇+比-駐)〔m 0式中：/r = f HUK=fi^y n)孤三/(匚十可k儿十—^1)■■俎二/久卄片乩儿亠:展丿■*h=f(f. +九儿十碣｝如果微分方程是如下丿枚式的向馆微分方程, Jgf) = d⑴少⑴)U(O)=兀M中＜X(t)为E维向量，u⑴均为标m .则在匸处(gfgj的近似(f［为：兀+】二兀+ £ ［£ +皿4 2心+瓦］( 7-4)O(7-1)(7-3)式中：“也K严F（r”Kr』）K严弘+杯兀+*,心））AB 亠K3=F（r”+£・X”+£KyM（Fj）瓦=尸亿+力丄”+也3・?心））n = 01 .........2.控制系统数字仿真设系统的闭环传递函数为^如=凹=**宀…％ZM（$）S n+“s"T 十…+ 心_] +a”引入中间变量7(s)则上式叮化为：如二凹M(S) v(s) 令：型= ___________ ! ________H(5) s"十as"T+・..a”]S + a”誥=5严+巾严+…c”4q由以上两式吋得如下两个微分方程v w(r) + av^(『) + ••• + d”_p(r) + a n v(t) = w(r)W) = epi (r) + C2V(W_2) (/) + •■• + c』(r) + e…v(r)令:v(B_1)(0) = v(n_2)(0) = - = v(0) = v(0) = 0H (0 =啲，心(0 = "(『)• •，耳(0 = E (0则(7-8)式可化为如下一阶微分方程组：AW = ^2(r)右”) = x3(r)亢-1 ⑴=X” (0九(0 = 一讣(r) -务丙⑴------ 叭(r) + u(t) (7~9)式口J丐成：J(0 = c”“(r) +存辺⑴+…中左)(7-5)v(s)M(S)(7-6)(7-7)(7-8)(7-9) (7-10)(7-11)方程(了-10)和(7-11)可写成如下向量形武：x(r) = -IV(r)+ iu(r)y(t) = cX(t)”(2-12)x(0)二0这里疋⑴为H錐列向量F M(f)为标量,/为?常数矩阵. b为"维列向量，£为鬥维列向量，并分别具育如下形式・x(t) =巧(0伙)* b =<1耳(f)I■0 1 0… o00 1 0A =*曲i占bQ0 0 (1)5 -7一1 一口n-2 0e =l57.…订对比(7~3)式口J得: F(t, JT(t), u(0) = + 6w(O三冒实验内容已知系统结构如图7-1R q i c辄即” ~r*图7_1若输入为单位阶跃帼数.计算肖超调蜀分别为非作叭，25%, 和50%时K的取值(用主导极点方法佔算)，并根据确定的K值在计算机上进行数字仿真。
四、实验设备硬件：PC机一台软件：MATLAB^件，Microsoft Windows XP 。
五、实验报告 1、由公式% e 1 2100%计算出 %分别为5%、25% 50%时的 分别为0.690、0.404、0.216K画出G s2在K 为0到，时的闭环根轨迹，如下图s s 5再画出 分别为0.690、0.404、0.216的阻尼线，求出阻尼线与根轨迹的3个交点。
则可求出K 分别为30.88、59.25和103.55。
K 也可以这样算：若系统有超调量，则由主导极点法可知原系统 可简化为二阶系统，两个闭环极点共轭靠近虚轴，另一个闭环极恥）=-_处——y 旦点远离虚轴，分别设为 九为厂也，贝V ' 70胪+笑$+巧 .2 2 3 2 2、2©3己柑匚-三一flcot Locus-5 0 5 1CRea A XB(S 2 n S n )(S S 3) S (S 3 2 n )S (2 n S 3n)S n S 3，故2 2S 32 n 10,2 nS3n 25, n S3 k，即可算出 K o2、根据仿真结果，绘制阶跃响应曲线并求出 & (the settlingtime) 和彷 %(the over Shoot) ①当K = 30.88时， 以矩阵形式输入a: [10,25,30.88] 以矩阵形式输入c: [0,0,30.88]请输入步长h:0.025 请输入打印步长 mh 之m:8 请输入迭代次数N*m 之N:60 the over shoot % =4.425886 % the Settli ng time tS = 3.000000 S .②当K=59.25时，以矩阵形式输入a:[10,25,59.25]以矩阵形式输入c:[0,0,59.25]请输入步长h:0.025请输入打印步长mh之m:8请输入迭代次数N*m之N:80the over shoot %=23.117615 %. the settling time t s=3.150000 s.③当K=103.55时，以矩阵形式输入a:[10,25, 103.55]以矩阵形式输入c:[0,0,103.55]请输入步长h:0.025请输入打印步长mh之m:8请输入迭代次数N*m之N:80the over shoot is % =45.722310 %the settli ng time ts=4.950000 s六、实验结论可以看出，当系统其它参数不变，根轨迹增益K的值增加时，一对主导极点起作用；调节时间增大，响应速度变慢，快速性降低; 超调量增加，振荡加剧，稳定性变坏。
附录程序：function y=runge_kutta(f)%输入ai,ci,h,m,N.a=input(' 以矩阵形式输入a:\n');% 闭环传递函数分母的系数，除最高项系数 1 外。
cc=input(' 以矩阵形式输入c:\n');% 闭环传递函数分子的系数，元素数与 a 的相同。
h=input(' 请输入步长h:');m=input('请输入打印步长mh之m:');N=input('请输入迭代次数N*m之N:');%计算A,b,cA=[0 1 0;0 0 1;-a(3) -a(2) -a(1)];b=[0 0 1]';c=[cc(3) cc(2) cc(1)];u=1; % 时域形式的u(t) ，单位阶跃输入时u(t)=1.x=zeros(3,N*m);% 给X 初值t0=0; %t 初值t=t0+[0:N*m]*h;y0=0; %y 初值y=zeros(N*m,1);y(1,:)=y0;fprintf('t 的值为%f\n',t0); % 输出t,y 的初值。
fprintf('y 的值为%f\n',y0);f=fun;% 函数句柄for i=1:Nfor j=1:m k1=h*feval(f,t(j+m*(i-1)),x(:,j+m*(i-1)),u,A,b);k2=h*feval(f,t(j+m*(i-1))+h/2,x(:,j+m*(i-1))+k1/2,u,A,b); k3=h*feval(f,t(j+m*(i-1))+h/2,x(:,j+m*(i-1))+k2/2,u,A,b); k4=h*feval(f,t(j+m*(i-1))+h,x(:,j+m*(i-1))+k3,u,A,b); x(:,j+1+m*(i-1))=x(:,j+m*(i-1))+(k1+2*k2+2*k3+k4)/6;y(j+1+m*(i-1),:)=c*x(:,j+1+m*(i-1));end % 四阶龙格库塔算法fprintf('t 的值为%f\n',t(m*i+1));fprintf('y 的值为%f\n',y(m*i+1)); %打印(m*i+1)*h 步长时的t,y end[Y,k]=max(y); % 将一系列y 中的最大值赋给Y percentovershoot=100*(Y-1)/1; % 计算超调量 (对稳定系统的单位阶跃响应，终值为1)fprintf('the over shoot is %f %%.\n',percentovershoot); 打印超调量值i=length(t);while (y(i,:)>0.98)&(y(i,:)<1.02)i=i-1;endsettlingtime=t(i); % 找出2%误差带时的调节时间fprintf('the settling time is %f s.\n',settlingtime); % 打印调节时间值plot(t,y); % 利用数值解绘制单位阶跃响应曲线grid;title('step response');xlabel('t/s');ylabel('y/v');function f=fun(t,x,u,A,b)%给出函数关系式f=A*x+b*u;程序流程图:。