t=[t,t(i)+h]; y=[y,c*x]; end
plot(t,y); grid
2，二阶龙格-库塔法仿真（还有一个例题见试卷，注意次仿真方法的表达式不唯一，注意看清题目）
y(k+1)=y(k)+h/2*(k1+k2),k1=f(t(k),y(k)),k2=f(t(k)+h,y(k)+hk1) 例题：a=[-14,9,10;12,-9,10;24,-24,-18]; b=[2.8;4;12]; c=[1,1,1]; x=[0;0;0];Tf=10; h=0.01; u=100;t=0;m=Tf/h; y=0; for i=1:m k1=a*x+b*u; k2=a*(x+h*k1)+b*u; x=x+h*(k1+k2)/2; y=[y,c*x]; t=[t,t(i)+h]; end
plot(t,y); Grid
3，四阶龙格-库塔法仿真
a=[-14,9,10;12,-9,10;24,-24,-18]; b=[2.8;4;12];c=[1,1,1];d=0; x=[0;0;0];tf=10;h=0.01;t=0; u=100;m=tf/h;y=0; for i=1:m k1=a*x+b*u;
k2=a*(x+h*k1/2)+b*u; k3=a*(x+h*k2/2)+b*u; k4=a*(x+h*k3)+b*u;
x=x+h*(k1+2*k2+2*k3+k4)/6; y=[y,c*x]; t=[t,t(i)+h]; end
plot(t,y); grid
6
六，可控标准型的状态空间参数
3s2?2已知传递函数为：G(s)?，编写程序求其可控标准型的状态空
4323s?5s?2s?1间参数A,B,C,D。 >> a=[3 0 2];
>> b=[3 5 2 1]; （分子和分母的输入） >> b=b/a(1);a=a/a(1); >> n=length(a)-1; >> A=a(2:n+1);
>> A=[rot90(rot90(eye(n-1,n)));-fliplr(A)] >> B=[zeros(1,n-1),1]' >> m1=length(b);
>> C=[fliplr(b),zeros(1,n-m1)]
七，面向数字结构图的数字仿真
1..典型闭环系统仿真实现
1.输入数据
a=[a0,a1,…,an]; b=[b0,b1,…,bm]; X=[0,0,…,0]’;（n-1个） V=v0; T0=0;Tf=10; h=0.01; R=r; y=0;t=T0;
2.闭环系数矩阵
b=b/a(1);a=a/a(1); n=length(a)-1; A=a(2:n+1);
A=[rot90(rot90(eye(n-1,n))); -fliplr(A)];
B=[zeros(1,n-1),1]’; m1=length(b);
7
C=[fliplr(b),zeros(1,n-m1)]; Ab=A-B*C*V;
3.运用四阶龙格－库塔公式求解 N=round((Tf-T0)/h); for i=1:N
k1=Ab*X+B*R;
k2=Ab*(X+h*k1/2)+B*R; k3=Ab*(X+h*k2/2)+B*R; k4=Ab*(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
plot(t,y);
grid
2.，复杂联接的闭环系统数字仿真
例 已知系统结构图如下图所示，写出系统联结矩阵和输入联结矩阵。wo,wc,w.
8
注意：方框中每个都为传递函数，而不是数值 还有一个例题见试卷最后一题。
八，系统校正
1，某系统开环传函如右，求它的幅值裕量和相位裕量，并求闭环阶跃响应。 num=3.5;den=[1 2 3 2]; G=tf(num,den); bode(G);grid;
[Gm,Pm,Wcg,Wcp]=margin(G) G_close=feedback(G,1); step(G_close)
2，已知系统开环传递函数如下，绘出伯德图，并讨论其稳定性 >> k=1.5;num=1; den=poly([0,-1,-2]); w=logspace(-1,1,100); [m,p]=bode(k*num,den,w); subplot(2,1,1);
semilogx(w,20*log10(m)); grid;
ylabel('增益(dB)');
9
subplot(2,1,2); semilogx(w,p); grid;
xlabel('频率(rad/s)'); ylabel('相角(deg)');
3，超前校正
k0=1000；要求43 做原系统的Bode图与阶跃响应曲线，检查是否满足题目要求 k0=1000; num=1; den=conv(conv([1,0],[0.1,1]),[0.001,1]); figure(1); margin (k0*num,den);hold on figure(2);s1=tf(k0*num,den); sys=feedback(s1,1);step(sys);hold on 求超前校正器的传递函数 >> k0=1000;num=1; den=conv(conv([1,0],[0.1,1]),[0.001,1]); g0=tf(k0*num,den); [m,p,w]=bode(g0); [gm0,pm0,wcg0,wcp0]=margin(g0); r=45; gama=r-pm0+5; gama=gama*pi/180; a=(1+sin(gama))/(1-sin(gama)); wc=spline(m,w,1/sqrt(a)); T=1/(wc*sqrt(a)); gc=tf([a*T,1],[T,1]) 校验系统校正后是否满足题目要求。 >> g=g0*gc; >> figure(1); >> margin(g); 计算系统校正后阶跃响应曲线和性能指标。 >> sys=feedback(g,1); >> figure(3); >> step(sys); 4，滞后校正 k0=30；r>40;wc>2.3 10