实验二 面向系统结构图的连续系统数字仿真实验 matlab程序

4.离散相似法编程
时域离散编程求解 程序：
clear num=[8 10]; den=[0.1 1 0 0]; [A,B,C,D]=tf2ss(num,den) sysc=ss(A,B,C,D); T=0.05; sysd=c2d(sysc,T); Ad=sysd.a; Bd=sysd.b;
0 0.5 1
程序：
clear all A = [-1000.25 999.75 0.5;999.75 -1000.25 0.5;0 0 0]; t = 0; x = [1 -1 1]; h = 0.05; M = round(10/h); for k = 1:M tt = t(k) + h; x1 = -exp(-0.5*tt)+exp(-2000*tt)+1; x2 = -exp(-0.5*tt)-exp(-2000*tt)+1; x3 = 0; temp = [x1 x2 x3]; x = [x; temp]; t = [t; tt]; end T = 0.003125; m = 4; t_frog = 0; x_frog = [1 -1 1]; start_time_frog = clock; F = expm(A*T); for k = 0:m tt = 2^(k) * T; temp = F*x_frog(1,:)'; temp = temp'; x_frog = [x_frog; temp]; t_frog = [t_frog; tt]; if k<m F=F*F; end end qT = 2^(m)*T; for k = 1:M-1 tt = t_frog(m+k+1) + qT; temp = F*x_frog(m+k+1,:)'; temp = temp'; x_frog = [x_frog; temp]; t_frog = [t_frog; tt]; end pass_time_frog = etime(clock,start_time_frog); hh = 0.001; end pass_time_rk4=etime(clock,start_time_rk4); t_display = t(2:31,:); x1_display = x(2:31,1); x1_frog_display = x_frog(6:35,1); x1_rk4_display = x_rk4(51:50:1501,1); disp(' 时间 解析解 蛙跳法 RK4法 ') t_rk4 = 0; x_rk4 = [1 -1 1]; start_time_rk4=clock; for k = 1:50*M tt = k*hh; K1 = A*x_rk4(k,:)'; K2 = A*(x_rk4(k,:)'+hh*K1/2); K3 = A*(x_rk4(k,:)'+hh*K2/2); K4 = A*(x_rk4(k,:)'+hh*K3); temp = x_rk4(k,:)'+hh*(K1/2+2*K2+2*K3+K4)/6; temp = temp'; x_rk4 = [x_rk4; temp]; t_rk4 = [t_rk4; tt];
连续系统离散相似法的数字仿真实验
2.RK4 法不同 C 系统阶跃响应
1.4
1.4
1.2
1.2
1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
0
1
2
3
4
5
6
7
8
9
10
0
0
1
2
3
4
5
6
7
8
9
10
C=0 时系统阶跃响应曲线
1.5
C=0.1 时系统阶跃响应曲线
1.6 1.4 1.2
1
1 0.8 0.6
0
5
10
15
20
25
30
35
40
0
5
10
15
20
25
30
35
40
放大倍数
1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0
5
1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0
0
5
10
15
20
25
30
35
40
45
50
0
5
10
15
20
25
30
35
40
45
50
实验三 1.simulink 模型
T=0.5;
% 输入各环节传递函数参数并利用tf函数表示它们 num1=[1]; den1=[1 4 8]; sys1=tf(num1,den1);
0.5
y(t)
num2=[1]; den2=[1 4 8]; sys2=tf(num2,den2); num3=[1]; den3=[1 4]; sys3=tf(num3,den3);
运行结果： C=0
1.5
Cd=sysd.c; Dd=sysd.d; X=[0;0;0]; yt=0;tt=0; R=1; M=10/T; for k=1:M E=R-yt(end); X=Ad*X+Bd*E; Y=Cd*X+Dd*E; tt=[tt,k*T]; c=0.1 if (-c<Y<c) Y=0; else if (Y<-c) Y=Y+c; else (Y>c) Y=Y-c; End end yt=[yt,Y]; end ¿
0
1
2
3
4
5
6
7
8
9
10
C=0.1
1.5
1
0.5
0
0
1
2
3
4
5
6
7
8
9
10
C=0.5
1.5
1
0.5
0
0
1
2
3
4
5
6
7
8
9
10
C=1.0
plot(tt,yt,'k'); grid on;
1.4 1.2
1
0.8
0.6
0.4
0.2
0
0
1
2
3
4
5
6
7
8
9
10
实验四
连续系统的快速数字仿真实验
1.增广矩阵法与 RK4 法仿真比较
14
3
1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0
1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0
0
5
10
15
20
25
30
0
5
10
15
20
25
30
放大倍数
1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0
4
1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0
0.4 0.3
% 使用双线性变换法将各环节的传递函数转换为脉冲传递函数 sysd1=c2d(sys1,T,'tustin'); sysd2=c2d(sys2,T,'tustin'); sysd3=c2d(sys3,T,'tustin'); % 求出系统的闭环离散模型 K=256; sysd=K*sysd1*sysd2*sysd3; % 求出系统的开环脉冲传递函数 sysd_close=feedback(sysd,1); % 求出系统的闭环脉冲传递 函数，注意默认的是负反馈 A=sysd_close.den{1}; % 闭环脉冲传递函数分母多项式的系数 行向量 L=length(A); % 闭环脉冲传递函数分母多项式的系数行向量的维 数 B=sysd_close.num{1};% 闭环脉冲传递函数分子多项式的系数行 向量 B=B/A(1); % 对分子多项式系数进行归一化处理 A=A(2:L)/A(1); % 对分母多项式系数进行归一化处理 R=zeros(L,1); % 建立参与递推运算的输入信号序列存储列向量 R(1)=1; % 在0时刻的输入信号为1，0时刻之前的输入信号为0 Y=zeros(L-1,1); % 建立参与递推运算的输出信号序列存储列向 量 M=20/T; % 计算仿真递推次数 yt=0;tt=0; % 建立仿真计算时刻和仿真输出的初值 for k=1:M r=1; % r(t)=1(t)的离散序列，即在任意一个采样时刻，输 入值均为1 R=[r;R(1:L-1)]; % 刷新参与递推运算的输入信号序列 yk=-A*Y+B*R; % 递推计算 Y=[yk;Y(1:L-2)]; % 刷新参与递推运算的输出信号序列 yt=[yt,yk]; % 记录各仿真时刻的系统输出值 tt=[tt,k*T]; % 记录各仿真时刻 end sim('x4_2.mdl'); % 运行Simulink仿真模型，得到数值积分解 y plot(t,y,'k',tt,yt,'k*'); % 绘制仿真结果曲线 grid on; legend('RK4法标准解', '双线性变换法的解'); xlabel('time(s)'); ylabel('y(t)');
data = [t_display, x1_display, x1_frog_display, x1_rk4_display]; disp(data);