控制新版系统计算机仿真作业.docx

兰 州 理 工 大 学

《控制系统计算机仿真》
上机汇报Ⅰ

院系： 电气工程和信息工程学院
班级： 14级自动化3班
姓名： 孙悦
学号：
时间： 年 6 月 15 日

电气工程和信息工程学院

《控制系统计算机仿真》上机试验任务书Ⅰ（）

一、上机试验内容及要求
1. matlab软件
要求利用课余时间熟悉掌握matlab软件基础数值运算、基础符号运算、基础程序设计方法及常见图形命令操作。
2.各章节仿真试验内容及要求
具体试验内容及要求请详见上机试验汇报。

二、上机试验时间安排及相关事宜
1. 依据课程教学纲领要求，上机试验课时共16课时，学生须在每次上机之前做好对应准备工作，以确保在有限机时内完成仿真试验要求内容；
2. 试验完成后按要求完成相关仿真试验汇报；
3. 仿真试验汇报请按相关样本制作并A4打印，侧面装订，作为成绩评定一部分。

自动化系《控制系统计算机仿真》课程组
3月

一、Matlab基础操作
1-1用MATLAB语言求下列系统状态方程、传输函数、零极点增益和部分分式形式模型参数，并分别写出其对应数学模型表示式：
（1）
程序以下:
num=[7,24,24]

den=[10,35,50,24]

[A,B,C,D]=tf2ss(num,den)
系统状态方程：
A =

-3.5000 -5.0000 -2.4000
1.0000 0 0
0 1.0000 0

B =

1
0
0

C =

0.7000 2.4000 2.4000

D =

0
零极点增益形式：
[Z,P,K]=tf2zp(num,den)

Z =

-1.7143 + 0.6999i
-1.7143 - 0.6999i

P =

-1.2973 + 0.9838i
-1.2973 - 0.9838i
-0.9053

K =

0.7000
部分分式：
[R,P,H]=residue(num,den)

R =

-0.0071 - 0.2939i
-0.0071 + 0.2939i
0.7141

P =

-1.2973 + 0.9838i
-1.2973 - 0.9838i
-0.9053

H =

[]
（2）

A=[2.25 -5 -1.25 -0.5;2.25 -4.25 -1.25 -0.25;0.25 -0.5 -1.25 -1;1.25 -1.75 -0.25 -0.75]

A =

2.2500 -5.0000 -1.2500 -0.5000
2.2500 -4.2500 -1.2500 -0.2500
0.2500 -0.5000 1.2500 -1.0000
1.2500 -1.7500 -0.2500 -0.7500

>> B=[4;2;2;0]

B =

4
2
2
0

>> C=[0 2 0 2]

C =

0 2 0 2
>> D=0

D =

0
零极点增益形式：
>> [Z,P,K]=ss2zp(A,B,C,D)

Z =

-1.0000 + 1.2247i
-1.0000 - 1.2247i
-1.5000

P =

-0.5000 + 0.8660i
-0.5000 - 0.8660i
-1.5000 + 0.0000i
-1.5000 - 0.0000i

K =

4.0000
传输函数形式：
>> num=[0 4 14 22 15]

num =

0 4 14 22 15

>> den=[1 4 6.25 5.25 2.25]

den =

1.0000 4.0000 6.2500 5.2500 2.2500
部分分式：

>> [R，P，H]=residue(num,den)

R =

4.0000
-0.0000
0.0000 - 2.3094i
0.0000 + 2.3094i

P =

-1.5000
-1.5000
-0.5000 + 0.8660i
-0.5000 - 0.8660i

H =

[]

1-2 用殴拉法matlab编程实现下列系统输出响应在上，时数值解。
，
要求保留4位小数，并将结果以图形方法和真解比较。
t=0:0.1:1

t =

Columns 1 through 9

0 0.1000 0. 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000

Columns 10 through 11

0.9000 1.0000

h=0.1;
y(1)=1;
t=0:0.1:1;
h=0.1;
y(1)=1;
for i=1:10
y(i+1)=y(i)+h*(-1*y(i));
end
plot(t,y,'r')
hold on
m=exp(-1*t)

