龙格库塔法

BUAA

控制系统数字仿真

丁健39172223

2011/12/14

一、实验目的

通过本实验掌握利用四阶龙格-库塔法进行控制系统数字仿真的方法，并分析系统参数改变对系统性能的影响

二、实验方法

1、四阶龙格-库塔法

若一阶微分方程如下：y t=f(t,y(t))

y t=y0则在t n+1

(t n+1>t0)处，y(t n+1)的近似值为：y n+1=y n+

h

6

(k1+2k2+2k3+k4)

式中：h=t n+1?t n k1= f(t n,y n)

k2= f(t n+1

2h,y n+1

2

hk2)

k4= f(t n+1

2h,y n+1

2

hk1)

k3= f((t n+1

2h,y n+1

2

hk3)

n=0,1,2,3......

如果微分方程是如下形式的向量微分方程：X t=F(t,x t,u(t))

X t=X0其中，X t为m维向量，t,u(t)均为标

量，则在t n+1处（t n+1>t n）, X t的近似值为：

X n+1=X n+h

6

(K1+2K2+2K3+K4)

式中：h=t n+1?t n K1= F(t n,X n,u(t n))

K2= F(t n+1

2h,X n+1

2

hK1,u(t n))

K4= F(t n+1

2h,X n+1

2

hK2,u(t n))

K3= F((t n+1

2h,X+1

2

hK3,u(t n))

n=0,1,2,3,……

2、控制系统数字仿真

设系统的闭环传递函数为：φs=y(s)

u(s)=c1s n+1+c2s n+2+?+c n+1+c n s n+as n?1+?+a n?1+a n

引入中间变量V(s)则上式可化为：y(s)

u(s)=y(s)

v(s)

×v(s)

u(s)

令v(s)

u(s)=1

s+as+?+a n?1s+a n

, y(s)

v(s)

=c1s n+1+c2s n+2+?+c n+1s+c n

经过运算、整理可得如下向量形式：x t=AX t+bu(t)

y t=cX(t)

x0=0

所以，可得F(t,X(t),u(t))=AX(t)+bu(t)

三、实验内容

已知系统结构如图1所示：

若输入为单位阶跃函数，计算当超调量分别为非作5%，25%和50%时K的取值，（用主导极点方法估算），并根据确定的K 值在计算机上进行数字仿真。

四、实验步骤

图1

1、确定K值

（1）、σ%=5% 时：σ%=e?πζ/1?ζ2×100%=5% 得ζ=±0.69（可以利用Matlab的“solve”指令快速得到结果）

β=cos?1ζ≈45°

所以，设主导极点为：s1,2=a±ja 代入D(S)=S3+10S2+25S+K=0 得

D(a+ja)=(a+ja)3+10(a+ja)2+25(a+ja)+K=0 解得

K=31, a=-1.43

（2）、σ%=25% 时：同理可得：ζ=±0.403K=60,s1,2=-1.1±j2.5

（3）、σ%=50% 时：同理可得：ζ=±0.215 K=102, s1,2=-0.75±j3.4

2、根据仿真结果，绘制阶跃响应曲线

仿真过程通过Matlab实现,编写的程序代码分为两部分，如下：

（1）、主程序：

%----------------------%主程序---------------------

%系统的开环函数为：

% k

% -----------------------

% s(s+5)^2

num1=[31]; %根据之前计算出的k值对num1进行赋值

den1=conv([1,0],conv([1,5],[1,5]));

[num,den]=feedback(num1,den1,1,1);

[A,B,C]=tf2ss(num,den); %通过tf2ss指令将传递函数变为矩阵形式

x0=[0;0;0];

v=1;

tf=10;

t0=0;

h=0.001;

r=1;

[t,y]=Runge_Kutta(A,B,C,x0,h,r,v,t0,tf);

plot(t,y),grid

（2）、Runge_Kutta.m函数：

%------------------ Runge_Kutta.m----------

%该函数应该保存为M文件

%定义函数[t,y]=Runge_Kutta (A,B,C,x0,h,r,v,t0,tf);其中，（A，B，C）为系统的系数矩阵，x0为输入，h为仿真步长，r为输入信号幅值，t0为仿真的起始时间，tf为终止时间;t为仿真时间，y为系统输出function [t,y]=Runge_Kutta(A,B,C,x0,h,r,v,t0,tf);

%Runge_kutta.m

x=x0;

y=0;

t=t0;

fori=1:tf/h

K1=A*x+B*r;

K2=A*(x+(h/2)*K1)+B*r;

K3=A*(x+(h/2)*K2)+B*r;

K4=A*(x+h*K3)+B*r;

x=x+(h/6)*(K1+2*K2+2*K3+K4); y=[y;C*x]; t=[t;t(i)+h]; end

本次仿真实验的环境为Matlab2011b ，仿真过程如下：

a 、新建m 文件，输入Runge_Kutta.m 函数：并保存，注意文件名应与函数名相同，如图2

b 、在command Windows 中输入主程序代码，如图3

图 3

C 、得到系统的阶跃响应图像，如图

4

图 2

图 4

d、通过修改主程序中的k值，即对“num1”的值进行相应的修改，可以得到其他超调量下的系统阶跃响应曲线。阶跃曲线分别如下：

σ%=25%的阶跃响应图形 k=31 σ%=25%的阶跃响应图形 k=60

σ%=50%的阶跃响应图形 k=102

3、根据仿真结果，求ts 和σ%

取±5%误差带：在Figure 中进行分析取值，得到如下结果：

五、总结与分析

本次实验是通过Matlab 进行控制系统的数字仿真，通过对Matlab 的学习，我了解到Matlab 这款强大的软件在自动控制、信号处理、系统分析、数值运算、图形处理等方面有着无与伦比的优势，在系统的数字仿真方面，Matlab 也有这很多种方法，龙格-库塔法只是其中的一种。

在工程上面，我们遇到的很多控制系统都是高阶系统，用解微分方程的方法求得高阶系统的时间响应是十分困难的，但是通过龙格-库塔法可以对高阶系统进行有效的处理。

从实验结果中发现，通过龙格-库塔法进行的数字仿真实验也出现了一定的误差，但是通过分析发现，这种误差主要并不是由龙格-库塔法造成的，造成这一误差的主要来源是对主导极点的求解方法存在较大的误差，因此对于K 值的求解存在较大的误差。

实验后，为了解决由于主导极点的求解误差而产生的误差，发现了新的计算K 值的方法，具体步骤如下：

（1）、利用Matlab 中的“rlocus ”指令来绘制系统的根轨迹，代码如下： num=[1];

den=[1,10,25,0]; rlocus(num,den); 得到如图5所示的根轨迹

（2）、在根轨迹图中右键-grid,得到标有阻尼比数值的直线，如图6所示

（3）、利用zoom out 、zoom in 和data cursor 找到响应的超调量“overshoot ”，即可得到相应的K 值，如图7所示

图 7

图5

图6

（4）、改进后方法得到的结果如下表所示

这样的方法可以较为准确的得到K值，便于精确的计算，当然，这种方法也存在一定的误差，主要是由于作图过程中的步长选择造成的（步长若过大，在把图形放大后原本光滑的曲线会变成折线），在具体操作过程中，我们可以通过设置不同的步长来达到不同的精度要求

除此之外，通过实验我发现，高阶系统可以在特定情况下近似为一、二阶系统，并且在工程上保证了一定的精度，这样就可以利用一、二阶系统的分析方法来分析高阶系统，同样可以达到很高的精度，在次不再赘述。

总之，本次实验收获很大，我对于Matlab的掌握程度有了前所未有的提升，以前做实验的时候，总是利用别人已经编写并调试好的程序，从来没有思考过为什么要这样写。这次实验的过程中，我首先对Matlab 进行了深入的学习和练习，并且自己一句一句的编写程序，对于以前从来没有编写过Matlab程序的我来说，这个过程充满了艰辛，一遍遍的失败几乎让人崩溃，通过不断的查阅书籍和资料，当电脑屏幕上终于显示出正确的阶跃响应曲线的时候，我的心情异常激动。Matlab是一款很好很强大的软件，在我后续的学习过程中发挥了重要的作用，掌握好这个软件，对于之后的学习有很大的帮助。