自动控制上机实验

wd=wn*sqrt(1-zeta^2) %计算阻尼振荡频率wd a=atan(sqrt(1-zeta^2)/zeta) tr=(3.14-a)/wd %计算上升时间 tp=3.14/wd%计算峰值时间

overshoot=exp(-3.14*zeta/sqrt(1-zeta^2)) %计算超调量 ts1=3/(zeta*wn)%计算过渡过程时间误差带为5% ts2=4/(zeta*wn)%误差带为2% N1=ts1/(2*tp)%振荡次数误差带为5% N2=ts2/(2*tp)%误差带为2%

结果：wd =4,a =0.9273,tr =0.5532,tp = 0.7850,overshoot = 0.0949,ts1 = 1,ts2= 1.3333,N1 =0.6369,N2 =0.8493

法二：直接从响应曲线读取：

wn25有题设可得出该系统的闭环传递函数为：Gs2 22s6s25s2wnswnMatlab程序为：n=25

d=[1 6 25] g=tf(n,d) step(g)

2Step Response1.4System: gPeak amplitude: 1.09Overshoot (%): 9.47At time (sec): 0.7911.2System: g1Rise Time (sec): 0.3710.8System: gSettling Time (sec): 1.19Amplitude0.60.40.2000.20.40.60.81Time (sec)1.21.41.61.82

习题3-13已知控制系统的特征方程为 （1）s2ss2s10

（2）s2s8s12s20s16s160 试分析系统的稳定性。

解：由matlab求解方程的特征根： （1）d=[1 2 1 2 1] roots(d)

ans =

-1.8832 0.2071 + 0.9783i 0.2071 - 0.9783i -0.5310

即该特征方程有两个正实部特征根，则系统不稳定。 （2）d=[1 2 8 12 20 16 16]

roots(d) ans =

0.0000 + 2.0000i 0.0000 - 2.0000i -1.0000 + 1.0000i -1.0000 - 1.0000i -0.0000 + 1.4142i -0.0000 - 1.4142i

65432432即该特征方程有两对纯虚数根，则系统不稳定（临界稳定）。 习题3-16试分析3-16图（a）、（b）所示系统的稳定性。 解：(a) 法一：绘制系统零极点分布图 g1=tf(10,[1 1]); %子系统g1 h1=tf([10 1],1); %子系统h1

gc=feedback(g1,h1); %得到闭环系统传递函数 pzmap(gc) %绘制系统零极点分布图

Pole-Zero Map10.80.60.40.20-0.2-0.4-0.6-0.8-1-0.12Imaginary Axis (seconds-1)-0.1-0.08-0.06Real Axis (seconds)-1-0.04-0.020

法二：g1=tf(10,[1 1]);%子系统g1 h1=tf([10 1],1);%子系统h1

gc=feedback(g1,h1);%得到闭环系统传递函数 p=eig(gc)%求系统的特征根 结果：p =-0.1089

分析：由上图和系统特征根知，由于特征根全部在s平面的左半平面，所以此负反馈系统是稳定的。

(b) g1=tf(10,[1 1]);%子系统g1 h1=tf([2 0],1);%子系统h2

g2=feedback(g1,h1);%副回路闭环传递函数 g3=tf([1 1],[1 0]);%子系统g3

g4=feedback(g2*g3,1);%主回路闭环传递函数 pzmap(g4);%绘制系统零极点分布图

Pole-Zero Map0.80.60.4Imaginary Axis (seconds-1)0.20-0.2-0.4-0.6-0.8-1-0.9-0.8-0.7-0.6-0.5-0.4-1-0.3-0.2-0.10Real Axis (seconds)

分析：由上图知，由于特征根全部在s平面的左半平面，所以此负反馈系统是稳定的。 5.8.7 135页 %观察原系统 t=0:0.01:50; n=1;

d=conv([1.4 1 0],[0.002 1]);

G0=tf(n,d);%系统固有部分的传递函数

step(feedback(G0,1),t);%单位负反馈系统单位阶跃响应

Step Response1.41.21Amplitude0.80.60.40.2005101520253035404550Time (seconds)

观察原系统响应曲线可见系统存在较大的稳态误差。 %等幅振荡法

%求取系统临界稳定时的参数，做系统根轨迹图 num=1;

den=conv([1.4 1 0],[0.002 1]); G0=tf(num,den); rlocus(G0);

Root Locus15001000System: G0Gain: 31.8Pole: -0.334 - 4.76iDamping: 0.07Overshoot (%): 80.2Frequency (rad/s): 4.77Imaginary Axis (seconds-1)5000-500-1000-1500-2000-1500-1000-500Real Axis (seconds-1)05001000