m =

Columns 1 through 9

1.0000 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493

Columns 10 through 11

0.4066 0.3679

plot(t,m,'bo')

1-3 用四阶龙格－库塔梯形法matlab编程实现1-2题数值解，要求以图形方法经过和真值及殴拉法比较，分析其精度。
h=0.1;
y(1)=1;
t=0:0.1:1

t =

Columns 1 through 9

0 0.1000 0. 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000

Columns 10 through 11

0.9000 1.0000

for i=1:10
k1=-1*y(i)
k2=-1*(y(i)+k1*h/2)
k3=-1*(y(i)+k2*h/2)
k4=-1*(y(i)+h*k3)
y(i+1)=y(i)+(k1+2*k2+2*k3+k4)*h/6
end

k1 =

-1

k2 =

-0.9500

k3 =

-0.9525

k4 =

-0.9047

y =

1.0000 0.9048

k1 =

-0.9048

k2 =

-0.8596

k3 =

-0.8619

k4 =

-0.8187

y =

1.0000 0.9048 0.8187

k1 =

-0.8187

k2 =

-0.7778

k3 =

-0.7798

k4 =

-0.7407

y =

1.0000 0.9048 0.8187 0.7408

k1 =

-0.7408

k2 =

-0.7038

k3 =

-0.7056

k4 =

-0.6703

y =

1.0000 0.9048 0.8187 0.7408 0.6703

k1 =

-0.6703

k2 =

-0.6368

k3 =

-0.6385

k4 =

-0.6065

y =

1.0000 0.9048 0.8187 0.7408 0.6703 0.6065

k1 =

-0.6065

k2 =

-0.5762

k3 =

-0.5777

k4 =

-0.5488

y =

1.0000 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488

k1 =

-0.5488

k2 =

-0.5214

k3 =

-0.5227

k4 =

-0.4965

y =

1.0000 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966

k1 =

-0.4966

k2 =

-0.4718

k3 =

-0.4730

k4 =

-0.4493

y =

1.0000 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493

k1 =

-0.4493

k2 =

-0.4269

k3 =

-0.4280

k4 =

-0.4065

y =

Columns 1 through 9

1.0000 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493

Column 10

0.4066

k1 =

-0.4066

k2 =

-0.3862

k3 =

-0.3873

k4 =

-0.3678

y =

Columns 1 through 9

1.0000 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493

Columns 10 through 11

0.4066 0.3679

plot(t,y,'o')
hold on
m=exp(-1*t)
plot(t,m,'r*')
lea=y-m
plot(t,lea,'g')
hold off

m =

Columns 1 through 8

1.0000 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966

Columns 9 through 11

0.4493 0.4066 0.3679

lea =

1.0e-006 *

Columns 1 through 8

0 0.0820 0.1483 0. 0.2429 0.2747 0.2983 0.3149

Columns 9 through 11

0.3256 0.3315 0.3332

1-4采取matlab语言编程实现。
程序:
disp('y=')
h=1;
y=0;
for i=0:1:63
y=y+2^i;
end
disp(y);
运行结果：
y=
1.8447e+019

y=
1.8447e+019
1-5编写matlabM-函数，以实现。要求在函数中给出必需解释和说明，同时检测输入和返回变量个数。

function sum=zuoye5(m)
sum=0;
format long
for i=0:m
t=1;
for j=1:i
t=t*2;
end
sum=sum+t;
end

二、控制系统分析
2-1设经典闭环结构控制系统以下图所表示，当阶跃输入幅值时，用sp3_1.m求取输出y(t)响应。深入考虑：当反馈通道为时，怎样经过对sp3_1.m程序修改，以实现此时所对应y(t)。

y(t)
r(t)
_
2

