自动控制原理控制新版系统数字仿真.doc

自控原理实验四：控制系统数字仿真
一、实验目
通过本实验掌握运用四阶龙格-库塔法进行控制系统数字仿真办法，并分析系统参数变化对系统性能影响。
二、实验办法
1、四阶龙格——库塔法

四、实验设备
硬件： PC机一台
软件：MATLAB软件，Microsoft Windows XP。

五、实验报告
1、由公式

计算出分别为5%、25%、50%时分别为0.690、0.404、0.216。
画出在K为0到时闭环根轨迹，如下图

再画出分别为0.690、0.404、0.216阻尼线，求出阻尼线与根轨迹3个交点。
则可求出K分别为30.88、59.25和103.55。
K也可以这样算：若系统有超调量，则由主导极点法可知原系统可简化为二阶系统，两个闭环极点共轭接近虚轴，另一种闭环极点远离虚轴，分别设为，则
，故
，即可算出K。
2、依照仿真成果，绘制阶跃响应曲线并求出(the settling time)和σ%(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 settling time = 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:80
the over shoot =23.117615 %.
the settling time =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:80
the over shoot is =45.722310 %
the settling time =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,c
A=[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:N
for 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;
end
settlingtime=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;

开始
输入：ai、ci、h、m、N
计算 A、b、c
，
i=0，j=0，t=0，y=0
输出：t，y
i=i+1,j=j+1,t=t+h
计算、、、

j=m
N
计算
输出：t、y
Y
j=0

结束
Y
N
程序流程图：