西电电院自动控制技术上机报告 本文关键词：上机，自动控制，西电，报告，技术
西电电院自动控制技术上机报告 本文简介：自动控制技术上机实验报告班级：021215学号：021214姓名：时域分析程序源代码：closeall;clearall;ft=30;M=1;B=5;K=20;%系统参数t0=0;tfinal=5;tspan=[t0tfinal];%设置仿真开始和结束时间x0=[0,0];%系统初始值，零初始条件o
西电电院自动控制技术上机报告 本文内容：
自动控制技术上机实验报告
班级：021215
学号：021214
姓名：
时域分析
程序源代码：
close
all;
clear
all;
ft
=
30;
M=1;
B=5;
K=20;
%系统参数
t0=0;
tfinal
=
5;
tspan
=
[t0
tfinal];
%设置仿真开始和结束时间
x0
=
[0,0];
%系统初始值，零初始条件
options
=
odeset(
AbsTol,[1e-6;1e-6]);
%设置仿真计算精度
[t,x]
=
ode113(
xt4odefile,tspan,x0,options);
%微分方程求解，计算位移x(:,1)和速度x(:,2)
a
=
1/M*(ft-B*x(:,2)-K*x(:,1));
%计算加速度
i
=
1;
while
(abs(a(i))>0.0001|(abs(x(i,2))>0.0001))
i
=
i+1;
end
disp(
稳态时系统的位移、速度和加速度及对应的时间分别为：
);
%显示计算结果
result
=
sprintf(
位移
d=%6.4f/n,x(i,1));
disp(result);
result
=
sprintf(
速度
v=%8.6f/n,x(i,2));
disp(result);
result
=
sprintf(
加速度
a=%9.6f/n,a(i));
disp(result);
result
=
sprintf(
时间
t=%4.2f/n,t(i));
disp(result);
d
=
x(:,1);
subplot(1,3,1),plot(t,d);
%绘制时间-位移曲线
xlabel(
时间（秒）
);ylabel(
位移（米）
);
title(
时间-位移曲线
);grid;
v
=
x(:,2);
subplot(1,3,2),plot(t,v);
%绘制时间-速度曲线
xlabel(
时间（秒）
);ylabel(
速度（米/秒）
);
title(
时间-速度曲线
);grid;
subplot(1,3,3),plot(d,v);
%绘制位移-速度曲线
xlabel(
位移（米）
);ylabel(
速度（米/秒）
);
title(
位移-速度曲线
);grid;
运行结果：
>>
EX11
稳态时系统的位移、速度和加速度及对应的时间分别为：
位移
d=1.5000
速度
v=-0.000086
加速度
a=-0.000084
时间
t=4.46
源程序代码：
close
all;
clear
all;
num
=
[0,2,5,7];
den
=
[1,6,10,6];
[z1,p1,k1]
=
tf2zp(num,den)
[r2,p2,k2]
=
residue(num,den)
运行结果：
>>
EX12
z1
=
-1.2500
+
1.3919i
-1.2500
-
1.3919i
p1
=
-3.7693
+
0.0000i
-1.1154
+
0.5897i
-1.1154
-
0.5897i
k1
=
2
r2
=
2.2417
+
0.0000i
-0.1208
-
1.0004i
-0.1208
+
1.0004i
p2
=
-3.7693
+
0.0000i
-1.1154
+
0.5897i
-1.1154
-
0.5897i
k2
=
[]
源程序代码：
clc
a=[6.3223
18
12.811];
b=[1
6
1.322
18
12.811];
sys=tf(a,b);
t=0:.005:35;
step(sys)
title(
系统的单位阶跃响应
)
[y,t]=step(a,b,t);
r=1;
while(y(r)0.98
dun=[0.25,0.01,1,0,0];
sys=tf(num,dun)
figure(1)
bode(sys)
figure(2)
sys2=feedback(sys,1)
bode(sys2)
运行结果：
sys
=
0.01
s^2
+
0.0001
s
+
0.01
--------------------------
0.25
s^4
+
0.01
s^3
+
s^2
Continuous-time
transfer
function.
sys2
=
0.01
s^2
+
0.0001
s
+
0.01
------------------------------------------------
0.25
s^4
+
0.01
s^3
+
1.01
s^2
+
0.0001
s
+
0.01
Continuous-time
transfer
function.
源程序代码：
close
all;
clear
all;
num
=
[0
20
20
10];
%开环传递函数分子
den
=
conv([1
1
0],[1
10]);
%开环传递函数的分母
nyquist(num,den)
%v
=
[-2
3
-3];
axis([-2
2
-3
3])
grid
x
=
-pi:0.01:pi;
plot(x,sin(x)),grid
on
运行结果：
源程序代码：
close
all;
clear
all;
num
=
[2000,2000];
%开环传递函数的分子
den
=
conv([1
0.5
0],[1
14
400]);
%开环传递函数的分母
nichols(num,den)
%绘制nichols图
v
=
[-270
-90
-40
40];
axis(v)
ngrid
%标出nichols图线
运行结果：
源程序代码：
num=[0
2000
2000];
den=conv([1
0.5
0],[1
14
400]);
h=tf(num,den);
[gm,pm,wg,wc]=margin(h);
gm,pm,wg,wc
运行结果：
>>
EX24
gm
=
2.7493
pm
=
73.3527
wg
=
19.8244
wc
=
5.3477
源程序代码：
clc
num=[1];
den=[0.5
1.5
1
0];
sys=tf(num,den)
sys2=feedback(sys,1)
bode(sys2)
[gm
pm
wg
wp]=margin(sys)
运行结果：
-------------------------
0.5
s^3
+
1.5
s^2
+
s
+
1
Continuous-time
transfer
function.
gm
=
3.0000
pm
=
32.6133
wg
=
1.4142
wp
=
0.7494
现代控制理论
源程序代码：
close
all;
clear
all;
%3.1_A
num
=
[1
2
3];
%传递函数分子多项式的系数
den
=
[1
3
3
1];
%传递函数分母多项式的系数
[A,B,C,D]
=
tf2ss(num,den)
%3.1_B
z
=
[-1;-3];
%传递函数的零点
p
=
[0;-2;-4;-6;];
%传递函数的极点
k
=
4;
[A,B,C,D]
=
zp2ss(z,p,k)
%3.1_C
A
=
[0,1;-1,-2];
B
=
[0;1];
C
=
[1,3];
D
=
[1];
[num,den]
=
ss2tf(A,B,C,D)
printsys(num,den,s
)
[z,p,k]
=
ss2zp(A,B,C,D)
运行结果：
>>
EX31
A
=
-3
-3
-1
1
0
0
0
1
0
B
=
1
0
0
C
=
1
2
3
D
=
0
A
=
-10.0000
-4.8990
0
0
4.8990
0
0
0
-6.0000
-4.2866
-2.0000
0
0
0
1.0000
0
B
=
1
0
1
0
C
=
0
0
0
4
D
=
0
num
=
1.0000
5.0000
2.0000
den
=
1
2
1
num/den
=
s^2
+
5
s
+
2
-------------
s^2
+
2
s
+
1
z
=
-0.4384
-4.5616
p
=
-1
-1
k
=
1
源程序代码：
close
all;
clear
all;
A1
=
[0,1;-1,-2];
B1
=
[0;1];
C1
=
[1,3];
D1
=
[1];
A2
=
[0,1;-1,-3];
B2
=
[0;1];
C2
=
[1,4];
D2
=
[0];
[A,B,C,D]
=
series(A1,B1,C1,D1,A2,B2,C2,D2)
[A,B,C,D]
=
parallel(A1,B1,C1,D1,A2,B2,C2,D2)
[A,B,C,D]
=
feedback(A1,B1,C1,D1,A2,B2,C2,D2)
[A,B,C,D]
=
feedback(A1,B1,C1,D1,A2,B2,C2,D2,+1)
运行结果：
>>
EX32
A
=
0
1
0
0
-1
-3
1
3
0
0
0
1
0
0
-1
-2
B
=
0
1
0
1
C
=
1
4
0
0
D
=
0
A
=
0
1
0
0
-1
-2
0
0
0
0
0
1
0
0
-1
-3
B
=
0
1
0
1
C
=
1
3
1
4
D
=
1
A
=
0
1
0
0
-1
-2
-1
-4
0
0
0
1
1
3
-2
-7
B
=
0
1
0
1
C
=
1
3
-1
-4
D
=
1
A
=
0
1
0
0
-1
-2
1
4
0
0
0
1
1
3
0
1
B
=
0
1
0
1
C
=
1
3
1
4
D
=
1
源程序代码：
clc;
close
all;
clear
all;
A
=
[0
-2;1
-3];
t
=
0.2;
Phi
=
expm(A*t);
%求状态转移矩阵
B
=
[2;0];
C
=
[0
3];
D
=
[0];
x0
=
[1
1];
t
=
[0
0.2];
u
=
0*t;
[y,x]
=
lsim(A,B,C,D,u,t,x0)
%求系统响应
运行结果：
y
=
3.0000
2.0110
x
=
1.0000
1.0000
0.6703
0.6703
源程序代码：
clc
A=[-3
1;1
-3];
B=[1
1;1
1];
C=[1
1;1
-1];
D=[0
0;0
0];
N=size(A);
n=N(1);
[num,den]=ss2tf(A,B,C,D,2);
disp(
可控矩阵：
)
S=ctrb(A,B)
f=rank(S)
if(f==n)
disp(
系统是可控的
)
else
disp(
系统是不可控的
)
end
disp(
)
disp(
可观测矩阵：
)
V=obsv(A,C)
m=rank(V)
if(f==n)
disp(
系统时可观测的
)
else
disp(
系统是不可观测的
)
end
运行结果：
可控矩阵：
S
=
1
1
-2
-2
1
1
-2
-2
f
=
1
系统是不可控的
可观测矩阵：
V
=
1
1
1
-1
-2
-2
-4
4
m
=
2
系统是不可观测的
源程序代码：
clc
A=[0
1;-2
-3];
B=[0;1];
C=[2
0];D=0;
P_S=[-1
-2];
k=acker(A,B,P_S);
P_O=[-3
-3];
h=(acker(A,C,P_O));
A1=[A,-B*k;h*C
A-B*k-h*C];
B1=[B;B];
C1=[C,zeros(1,2)];
D1=0;
sys=ss(A1,B1,C1,D1)
tf(sys)
运行结果：
校正设计
源程序代码：
ts=0.001;
sys=tf([400],[1,50,0]);
dsys=c2d(sys,ts,z
);
[num,den]=tfdata(dsys,v
);
u_1=0.0;u_2=0.0;
y_1=0.0;y_2=0.0;
x=[0,0,0]
;
error_1=0;
error_2=0;
for
k=1:1:1000
time(k)=k*ts;
kp=8;ki=0.1;kd=10;
rin(k)=0.5*sin(2*pi*k*ts);
du(k)=kp*x(1)+kd*x(2)+ki*x(3);
u(k)=u_1+du(k);
yout(k)=-den(2)*y_1-den(3)*y_2+num(2)*u_1+num(3)*u_2;
error(k)=rin(k)-yout(k);
u_2=u_1;u_1=u(k);
y_2=y_1;y_1=yout(k);
x(1)=error(k)-error_1;
x(2)=error(k)-2*error_1+error_2;
x(3)=error(k);
error_2=error_1;
error_1=error(k);
end
figure(1);
plot(time,rin,b,time,yout,r
),grid
on
gtext(
rin/rightarrow
)
gtext(
/leftarrowyout
)
title(
系统输出曲线
)
xlabel(
time(s)
),ylabel(
rin,yout
);
figure(2);
plot(time,error,r
),grid
on
title(
误差曲线
)
xlabel(
time(s)
);ylabel(
error
);
注：输入信号为正弦，采样信号为1ms。
运行结果：
源程序代码：
clc;
clear
all;
ts=20;
sys=tf([1],[60,1],inputdelay,80);
dsys=c2d(sys,ts,z
);
[num,den]=tfdata(dsys,v
);
u_1=0;u_2=0;u_3=0;u_4=0;u_5=0;
y_1=0;y_2=0;y_3=0;
error_1=0;
ei=0;
for
k=1:1:200
time(k)=k*ts;
rin(k)=40;
kp=0.80;
ki=0.005;
kd=3.0;
yout(k)=-den(2)*y_1+num(2)*u_5;
error(k)=rin(k)-yout(k);
M=2;
if
M==1
if
abs(error(k))>=30
end
if
u(k)=um
u(k)=um;
end
if
u(k)=num
if
error(k)>0
alpha=0;
else
alpha=1;
end
elseif
u(k)0
alpha=1;
else
alpha=0;
end
else
alpha=1;
end
elseif
M==2
alpha=1;
end
u_3=u_2;u_2=u_1;u_1=u(k);
y_3=y_2;y_2=y_1;y_1=yout(k);
error_1=error(k);
x(1)=error(k);
x(2)=(error(k)-error_1)/ts;
x(3)=x(3)+alpha*error(k)*ts;
xi(k)=x(3);
end
figure(1);
subplot(3,1,1);
plot(time,rin,b,time,yout,r
);
xlabel(
time(s)
);ylabel(
Controller?output
);
subplot(3,1,3);
plot(time,xi,r
);
xlabel(
time(s)
);ylabel(
Integration
);
运行结果：
M=2时，无抗积分饱和：
M=1时，有抗积分饱和：
源程序代码：
clc;
clear
all;
ts=0.001;
sys=tf(5.235e005,[1,87.35,1.047e004,0]);
dsys=c2d(sys,ts,z
);
[num,den]=tfdata(dsys,v
);
u_1=0;u_2=0;u_3=0;u_4=0;u_5=0;
y_1=0;y_2=0;y_3=0;
yy_1=0;
error_1=0;error_2=0;ei=0;
sys1=tf([1],[0.04,1]);
dsys1=c2d(sys1,ts,tucsin
);
[num1,den1]=tfdata(dsys1,v
);
f_1=0;
for
k=1:1:2000
time(k)=k*ts;
rin(k)=1;
yout(k)=-den(2)*y_1-den(3)*y_2-den(4)*y_3+num(2)*u_1+num(3)*u_2+num(4)*u_3;
D(k)=0.50*rands(1);
yyout(k)=yout(k)+D(k);
filty(k)=-den1(2)*f_1+num1(1)*(yyout(k)+yy_1);
error(k)=rin(k)-filty(k);
if
abs(error(k))=10
u(k)=10;
end
if
u(k)<=-10
u(k)=-10;
end
rin_1=rin(k);
u_3=u_2;u_2=u_1;u_1=u(k);
y_3=y_2;y_2=y_1;y_1=yout(k);
f_1=filty(k);
yy_1=yyout(k);
error_2=error_1;
error_1=error(k);
end
figure(1);
subplot(2,1,1);
plot(time,rin,r,time,filty,b
);
xlabel(
time(s)
);ylabel(
rin,yout
);
subplot(2,1,2);
plot(time,u,r
);
xlabel(
time(s)
);ylabel(
u
);
figure(2);
plot(time,D,r
);
xlabel(
time(s)
);ylabel(
Disturbance?signal
);
运行结果：
M=1时，不带死区结果：
M=2时，带死区结果：