解：程序1：
clear all;close all;
a=[0.016 0.864 3.27 3.42 1];
b=[30 25];V=2;R=20;
X0=[0 0 0 0];n=4;T0=0;Tf=10;h=0.01;
b=b/a(1);a=a/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)];
Ab=A-B*C*V;
X=X0';y=0;t=T0;
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+K2*2+K3*2+K4)/6;
y=[y,C*X];
t=[t,t(i)+h];
end
[t',y']
plot(t,y,'b-')
title('反馈系数为2时阶跃响应曲线');
hold on; grid on
程序2：
clear all;close all;
a=[0.016 0.864 3.27 3.42 1];b=[30 25];
a1=conv(a,[0.1 1]);b1=b;
b1=b1/a1(1);a1=a1/a1(1);m1=length(b1);
b=b/a(1);a=a/a(1);m=length(b);
V=1;R=20;
T0=0;Tf=15;h=0.01;
X10=[0 0 0 0 0];n1=5;
X0=[0 0 0 0];n=4;
A1=a1(2:n1+1);
A=a(2:n+1);
C1=[fliplr(b1),zeros(1,n1-m1)];
C=[fliplr(b),zeros(1,n-m)];
A1=[rot90(rot90(eye(n1-1,n1)));-fliplr(A1)];
A=[rot90(rot90(eye(n-1,n)));-fliplr(A)];
B1=[zeros(1,n1-1),1]';
B=[zeros(1,n-1),1]';
Ab=A1-B1*C1*V;
X1=X10';y1=0;t=T0;
N=round(Tf-T0)/h;
for i=1:N;
K1=Ab*X1+B1*R;
K2=Ab*(X1+h*K1/2)+B1*R
K3=Ab*(X1+h*K2/2)+B1*R;
K4=Ab*(X1+h*K3)+B1*R;
X1=X1+h*(K1+K2*2+K3*2+K4)/6;
y1=[y1,C1*X1];
t=[t,t(i)+h];
end
X=X0';y=0;t=T0;
for i=1:N
K1=A*X+B*(R-y1(i));
K2=A*(X+h*K1/2)+B*(R-y1(i));
K3=A*(X+h*K2/2)+B*(R-y1(i));
K4=A*(X+h*K3)+B*(R-y1(i));
X=X+h*(K1+K2*2+K3*2+K4)/6;
y=[y,C*X];
t=[t,ti)+h];
end
[t',y']
plot(t,y,'r-');
title('反馈通道为惯性步骤时阶跃响应曲线');
hold on;grid on

2-2 下图中，若各步骤传输函数已知为：，，，，，，，，， ；试列写链接矩阵W、W0和非零元素阵WIJ，将程序sp4_2完善后，应用此程序求输出响应曲线。

G2(s)
G5(s)
G10(s)

G4(s)
y0
_
y7
_
G1(s)
G3(s)
G6(s)
G7(s)
G8(s)
G9(s)

解：程序
clear all;close all
a=[1 0 1 0 1 1 0 1 1]';
b=[0.01 0.085 0.01 0.051 0.0067 0.15 1 0.01 0.01]';
c=[1 1 1 1 70 0.21 130 0.1 0.0044]';
d=[0 0.17 0 0.15 0 0 0 0 0]';
P=[a,b,c,d];
WIJ=[1 0 1;2 1 1;2 9 -1;3 2 1;4 3 1;
4 8 -1;5 4 1;6 5 1;6 7 -0.212;7 6 1;
8 6 1;9 7 1];
n=9;Y0=1;
Yt0=[0 0 0 0 0 0 0 0 0];
h=0.001;L1=10;T0=0;Tf=3;
nout=7;
A=diag(P(:,1));B=diag(P(:,2));
C=diag(P(:,3));D=diag(P(:,4));
m=length(WIJ(:,1));
W0=zeros(n,1);W=zeros(n,n);
for k=1:m
if(WIJ(k,2)==0)
W0(WIJ(k,1))=WIJ(k,3);
else
W(WIJ(k,1),WIJ(k,2))=WIJ(k,3);
end
end
Q=B-D*W;
R=C*W-A;V1=C*W0;
Ab=Q\R;b1=Q\V1;
Y=Yt0';y=Y(nout);t=T0;
N=round((Tf-T0)/(h*L1));
for i=1:N
for j=1:L1
K1=Ab*Y+b1*Y0;
K2=Ab*(Y+h*K1/2)+b1*Y0;
K3=Ab*(Y+h*K2/2)+b1*Y0;
K4=Ab*(Y+h*K3)+b1*Y0;
Y=Y+h*(K1+K2*2+K3*2+K4)/6;
end
y=[y,Y(nout)];
t=[t,t(i)h*L1];
end
[t',y']
plot(t,y,'r');
title('单位阶跃响应曲线');
grid on
2-3用离散相同法仿真程序sp4_3.m重求上题输出数据和曲线，并和四阶龙格－库塔法比较精度，同时分析两种方法对仿真步长选择区分。