受控对象根轨迹图

由上图可得原系统在临界稳定时，Kp32，P2/4.771.32。

%施加不同比例控制效果 n=[30:10:100];

d=conv([1.4 1 0],[0.002 1]); for ii=1:8

G0=tf(n(ii),d);%施加不同比例控制时的系统开环传递函数 step(feedback(G0,1));%系统的单位阶跃响应曲线 hold on;%保持，循环绘制曲线全部位于一图上 end

Step Response21.81.61.41.2Amplitude10.80.60.40.200246810Time (seconds)1214161820

由比例控制曲线，取Kp=80 %不同微分控制效果 Kp=80;

tau=[0.1 0.125 1];

d=conv([1.4 1 0],[0.002 1]); figure; for ii=1:3

G0=tf([Kp*tau(ii),Kp],d); step(feedback(G0,1)); hold on; end hold off;

gtext('tau=0.1'); gtext('tau=0.125'); gtext('tau=1');

Step Response1.4tau=0.1tau=0.1251.21tau=1Amplitude0.80.60.40.2000.20.40.60.811.21.41.61.8Time (seconds)

由比例微分控制曲线，取tau=0.125，即取Kd=10 %不同积分控制效果 Kp=80; Kd=10;

Ki=[ 0.001 5 50];

d=conv([1.4 1 0],[0.002 1]); n=1;

G0=tf(n,d); for ii=1:3

Gc=tf([Kd Kp Ki(ii)],[1 0]);%受控对象 G=G0*Gc;%开环传递函数 step(feedback(G,1)); hold on; end hold off; gtext('Ki=0.1'); gtext('Ki=5'); gtext('Ki=10');

Step Response1.31.251.21.151.11.0510.950.90.850.20.30.40.50.60.70.80.911.11.2Ki=10Ki=5Ki=0.1AmplitudeTime (seconds)

由积分控制曲线，取Ki=5

总结分析：比例控制会导致系统的相对稳定性变差，甚至不稳定。当增大开环放大倍数来改善系统稳态性能的同时，也牺牲了系统的相对稳定性。综上，PID通过积分作用消除误差，而微分作用降低超调量、加快系统响应速度，综合了PI和PD控制各自的长处。实际工程中，PID控制器被广泛运用。 P137 5-9-1

delta=6;%最大超前相位调节参数 s=tf('s');

G=20/(s*(0.5*s+1));%得到原系统传递函数 margin(G);%原系统Bode图

[gm,pm]=margin(G)%原系统的相位裕量和幅值裕量 phim1=50;%期望相位裕度

phim=phim1-pm+delta;%需补偿的相位裕度 phim=phim*pi/180;%phim单位转换

alfa=(1+sin(phim))/(1-sin(phim));%求取校正参数a a=10*log10(alfa);%校正器在最大超前相位处的增益 [mag,phase,w]=bode(G);%返回Bode图参数 adB=20*log10(mag);%原系统幅值单位转换 wm=spline(adB,w,-a);%得到最大超前相位处的频率 t=1/(wm*sqrt(alfa));%求取校正参数t Gc=(1+alfa*t*s)/(1+t*s)%得到补偿后的校正器

[gmc,pmc]=margin(G*Gc)%求取校正后的系统稳定裕度参数 figure;

margin(G*Gc)%得到校正后的系统Bode图

Bode DiagramGm = Inf dB (at Inf rad/s) , Pm = 18 deg (at 6.17 rad/s)50System: GFrequency (rad/s): 6.24Magnitude (dB): -0.1930Magnitude (dB)Phase (deg)-50-90-135System: GFrequency (rad/s): 5.87Phase (deg): -161-18010-110010Frequency (rad/s)1102

校正前系统Bode图

Bode DiagramGm = Inf dB (at Inf rad/s) , Pm = 50.6 deg (at 8.95 rad/s)50System: untitled1Frequency (rad/s): 9.04Magnitude (dB): -0.122Magnitude (dB)0-50-100-90System: untitled1Frequency (rad/s): 8.68Phase (deg): -129-135Phase (deg)-18010-1100101102103Frequency (rad/s)

delta=6时校正后系统的Bode图

运行结果

%校正前的稳定裕度参数 gm = Inf pm =17.9642 %补偿后的校正器 Transfer function: 0.2293 s + 1 ------------- 0.05445 s + 1

%校正后的稳定裕度

gmc =Inf pmc = 50.6324 即相位裕度符合要求。 P140 例5-10-1

delta=6;%最大超前相位调节参数 s=tf('s');

G=5/(s*(s+1)*(0.5*s+1));%得到原系统传递函数 figure(1); margin(G); figure(2);

