系统建模与仿真期末考试试卷

(2014年12月)

院系： 年级： 班级： 学号： 姓名： 题号 1 2 3 4 5 总成绩 评卷人 得分 说明：请保留题目，在每个题目解答部分的空白处依次作答，并写清楚每个小题的题号；作答要给出程序代码、仿真结果；为了节约纸张环保，请缩小贴图、合理排版、双面打印。

1. （30分）已知系统的传递函数模型为：

（2s2)G(s)

（s3)(s4)（1）利用zp2ss()函数将该传递函数模型转化为状态空间模型；（5分） （2） 假设系统的输入为e-t：

①利用状态空间模型，假设状态的初始条件为[1 ;2]，t=0:0.1:4，求在 e-t输入下的状态响应、输出响应（利用subplot()函数将仿真曲线作在同一个窗口中）。（5分）

②利用laplace()函数求e-t的拉普拉斯变换；（5分）

③利用拉普拉斯反变换函数ilaplace()求系统输出的解析解，并根据此解析解仿真t=0:0.1:6系统输出响应；（5分）

④利用lsim()函数仿真t=0:0.1:6系统输出响应。（5分）

h (t ) ，利用y⑤假设系统的脉冲响应为 ( t )  h ( t ) * u ( t) 仿真t=0:0.1:6系统输出

响应。（5分） 解： (1):

[A,B,C,D]=zp2ss(-2,[-3,-4],2) A =

-7.0000 -3.41 3.41 0 B = 1 0 C =

2.0000 1.17 D = 0 (2):

①:

clc,clear;

[A,B,C,D]=zp2ss(-2,[-3,-4],2); G=ss(A,B,C,D); t=0:0.1:4; u=exp(-t);

[y,x]=lsim(G,u,t,[1:2]); subplot(2,1,1); plot(x,t); subplot(2,1,2); plot(y,t);

②

clc,clear; syms t; f=exp(-t); F=laplace(f); pretty(simple(F))

1 ----- s + 1 ③

clc,clear; syms s; f=exp(-s); F=laplace(f);

H=2*(s+2)/((s+3)*(s+4)); pretty(simple(ilaplace(F*H)))

2 4 -------- - --------

exp(3 t) exp(4 t) - ------------------- t + 1

clc,clear; t=0:0.1:6;

y=-((2./exp(3.*t)-4./exp(4.*t)))./(t+1); plot(t,y);

④

clc,clear;

[A,B,C,D]=zp2ss(-2,[-3,-4],2); G=ss(A,B,C,D); t=0:0.1:6; u=exp(-t); lsim(G,u,t);

⑤

2. （30分）假设系统的框图为：

其中k为系统的增益，用状态空间表示的系统G1、G2分别为：

122111 12xxx11u1x21u2G1:G2:4215

y12xuy24x11122（1）求系统的开环传递函数；（5分）

（2）绘制开环传递函数的根轨迹，并利用根轨迹确定闭环系统稳定的k的范围 ;

（5分）

（3）假设系统G1、G2的初始值均为0，分别取k=0.3、k=0.4，在t=0时刻加上

阶跃为1的输入，利用simulink仿真系统的输出响应(t=0-100)；（10分） (4) 假设系统G1、G2的初始值均为0，分别取k=0.3、k=0.4，在t=0时刻加上阶跃为1的输入，利用step( )仿真系统的输出响应(t=0-100)。（10分） 解： (1)

clc,clear;

a1=[-1 -2; 4 -2]; b1=[2; 1]; c1=[1 2];

d1=[1];

a2=[1 -1; 1 -5]; b2=[-1; 1]; c2=[-2 4]; d2=[0];

[num1, den1]=ss2tf(a1, b1, c1, d1); [num2, den2]=ss2tf(a2, b2, c2, d2); G1=tf(num1, den1); G2=tf(num2, den2); G=G1*G2

Transfer function:

6 s^3 + 46 s^2 + 208 s + 120 -------------------------------- s^4 + 7 s^3 + 18 s^2 + 28 s – 40 (2)

clc,clear;

a1=[-1 -2; 4 -2]; b1=[2; 1]; c1=[1 2]; d1=[1];

a2=[1 -1; 1 -5]; b2=[-1; 1]; c2=[-2 4]; d2=[0];

[num1, den1]=ss2tf(a1, b1, c1, d1); [num2, den2]=ss2tf(a2, b2, c2, d2); G1=tf(num1, den1); G2=tf(num2, den2); G=G1*G2; rlocus(G); sgrid

从图上可得k的范围为：0.33-inf; (3)

(4)

clc,clear; a1=[-1 -2; 4 -2]; b1=[2; 1];

c1=[1 2]; d1=[1];

a2=[1 -1; 1 -5]; b2=[-1; 1]; c2=[-2 4]; d2=[0];

[num1, den1]=ss2tf(a1, b1, c1, d1); [num2, den2]=ss2tf(a2, b2, c2, d2); G1=tf(num1, den1); G2=tf(num2, den2); G=G1*G2; T=tf(0.3,1);

GG=feedback(G,T,-1); step(GG,100)

clc,clear; a1=[-1 -2; 4 -2]; b1=[2; 1]; c1=[1 2]; d1=[1];

a2=[1 -1; 1 -5]; b2=[-1; 1]; c2=[-2 4]; d2=[0];

[num1, den1]=ss2tf(a1, b1, c1, d1); [num2, den2]=ss2tf(a2, b2, c2, d2); G1=tf(num1, den1); G2=tf(num2, den2); G=G1*G2; T=tf(0.4,1);

GG=feedback(G,T,-1); step(GG,100)

3. （15分）已知时间微分方程