解：程序
clear all;close all;
a=[1 0 1 0 1 1 0 1 1]';
b=[0.01 0.085 0.01 0.051 0.0067 0.15 1 0.01 0.01]';
c=[1 1 1 1 70 0.21 130 0.1 0.0044]';
d=[0 0.17 0 0.15 0 0 0 0 0]';
P=[a,b,c,d];
WIJ=[1 0 1;2 1 1;2 9 -1;3 2 1;4 3 1;
4 8 -1;5 4 1;6 5 1;6 7 -0.212;7 6 1;
8 6 1;9 7 1];
n=9;Y0=1;
h=0.001;L1=10;T0=0;Tf=3;
nout=7;
A=P(:,1);B=P(:,2);
C=P(:,3);D=P(:,4);
m=length(WIJ(:,1));
W0=zeros(n,1);W=zeros(n,n);
for k=1:m
if(WIJ(k,2)==0)
W0(WIJ(k,1))=WIJ(k,3);
else
W(WIJ(k,1),WIJ(k,2))=WIJ(k,3);
end
end

for i=1:n
if(A(i)==0);
FI(i)=1;
FIM(i)=h*C(i)/B(i);
FIJ(i)=h*h*C(i)/B(i)/2;
FIC(i)=1;FID(i)=0;
if(D(i)~=0)
FID(i)=D(i)/B(i);
else
end
else
FI(i)=exp(-h*A(i)/B(i));
FIM(i)=(1-FI(i))*C(i)/A(i);
FIJ(i)=h*C(i)/A(i)-FIM(i)*B(i)/A(i);
FIC(i)=1;FID(i)=0;
if(D(i)~=0)
FIM(i)=(1-FI(i))*D(i)/A(i)
FIJ(i)=h*D(i)/A(i)-FIM(i)*B(i)/A(i)
FIC(i)=C(i)/D(i)-A(i)/B(i);
FID(i)=D(i)/B(i);
else
end
end
end
Y=zeros(n,1);X=Y;y=0;Uk=zeros(n,1);Ub=Uk;
t=T0:h*L1:Tf;N=length(t);
for k=1:N-1
for l=1:L1
Ub=Uk;
Uk=W*Y+W0*Y0;
Udot=(Uk-Ub)/h;
Uf=2*Uk-Ub;
X=FI'.*XFIM'.*Uk+FIJ'.*Udot;
Y=FIC'.*X+FID'.*Uf;
end
y=[y,Y(nout)];
end
[t',y']
plot(t,y)

2-4求下图非线性系统输出响应y(t)，并和无非线性步骤情况进行比较。

y(t)
e(t)
r(t)=10

5
-5