step(feedback(G,1));%查看原系统的单位阶跃响应 ex_pm=40;%期望相位裕度

phi=-180+ex_pm+delta;%期望幅值穿越频率处的相位 [mag,phase,w]=bode(G);%由bode函数返回系统参数 wc=spline(phase,w,phi);%得到期望幅值穿越频率 mag1=spline(w,mag,wc);%期望剪切频率处的原系统幅值 magdB=20*log10(mag1);%幅值单位转换 beta=10^(-magdB/20);%求得校正器参数beta t=1/(beta*(wc/10));%求得校正器参数t Gc=(1+beta*t*s)/(1+t*s);%得到校正器模型 figure(3)

margin(Gc*G)%查看校正后的Bode图 figure(4)

step(feedback(Gc*G,1))%查看校正后的阶跃响应曲线 G2=feedback(G*Gc,1)%校正后系统闭环传递函数

Bode DiagramGm = -4.44 dB (at 1.41 rad/s) , Pm = -13 deg (at 1.8 rad/s)10050Magnitude (dB)Phase (deg)0-50-100-150-90-135-180-225-27010-210-1100101102Frequency (rad/s)

校正前bode图（不稳定）

264321x 10Step ResponseAmplitude0-1-2-3-4-5050100150200Time (seconds)250300350400

校正前系统的阶跃响应曲线

Bode DiagramGm = 12.7 dB (at 1.36 rad/s) , Pm = 40.9 deg (at 0.549 rad/s)10050Magnitude (dB)Phase (deg)0-50-100-150-90-135-180-225-27010-410-310-210-1100101102Frequency (rad/s)

校正后系统bode图

System: untitled1Peak amplitude: 1.33Overshoot (%): 32.8At time (seconds): 5.26Step Response1.41.2System: untitled1Settling time (seconds): 25.61Amplitude0.80.60.40.20010203040506070Time (seconds)

校正后系统的阶跃响应曲线

习题P154 5-5

(9)g=tf(50*[0.6 1],[4 1 0 0]); bode(g);

Bode Diagram150100Magnitude (dB)Phase (deg)500-50-100-180-210-24010-210-1100101102Frequency (rad/s)

（10）g=tf(7.5*conv([0.2 1],[1 1]),[1 16 100 0]); bode(g);

Bode Diagram4020Magnitude (dB)Phase (deg)0-20-40-600-45-9010-210-1100101102103Frequency (rad/s)

5-13

g=tf(10*[0.56 1],conv([0.1 1 0],conv([1 1],[0.028 1]))); [kg,r]=margin(g) 运行结果：

%幅值裕度 kg = 7.3289 %相位裕度 r =46.5215

习题P278 8-3

a=[-5 -1;3 -1]; b=[2;5]; c=[1 2]; d=[0];

g1=ss(a,b,c,d);%状态空间形式 g2=tf(g1)%传递函数

运行结果：

Transfer function: 12 s + 59 ------------- s^2 + 6 s + 8 8-4

a=[0 1;-2 -3]; t=sym('t'); eat=expm(a*t) x0=[1;-1]; eat*x0 %运行结果： eat =

[ 2/exp(t) - 1/exp(2*t), 1/exp(t) - 1/exp(2*t)] [ 2/exp(2*t) - 2/exp(t), 2/exp(2*t) - 1/exp(t)] ans =

1/exp(t) -1/exp(t)

8-14（6）

a=[1 3 2;0 2 0;0 1 3]; b=[2 1;1 1;-1 -1]; ca=ctrb(a,b)%求可控性矩阵 rca=rank(ca)%可控矩阵的秩

%运行结果 ca =

2 1 3 2 5 4 1 1 2 2 4 4 -1 -1 -2 -2 -4 -4 rca = 2

分析：由于可控矩阵的秩为2，小于3，则不可控。 8-16

A=[a 0 0 0;0 b 0 0;0 0 c 0;0 0 0 d]; B=[0;1;0;1]; C=[0 0 1 0];

D=[C;C*A;C*A^2]%可观测性矩阵 r=rank(D)%可观测性矩阵的秩

D = 0 0 1 0 0 0 3 0 0 0 9 0

r = 1

%运行结果

分析：由于可观测性矩阵的秩为1，小于3，则不可观；当a、b、c、d互不相等，由于A是对角线标准型，c中有全零的列，故不可观。 8-24

a=[1 0.1;0 1]; b=[0.005;0.1]; p=[0.6;0.8];

k=place(a,b,p)%状态反馈矩阵

%运行结果

k =8.0000 5.6000

P285 图9-2-1开环系统和闭环系统的稳定误差

图9-2-2系统的稳定性

图9-3-1系统时域响应仿真框图