控制系统数字仿真-基于MATLAB实现四阶龙格库塔法

系统结构图如图



则先画根轨迹图：

den=[1 10 25 0]
num=[0 0 0 1]
rlocus(num,den)



用四阶龙格库塔法进行仿真，分别求超调量为5%，25%和50%时的K值

y=[0 0];
k=1;
while max(y)<=1.5 %%1.05(5%超调),1.25(25%超调),1.5(50%超调)
num1=[k];
den1=[1 10 25 0];
[num,den]=feedback(num1,den1,1,1);
[A,B,C,D]=tf2ss(num,den);
x0=[0;0;0];
v=1;
tf=15;
t0=0;
h=0.1;
r=1;
x=x0;
y=0;
t=t0;
for i=1:tf/h
K1=A*x+B*r;
K2=A*(x+h*K1/2)+B*r;
K3=A*(x+h*K2/2)+B*r;
K4=A*(x+h*K3)+B*r;
x=x+h*(K1+2*K2+2*K3+K4)/6;
y=[y;C*x];
t=[t;t(i)+h];
end
k=k+1;
end
%stepvalue=1 [OSValue, OSIndex] = max(y); OverShoot = (OSValue - stepvalue)/stepvalue*100;PeakTime = t(OSIndex);text(PeakTime, (1 + OverShoot/100 + 0.05)*stepvalue, sprintf('超调量%.2f%%',OverShoot))
plot(t,y)
k
mp=max(y); %峰值时间
tp=spline(y, t, mp)
cs=length(t);%稳态值
yss=y(cs)%超调量
ct=(mp - yss)/yss % 超调量和峰值时间
最后可以在命令行窗口看到输出的K值、稳态值、超调量和峰值时间，并画出阶跃响应的图形