P=[0.1 1 0.5 1;0 1 20 0;2 1 1 0;10 1 1 0];
WIJ=[1 0 1;1 4 -1;2 1 1;3 2 1;4 3 1];
Z=[0 0 0 0];
S=[0 0 0 0];
h=0.01;
L1=25;
n=4;
T0=0;
Tf=20;
nout=4;
Y0=10;
sp4_4;
plot(t,y,'r')
hold on
Z=[4 0 0 0];
S=[5 0 0 0];
sp4_4;
plot(t,y,'g')
A=diag(P(:,1));B=diag(P(:,2));
C=diag(P(:,3));D=diag(P(:,4));
m=length(WIJ(:,1));
W0=zeros(n,1);W=zeros(n.n);
for k=1:m
if(WIJ(k,2)==0);W0(WIJ(k,1))=WIJ(k,3);
else W(WIJ(k,1),WIJ(k,2))=WIJ(k,3);
end;
end;
for i=1:n
if(A(i,i)==0);
FI(i,i)=1;
FIM(i,i)=h*C(i,i)/B(i,i);
FIJ(i,i)=h*h*C(i,i)/B(i,i)/2;
FIC(i,i)=1;FID(i,i)=0;
if(D(i,i)~=0);
FID(i,i)=D(i,i)/B(i,i);
else
end
else
FI(i,i)=exp(-h*A(i,i)/B(i,i));
FIM(i,i)=(1-FI(i,i))*C(i,i)/A(i,i);
FIJ(i,i)=h*C(i,i)/A(i,i)-FIM(i,i)*B(i,i)/A(i,i);
FIC(i,i)=1;FID(i,i)=0;
if(D(i,i)~=0);
FIC(i,i)=C(i,i)/D(ii)-A(i,i)/B(i,i);
FID(i,i)=D(i,i)/B(i,i);
else
end
end
end
Y=zeros(n,1);X=Y;y=0;Uk=zeros(n,1);Ubb=Uk;
t=T0:h*L1:Tf;N=length(t);
for k=1:N-1
for i=1:L1
Ub=Uk;
Uk=W*Y+W0*Y0;
for i=1:n
if(Z(i)~=0)
if(Z(i)==1)
Uk(i,i)=satu(Uk(i,i),S(i));
end
if(Z(i)==2)
Uk(i,i)=dead(Uk(i,i),S(i));
end
if(Z(i)==3)
[Uk(i,i),Ubb(i,i)]=backlash(Ubb(i,i),Uk(i,i),Ub(i,i),s(i));
end
end
end
Udot=(Uk-Ub)/h;
Uf=2*Uk-Ub;
X=FI'*X+FIM'*Uk+FIJ'*Udot;
Yb=Y;
Y=FIC'*X+FID'*Uf;
for i=1;n
if(Z(i)~=0)
if(Z(i)==4)
Y(i,i)=satu(Y(i,i),S(i));
end
if(Z(i)==5)
Y(i,i)=dead(Y(i,i),S(i));
end
if(Z(i)==6)
[Y(i,i),Ubb(i,i)]=backlash(Ubb(i,i),Y(i,i),Yb(i,i),s(i));
end
end
end
end
y=[y,Y(nout)];
end
backlash函数:
function[Uc,Ubb]=backlash(Urb,Ur,Ucb,S1)
if(Ur>Urb)
if((Ur-S1)>=Ucb)
Uc=Ur-S1;
else Uc=Ucb;
end
else if(Ur函数:
function
Uc=satu(Ur,S1)
if(abs(Ur)>=S1)
if(Ur>0)
Uc=S1;
else Uc=-S1;
end
else Uc=Ur;
end
dead函数:
function
Uc=dead(Ur,S1)
if(abs(Ur)>=S1)
if(Ur>0)
Uc=Ur-S1;
else Uc=Ur+S1;
end
else Uc=0;
end

2-5 已知系统，当分别取-1、0、1时，判别系统能控性和能观察性。
a=-1时
num=[-1]
den=[10,27,18]
[A,B,C,D]=tf2ss(num,den)
P=ctrb(A,B)
N=rank(P)
Q=obsv(A,C)
N=rank(Q)
A =

-2.7000 -1.8000
1.0000 0

B =

1
0

C =

0 -0.1000

D =

0

P =

1.0000 -2.7000
0 1.0000

N =

2

Q =

0 -0.1000
-0.1000 0

N =

2
a=0时
A =

-2.7000 -1.8000
1.0000 0

B =

1
0

C =

0 0

D =

0

P =

1.0000 -2.7000
0 1.0000

N =

2

Q =

0 0
0 0

N =

0
a=1时
A =

-2.7000 -1.8000
1.0000 0

B =

1
0

C =

0 0.1000

D =

0

P =

1.0000 -2.7000
0 1.0000

N =

2

Q =

0 0.1000
0.1000 0

N =

2
2-6 已知系统状态方程为，采取状态反馈，将系统几点配置到-1、-2、-3，求状态反馈阵K。
A=[-2 -1 1;1 0 1;-1 0 1]
B=[1 1 1]'
P=[-1,-2,-3]
K=acker(A,B,P)
K =

-1 2 4