tyy2t1 y(0)1（1） (5分)用Euler方法求解常微分方程初值问题，并将数值解和该问题的解析解（y(t)(t21)1/2）比较；

（2）（5分）利用四阶Runge-Kutta方法编程仿真y(t)； （3）（5分）利用ode45（）函数求解并仿真y(t)。 注：本题仿真时间取t0:0.02:5。

解： (1)

clc;clear; h=0.02; y(1)=1; t=0:h:5;

for n=1:length(t)-1 xn=t(n);yn=y(n);

y(n+1)=yn+h*(-yn*xn/(xn*xn+1)); end

t0=0:h:5;y0=(t0.*t0+1).^(-1/2); plot(t0,y0,'bo',t,y,'r*') legend('解析解','数值解')

(2)

clear;clc;

t0=0;tN=5;y0=1;h=0.02; t = t0: h : tN; N = length (t); for i = 1 : N-1 t1 = t0 + h; K1 = Runge(t0, y0);

K2 = Runge(t0 + h/2, y0 + h*K1/2); K3 =Runge(t0 + h/2, y0 + h*K2/2); K4 = Runge(t0 + h, y0 + h*K3);

y= y0 + (h/6)*(K1 + 2*K2 + 2*K3 + K4); t0=t1; y0=y; yy1(i)=y; end

plot (t, [1,yy1]);

function dy=Runge(t,y) dy=-t*y/(t*t+1); end

(3)

clc;clear; y0=1;

[t,y]=ode45('Runge',[0,5],y0); plot(t,y)

4. （15分）已知一个离散时间系统的输入输出数据如下表给出：

y(n)0.2y(n1)0.5y(n2)0.1y(n3)3x(n)5x(n1)4x(n2)4x(n3)（1）求它的级联结构形式；（5分） （2）求它的并联结构形式；（5分）

（3）分别利用直接型、并联型结构求阶跃输出响应（n=0:20），并比较响应曲线。

（5分）

解： (1)

clc;clear; n=[0: 20]; b=[3 5 4 4]; a=[1 -0.2 0.5 0.1]; [sos, g]=tf2sos(b, a)

sos =

1.0000 1.3958 0 1.0000 0.1765 0 1.0000 0.2709 0.9553 1.0000 -0.3765 0.5665 g =

3 (2)

直接型转换为并联型需要编写子程序dir2par.m、cplxcomp.m：

function [C,B,A]=dir2par(b,a) % 直接型转成并联型子程序 M=length(b); N=length(a);

[r1,p1,C]=residuez(b,a);

p=cplxpair(p1,10000000*eps); I=cplxcomp(p1,p); r=r1(I);

K=floor(N/2); B=zeros(K,2); A=zeros(K,3); if K*2==N

for i=1:2:(N-2)

Brow=r(i:1:(i+1),:); Arow=p(i:1:(i+1),:);

[Brow,Arow]=residuez(Brow,Arow,[]); B(fix((i+1)/2),:)=real(Brow);

A(fix((i+1)/2),:)=real(Arow); end

[Brow,Arow]=residuez(r(N-1),p(N-1),[]); B(K,:)=[real(Brow) 0];

A(K,:)=[real(Arow) 0]; else

for i=1:2:(N-1)

Brow=r(i:1:(i+1),:); Arow=p(i:1:(i+1),:);

[Brow,Arow]=residuez(Brow,Arow,[]); B(fix((i+1)/2),:)=real(Brow); A(fix((i+1)/2),:)=real(Arow); end end

function I=cplxcomp(p1,p2) % I=cplxcomp(p1,p2)

% 比较两个包含同样标量元素但（可能）有不同下标的复数对

% 本程序必须用cplxpair函数之后使用，以便重新排序频率极点向量 % 及其相应的留数向量； % p2=cplxpair(p1); % I=[];

for j=1:1:length(p2)

for i=1:1:length(p1)

if(abs(p1(i)-p2(j))<0.0001) I=[I,i]; end end end I=I';

clc;clear; n=[0: 20]; b=[3 5 4 4];

a=[1 -0.2 0.5 0.1]; [C, B, A]=dir2par(b, a) C =

40 B =

-7.7155 3.3355

-29.2845 0 A =

1.0000 -0.3765 0.5665 1.0000 0.1765 0 (3)

function y = parfiltr(C,B,A,x)

% PARALLEL form realization of IIR filters % [y] = parfiltr(C,B,A,x); % y = output sequence

% C = polynomial (FIR) part when M >= N

% B = K by 2 matrix of real coefficients containing bk's % A = K by 3 matrix of real coefficients containing ak's % x = input sequence [K,L] = size(B); N = length(x); w = zeros(K+1,N); w(1,:)=filter(C,1,x) ; for i = 1:1:K

w(i+1,:) = filter(B(i,:),A(i,:),x) ; end

y = sum(w);

clc;clear; n=[0: 20]; b=[3 5 4 4]; x=[n>=0];

a=[1 -0.2 0.5 0.1];

[C, B, A]=dir2par(b, a); y1=filter(b, a, x);

y2=parfiltr(C, B, A, x); subplot(2,1,1); stairs(n, y1); grid;

subplot(2,1,2); stairs(n, y2); grid;

5. （10分）下图是一个由理想运算放大器组成的电路，已知R11，R22，

C0.5F，电压ui和uo分别表示输入量和输出量。求：

（1）系统的传递函数；（5分）

（2）分别求其脉冲响应、阶跃响应的曲线；（5分） 解：

(1)

H(s)=-2/(1+s);

(2)

clc;clear; num=[-2]; den=[1 1]; G=tf(num,den); subplot(2, 1, 1); impulse(G); subplot(2, 1, 2); step(G);