Matlab大作业

2008 级

《MATAB语言与控制系统仿真》课程

大作业

姓 名 学 号 所在院系 班 级 日 期 同组人员

作业评分 评阅人

设计报告评分表

项 目 方案论证9分 设计过程应包括的主要内容或考核要点 性能指标分析；控制方法及实现方案 常见问题 未分析指标 没有说明方案理由 未进行方案比较 有分析指标，未用程序分析 有方案说明，未用程序说明 有设计过程，未运用程序 有简单程序应用，无算法 有程序，没有程序说明 未说明设计理由 扣分 -5 -1 -1 -1 -1 -20 -5 -5 -2 -9 -3 -2 -2 -2 -2 -2 -1 -1 -1 -1 -1 -9 -4 得分 控制器设计与参数计算 基本35分 部分 结果63分 分析10分 格式规范9分 第1没有误差分析 对设计结果的分析与核没有对结果验证 有量化指标，但没有分析 算，分析原因和改进 有分析验证，没有运用程序 没有总结或总结太虚 图号图名等问题 没有对图表说明 重点考查完整性，图表，中英文参杂，不一致 字体不一致 公式的规范性 图形截屏 参考文献规范问题 参考文献的引用问题 提出改进的性能指标，完成分析，设计并对结果予以验证 考虑参数变化，干扰影响未提出其他方案 有，未用程序验证 有程序，未对程序说明 提高项9分 第29 -2 -9 -4 部分 项27分 分 第3项9分 未对误差干扰进行分析 等其他因素，完成分析，有，未用程序验证 有程序，未对程序说明 设计并对结果予以验证 提出其他更完善的性能指标，完成分析，设计并对结果予以验证 未提出其他指标要求 有，未用程序验证 有程序，未对程序说明 -2 -9 -4 -2 报告 得分 90分 合计 特色报告的特色和难度系数，掌握程度予以评价 加分 0-30分 总分 报告得分+答辩/特色加分 2

目录

一、引言....................................................4

二、设计方法................................................5

三、结果分析...............................................18

四、深入探讨...............................................20

五、总结...................................................27

六、致谢和参考文献.........................................29

七、附录（程序）...........................................30

3

铣床（the milling machine）的系统设计

一、引言

本学期开始学习《MATLAB语言、控制系统分析与设计》，其实最早接触matlab是在学习概率论时。在老师的讲解与自己平时的练习中，不断地使用和学习matlab软件。其实一开始学习的目标就不局限于控制系统领域，希望了解该软件的各种用途。当然了，最基本的用途是绘图。经过半个学期的汲取，觉得它用途广泛，功能强大。我最喜欢的就是它的仿真（simulink），也许是对于初学者来说最简单的操作。学习matlab的同时，也同时加深了对控制系统知识的掌握。自控中的绘图可以在matlab下的响应分析的图形用户界面（ltiview）中建立。此外还有sisotool等自控方面的超强工具箱。

图10.1铣床模型

4

图10.2铣床系统方框图

铣床系统的开环传递函数

G(s)?2s(s?1)(s?5) （10.1）

设计要求：

重新设计滞后校正器参数

R(s)?a2a1）对于斜坡输入

s，稳态误差小于8

2）使校正后系统的阶跃响应超调量P.O.?20%

对于（1）参数a不影响系统的设计，在以下分析中作为‘1’处理

本章的核心是滞后校正系统，但之所以对铣床系统采用滞后校正，理由见（四）深入探讨中频域下的分析。另外我们超前校正器、超前-滞后校正器、PI校正器、PD校正器、PID校正器。

二、设计方法

以运用matlab编程为主，其他各种工具箱为辅助

根据极点、零点、增益输入原系统开环传递函数，并加入单位负反馈构成闭环系统 程序1

G=zpk([],[0,-1,-5],2);%原系统开环传递函数 sys=feedback(G,1);%单位负反馈，构成闭环 得到

G(s)?2s(s?1)(s?5) （10.2）

2(s?5.096s)?(2

sys? 第一个要求稳态误差小于a/8 程序2 tf1=10;

a=1;

t=[0:0.1:tf1]; u=a*t;%斜坡输入 y_t=lsim(sys,u,t);

0.s9?04205 ) （10.3） . 3 9 2

plot(t,y_t,'b-',t,u,'r:')

yss=y_t(length(t));

5

y_deta=t(length(t))-yss if y_deta<(a/8)%稳态误差小于a/8 disp('yes') else

disp('no') end

100009000800070006000500040003000200010000010002000300040005000600070008000900010000

图10.3原系统单位斜坡输入的稳态误差

得到原系统E(s)=2.500 no 根据闭环传递函数求原系统性能 程序3

step(sys);%加阶跃响应

6

Step Response1.4System: sys1.2Peak amplitude: 1.04Overshoot (%): 3.75At time (sec): 7.43System: sysSettling Time (sec): 9.6410.8udelitAmp0.60.40.20051015Time (sec)图10.4原系统阶跃响应

第二个要求超调量小于20%，则根据阻尼比与超调量的函数关系

???P.O.?100*e1??2 (10.4)

程序4

sigmaN=0.2;%要求超调上限 A=log(1/sigmaN);

zeta=sqrt((A^2)/(pi^2+A^2)); 得到阻尼比??0.4559 引入滞后校正系统

Gc(s)?Ks?z?s?p,p?z??z,

p 先考虑要求（1）

limsEas?0(s)?limsG?ac(s)G(s)8s?0 ?limsGc(s)G(s)?8 ?K2?K??z??2??z?vcomp?5??????5K????p??8????p??? 7

(10.5)

(10.6)

(10.7)

(10.8)

?1?K2s(s?1)(s?5)?0 (10.9)

画出其根轨迹

并且找到满足??0.4559的点， 程序5

n=[2];d=[1 6 5 0]; rlocus(n,d),hold on grid

z=0.4559; hs=12*z;

hc=12*sqrt(1-z^2);

plot([0 -hs],[0 hc],'--',[0 -hs],[0 -hc],'--') text(-9.5,1,'Desired performance region')

会显示Select a point in the graphics window之后在图上找符合条件的点

Root Locus150.86100.945Imaginary Axis0.760.640.50.340.160.98520017.51512.510Desired performance region7.552.50.985-50.94-100.860.76-15-20-150.64-100.50.34-5Real Axis0.160510

图10.5原系统的根轨迹来求取K的临界值

为了精确性采用Data Cursor，得到闭环极点p=-0.401+0.782i

K?2.02

8

Root Locus0.4580.785System: sysGain: 2.02Pole: -0.401 + 0.782iDamping: 0.456Overshoot (%): 20Frequency (rad/sec): 0.8790.8790.4560.4550.4540.4530.4520.880.7840.7830.878Imaginary Axis0.4590.7820.8770.7810.460.780.8760.8750.779-0.404-0.403-0.402-0.401Real Axis-0.4-0.399-0.3980.874-0.397

图10.6找寻K临界值

仅当引入增益K下，系统的阶跃响应

程序6

G=zpk([],[0,-1,-5],2);%原系统开环传递函数 G2=2.02*G;%引入K/alpha

sys2=feedback(G2,1);%单位负反馈，构成闭环 step(sys2)

9

Step Response1.41.2System: sys2Peak amplitude: 1.25Overshoot (%): 25.11At time (sec): 3.78e0.8duitlpmA0.60.40.20051015Time (sec)图10.7引入增益2.02下的阶跃响应

要满足要求一则：

?K2?K??z??2??z?vcomp?5????????p??K?5????p??8?? 根据K，并取

Kvcomp?10

程序7

Kvc=10;%Kvcopm>8

Kgang=2.02;%K/alpha Kvu=2*Kgang/5;

alpha=Kvc/Kvu%alpha即z/p 得到零点与极点的倍数 ??12.3762

得到了校正后系统开环传递函数

G2(s)?4.04s?0.01238s(s?1)(s?5)(s?0.001) 看系统是否符合要求（1）

程序8

G=zpk([],[0,-1,-5],2);%原系统开环传递函数

10

(10.10)

(10.11)

p=0.001;z=alpha*p; nc=[1 z];dc=[1,p]; Gc=2.02*tf(nc,dc); G2=G*Gc;%引入K/alpha

sys2=feedback(G2,1);%单位负反馈，构成闭环 系统建立后，输入斜坡信号 程序9 tf1=10000; a=1;

t=[0:0.1:tf1]; u=a*t;%斜坡输入

y_t=lsim(sys2,u,t);

plot(t,y_t,'b-',t,u,'r:') yss=y_t(length(t));

y_deta=u(length(t))-yss

if y_deta<(a/8)%稳态误差小于a/8 disp('yes')%符合要求（1） else

disp('no')%不符合 end 得到

E(s)?y_deta?0.100 yes

100009000800070006000500040003000200010000010002000300040005000600070008000900010000

图10.8第一次校正后系统斜坡响应

接着检验其阶跃响应

程序10

G=zpk([],[0,-1,-5],2);%原系统开环传递函数 p=0.001;z=alpha*p;

11

nc=[1 z];dc=[1,p]; Gc=2.02*tf(nc,dc);

G2=G*Gc;%引入K/alpha

sys2=feedback(G2,1);%单位负反馈，构成闭环 step(sys2)

Step Response1.41.2System: sys2Peak amplitude: 1.21Overshoot (%): 21.1At time (sec): 4.181Amplitude0.80.60.40.200102030Time (sec)405060

图10.9第一次校正后的系统阶跃响应

超调21.2%，需要继续调整

令参数稍微改变，K?2.02,1.90,2.1

程序11

G=zpk([],[0,-1,-5],2);%原系统开环传递函数 p=0.001;z=alpha*p; nc=[1 z];dc=[1,p]; Gca=2.02*tf(nc,dc);

Gcb=1.90*tf(nc,dc);

Gcc=2.1*tf(nc,dc);%引入K/alpha G2a=G*Gca; G2b=G*Gcb; G2c=G*Gcc;

sys2a=feedback(G2a,1);%单位负反馈，构成闭环 sys2b=feedback(G2b,1); sys2c=feedback(G2c,1); step(sys2a,'b--') hold on

12

step(sys2b,'r:') step(sys2c,'y-')

Step Response1.41.2System: sys2bPeak amplitude: 1.2Overshoot (%): 19.6At time (sec): 4.391Amplitude0.80.60.40.20010203040Time (sec)50607080

图10.10对参数稍作调整的系统阶跃响应

当K?1.90时才能使超调降低到20%以下 再次使用求K临界值的程序验证 程序12 hold off

clg

n=[2];d=[1 6 5 0]; rlocus(n,d),hold on grid z=0.4559;

hs=12*z;

hc=12*sqrt(1-z^2);

plot([0 -hs],[0 hc],'--',[0 -hs],[0 -hc],'--') text(-9.5,1,'Desired performance region') rlocfind(n,d)

13

Root Locus0.4640.4580.4520.4460.4410.4350.790.47System: sysGain: 2.02Pole: -0.401 + 0.782iDamping: 0.456Overshoot (%): 20Frequency (rad/sec): 0.8780.785Imaginary Axis0.780.475System: sysGain: 1.98Pole: -0.403 + 0.771iDamping: 0.463Overshoot (%): 19.4Frequency (rad/sec): 0.870.875-0.42-0.4150.87-0.41-0.4050.865-0.4Real Axis-0.395-0.390.86-0.385-0.380.7750.77

图10.11在根轨迹上调整参数探索验证

可以看出K?2.02的点满足要求（2） 对K?0.1?2.0进行探索

程序13

G=zpk([],[0,-1,-5],2);%原系统开环传递函数 Kvc=10;%Kvcopm>8

y=zeros(200,1);i=0;

for K=0.1:0.1:2%K/alpha Kvu=2*K/5;

alpha=Kvc/Kvu;%alpha即z/p p=0.001;z=alpha*p; nc=[1 z];dc=[1,p];

Gc=K*tf(nc,dc);%引入K/alpha

G2=G*Gc;

sys2=feedback(G2,1);%单位负反馈，构成闭环 t=[0:0.1:19.9]; i=i+1;

y(:,i)=step(sys2,t);

end

plot(y)

legend('K=0.1','K=0.2','K=0.3','K=0.4','K=0.5','K=0.6','K=0.7','K=0.

14

8','K=0.9','K=1.0','K=1.1','K=1.2','K=1.3','K=1.4','K=1.5','K=1.6','K=1.7','K=1.8','K=1.9','K=2.0')

1.5 K=0.1K=0.2K=0.3K=0.4K=0.5K=0.6K=0.7K=0.8K=0.9K=1.0K=1.1K=1.2K=1.3K=1.4K=1.5K=1.6K=1.7K=1.8K=1.9K=2.010.50 020406080100120140160180200

图10.12不同取值下的阶跃响应

求各曲线超调量，先将时间范围扩大到200 程序14

G=zpk([],[0,-1,-5],2);%原系统开环传递函数 Kvc=10;%Kvcopm>8

y=zeros(200,1);i=0;

for K=0.1:0.1:2%K/alpha Kvu=2*K/5;

alpha=Kvc/Kvu;%alpha即z/p p=0.001;z=alpha*p;

nc=[1 z];dc=[1,p];

Gc=K*tf(nc,dc);%引入K/alpha G2=G*Gc;

sys2=feedback(G2,1);%单位负反馈，构成闭环 t=[0:1:199];%改变时间长度 i=i+1;

y(:,i)=step(sys2,t); end

plot(y)

legend('K=0.1','K=0.2','K=0.3','K=0.4','K=0.5','K=0.6','K=0.7','K=0.

15

8','K=0.9','K=1.0','K=1.1','K=1.2','K=1.3','K=1.4','K=1.5','K=1.6','K=1.7','K=1.8','K=1.9','K=2.0')

K=0.1K=0.2K=0.3K=0.4K=0.5K=0.6K=0.7K=0.8K=0.9K=1.0K=1.1K=1.2K=1.3K=1.4K=1.5K=1.6K=1.7K=1.8K=1.9K=2.0200 1.81.61.41.210.80.60.40.20 020406080100120140160180

图10.13时间轴扩大100

观察分布情况 程序15

sigma=zeros(20,20);j=0;Kgan=zeros(20,20); for i=1:1:20;

[mp,tf]=max(y(:,i));%求最大幅值点 yss=y(length(t));%求稳态值 tp=t(tf);%峰值时间 j=j+1;

Kgan(j,1)=j/10;%K/alpha

sigma(j,1)=100*(mp-yss)/yss;%超调量 end

plot(Kgan,sigma,[0 2],[20 20]) grid

xlabel('K/alpha') ylabel('Overshoot')

16

80706050Overshoot40302010000.20.40.60.811.2K/alpha1.41.61.82

图10.14K与P.O.的关系

若找到K的取值范围，问题就比较清晰 程序16

n=find(sigma(:,1)<20);%找满足要求二的K/alpha Kn=Kgan(n,1)

0.5?K?1.9

再考虑K?0.1?2.0的情况下的单位斜坡输入的稳态误差 程序17

G=zpk([],[0,-1,-5],2);%原系统开环传递函数

Kvc=10;%Kvcopm>8

y=zeros(200,1);i=0;y_deta=zeros(20,1);Kgan=zeros(20,1); for K=0.1:0.1:2%K/alpha Kvu=2*K/5;

alpha=Kvc/Kvu;%alpha即z/p p=0.001;z=alpha*p; nc=[1 z];dc=[1,p];

Gc=K*tf(nc,dc);%引入K/alpha G2=G*Gc;

sys2=feedback(G2,1);%单位负反馈，构成闭环 t=[0:1:199]; i=i+1;

Kgan(i,1)=i/10;

17

u=t;

y(:,i)=lsim(sys2,u,t); yss=y(length(t),i);

y_deta(i,1)=u(length(t))-yss if y_deta(i,1)<(1/8)%稳态误差小于a/8 disp('yes')%符合要求（1） else

disp('no')%不符合 end

end

plot(Kgan,y_deta,[0 2],[0.125 0.125])%找到误差小于a/8的直线 grid

xlabel('K/alpha') ylabel('E(s)')

0.70.60.50.4E(s)0.30.20.1000.20.40.60.811.2K/alpha1.41.61.82

图10.15K与稳态误差的关系

程序18

m=find(y_deta<0.125);%找满足要求一的K/alpha Km=Kgan(m,1) 得到

0.2?K?1.1

综合上述两个范围0.5?K?1.1

18

三、结果分析

令K?1 程序19

Kvc=10;%Kvcopm>8 Kgang=1;%K/alpha

Kvu=2*Kgang/5;

alpha=Kvc/Kvu;%alpha即z/p 将参数带入校正后系统 程序20

G=zpk([],[0,-1,-5],2);%原系统开环传递函数 p=0.001;z=alpha*p; nc=[1 z];dc=[1,p]; Gc=1*tf(nc,dc);

G2=G*Gc;%引入K/alpha

sys2=feedback(G2,1);%单位负反馈，构成闭环

求单位斜坡响应 程序21 tf1=1000; a=1;

t=[0:0.1:tf1]; u=a*t;%斜坡输入

y_t=lsim(sys2,u,t);

plot(t,y_t,'b-',t,u,'r:') grid

yss=y_t(length(t));

y_deta=u(length(t))-yss; if y_deta<(a/8)%稳态误差小于a/8 disp('yes')%符合要求（1） else

disp('no')%不符合 end

得到E(s)?y_deta?0.100?0.125，yes，校正系统为

G25c(s)?1000s?1000s?1然后求单位阶跃响应

19

(10.12)

程序22 step(sys2)

Step Response1.4System: sys2Peak amplitude: 1.1Overshoot (%): 9.7At time (sec): 7.641.21Amplitude0.80.60.40.200204060Time (sec)80100120

图10.16K=1的系统阶跃响应

再讨论K?0.5 程序23

Kvc=10;%Kvcopm>8 Kgang=0.5;%K/alpha Kvu=2*Kgang/5;

alpha=Kvc/Kvu;%alpha即z/p

G=zpk([],[0,-1,-5],2);%原系统开环传递函数 p=0.001;z=alpha*p; nc=[1 z];dc=[1,p]; Gc=Kgang*tf(nc,dc);

G2=G*Gc;%引入K/alpha

sys2=feedback(G2,1);%单位负反馈，构成闭环 tf1=1000; a=1;

t=[0:0.1:tf1];

u=a*t;%斜坡输入

y_t=lsim(sys2,u,t);

plot(t,y_t,'b-',t,u,'r:')

20

grid

yss=y_t(length(t));

y_deta=u(length(t))-yss;

if y_deta<(a/8)%稳态误差小于a/8 disp('yes')%符合要求（1） else

disp('no')%不符合 end

得到E(s)?y_deta?0.100?0.125，yes

G

c(s)?500s?251000s?1 程序24

step(sys2)

Step Response1.41.2System: sys21Peak amplitude: 1.16Overshoot (%): 16.5At time (sec): 16.8de0.8tuliAmp0.60.40.200102030405060708090Time (sec)图10.17K=0.5的系统阶跃响应

四、深入探讨

从频域下分析 先做原系统的伯德图

程序25

G=zpk([],[0,-1,-5],2);%原系统开环传递函数

21

(10.13)

bode(G);%原系统伯德图 grid

Bode Diagram50System: GGain Margin (dB): 23.5At frequency (rad/sec): 2.240)Closed Loop Stable? YesBd(e du-50itngMa-100-150-90-135System: G)gePhase Margin (deg): 65.2d(Delay Margin (sec): 3.05e -180saAt frequency (rad/sec): 0.374PhClosed Loop Stable? Yes-225-27010-210-1100101102Frequency (rad/sec)图10.18原系统伯德图

程序26

G=zpk([],[0,-1,-5],2);%原系统开环传递函数 [mag,phase,w]=bode(G);%原系统伯德图 m=find(phase>=-135); w1=w(max(m))

得到当w?0.6877rad/sec，相位等于45deg 频率穿越点w?0.374rad/sec

系统已经满足要求二，而且有23.5dB的增益裕度 若利用增益K来调节系统使其满足要求：

G(s)?2Ks(s?1)(s?5) 对于要求一：

ess(t)?0.125a?lim1/ss?01?G(s)?5a2K 22

(10.14)

(10.15)

得到K=20 程序28

G=zpk([],[0,-1,-5],2);%原系统开环传递函数 bode(G);%原系统伯德图 hold on

bode(20*G)%引入增益K grid

legend('K=1','K=20')

Bode Diagram100System: K=20Gain Margin (dB): -2.550At frequency (rad/sec): 2.24)BClosed Loop Stable? Nod( e0duitng-50aM-100-150-90 K=1-135K=20)ged( e-180sahPSystem: K=20-225Phase Margin (deg): -6.02Delay Margin (sec): 2.4-270 At frequency (rad/sec): 2.5710-210-1100Closed Loop Stable? No101102Frequency (rad/sec)图10.19增益校正前后系统伯德图

虽然K=20满足条件一可是相位裕度为负，闭环系统不稳定 所以不能用使用单一的增益校正。 由于要求一等价于

??0.4559

而且PM?100*? (10.16) 于是将要求一等价于：相位裕度PM?50deg 在引入

K0?20的条件下，设之后校正器

23

Gc(s)?1?Ts1??Ts,??1 (10.17)

采用滞后校正的理由：

因为原系统引入增益后，相位裕度不足，但低频下的相位均大于-135deg，所以有充足的裕度来βL(?)1?1?sin(-?m)实现滞后校正。相位滞后校正是1通过衰减校正前系统伯德图的-sin(-?m)幅值部分来减小增益穿越频率，这样可以使系统产生必要的相位裕度 求滞后校正系统的步骤及原理 0?T?m1T?20??20lg??(?)0? ???m ?90(10.18)

?m??m1T? (10.19)

取得校正系统最小相位

图10.20相位滞后网络的伯德图

步骤：

(1) 根据稳态误差要求确定开环增益K。绘制未校正系统的伯德图，并求出其相位裕量和增益裕量。

(2) 确定校正后系统的增益剪切频率值为0dB

?c，即要求的相位裕度下的频率，校正后系统此处应该幅

24

(3) (3) 求?值。确定原系统频率特性在

?c处幅值下降到0dB时所必需的衰减量ΔM。由等式

?M?20lg? (10.20)

(4) 选取T值。为了使滞后校正装置产生的相位滞后对校正后系统的增益剪切频率足够小，应满足，取

T?10?c处的影响

?c (10.21)

程序29 Ko=20;

G=zpk([],[0,-1,-5],2*Ko);%原系统开环传递函数 sys=tf(G);

[mag,phase,w]=bode(sys) figure(1);margin(sys) grid hold on

gama=50;%频域下要求二 pha=gama+5-180;

wgc=spline(phase,w,pha);%找校正后的频率穿越点 na=polyval(sys.num{1},j*wgc); da=polyval(sys.den{1},j*wgc); g=na/da; g1=abs(g);

h=20*log10(g1);%h的取值等于20lg（beta），才能使幅值等于0 beta=10^(h/20);%求beta T=10/wgc;betat=beta*T;

Gc=tf([T 1],[betat 1])%滞后校正系统函数 bode(G*Gc)

legend('G','G*Gc')

25

Bode DiagramGm = -2.5 dB (at 2.24 rad/sec) , Pm = -6.02 deg (at 2.57 rad/sec)200)100dB( edut0ingaM-100-200-90 G-135G*Gc)geSystem: G*Gcd( -180Phase Margin (deg): 49.6esDelay Margin (sec): 1.57ahPAt frequency (rad/sec): 0.551-225Closed Loop Stable? Yes-270 10-410-310-210-1100101102103Frequency (rad/sec)图10.21滞后校正前后系统伯德图

得到校正滞后系统传递函数

Gc(s)?18.23s?1231.8s?1 不含增益K0 程序30

sys2=feedback(Gc*sys,1);%校正后闭环传递函数 step(sys2)%求阶跃响应

26

(10.22)

Step Response1.41.21System: sys2Peak amplitude: 1.22Overshoot (%): 22.2At time (sec): 5.12Amplitude0.80.60.40.200102030Time (sec)405060

图10.22滞后校正后系统阶跃响应

没有达到要求二 误差分析：

1）编程中很多赋值只是取得一些分散的点，因而求取参数时存在误差

2）预计相位裕度虽然理论足够，但相位滞后校正器的最大相位不为零，余量5可能不足 改变相位裕度余量 程序31

Ko=20;

G=zpk([],[0,-1,-5],2*Ko);%原系统开环传递函数 sys=tf(G);

[mag,phase,w]=bode(sys) figure(1);margin(sys) grid hold on

gama=50;%频域下要求二

pha=gama+10-180;

wgc=spline(phase,w,pha);%找校正后的频率穿越点 na=polyval(sys.num{1},j*wgc); da=polyval(sys.den{1},j*wgc); g=na/da;

g1=abs(g);

h=20*log10(g1);%h的取值等于20lg（beta），才能使幅值等于0 beta=10^(h/20);%求beta

27

T=10/wgc;betat=beta*T;

Gc=tf([T 1],[betat 1])%滞后校正系统函数 bode(G*Gc)

legend('G','G*Gc')

Gc(s)?21.71s?1340.9s?1 (10.23)

Step Response1.41.2System: sys2Peak amplitude: 1.17Overshoot (%): 17.1At time (sec): 6.051Amplitude0.80.60.40.200102030Time (sec)40506070

图10.23修改后系统的阶跃响应

Bode DiagramGm = -2.5 dB (at 2.24 rad/sec) , Pm = -6.02 deg (at 2.57 rad/sec)200Magnitude (dB)1000-100-200-90G-135Phase (deg) G*Gc-180-225-270 System: G*GcPhase Margin (deg): 54.6Delay Margin (sec): 2.06At frequency (rad/sec): 0.463Closed Loop Stable? Yes10-410-310-210-1100101102103Frequency (rad/sec)

图10.24修正后系统伯德图

加入增益后传递函数为

28

Gc(s)?434.1s?20340.9s?1 (10.24)

与式（10.12）（10.13）相比较，其介于两者之间。又因为K可以去0.5和1.0之间的值，则时域和频域下的结果可以得到统一

五、总结（暨心得）

本次设计是结合了英文教材Binder10和课件Chapter10来设计的，两者的方向截然不同。英文教材虽然提及滞后校正，但其完全是基于时域下，用理论推导出传递函数，似乎只有在运用根轨迹处依赖了电脑操作，而课件里侧重于频域下编程。所以设计分析是时域，深入探讨是频域，结合着线性控制的理论。

两种方法的优缺点。频域下的滞后校正也有很多不严格的地方，一些参数求取的近似。但其更加直观和易于操作。时域下省略了转化要求参数的步骤，但时域的下的高阶系统分析比较繁琐，若采用近似模型误差很大。

铣床系统是一型系统，在斜坡输入下的稳态误差不为零。是否工程上一定要求非零的稳态误差。如果没有此类要求是否可以通过改善系统的型来消除稳态误差，更加符合要求一。当然了这一章的侧重点是滞后校正系统的设计（Lag Compensator Design）。这也与该系统的频率特性有关，充足的相位可以采用滞后校正。然而本次设计并没有考虑，滞后校正对于频带带宽的影响，如果转折频率1/T和1/βT选取的过小，会造成带宽过小，是系统响应变得缓慢。那么由滞后校正引起的响应时间的增加，应如何再不影响斜坡输入稳态误差和阶跃响应超调量的前提下，得到改善。我分别从时域和频域下求取了滞后校正器，是否可以利用状态空间的方法。如果利用状态空间的话，要求又如何转变。

学习matlab还是觉得只是纯粹的数学问题，虽然每一章都有实例，但最根本的还是围绕着传递函数来解决问题的。不理解系统原理的话，仍然可以掌握每一章节的精髓。个人比较喜欢sisotool，虽然用的很少，但对于根轨迹和零极点配置的理解和掌握有很大帮助。这次也用到了借助根轨迹来限制阻尼比，也使得自控上的这方面知识得到应用。

六、致谢和参考文献

《线性控制系统工程》（《Linear Control Systems Engineering》）

七、附录（程序）

程序1：

29

G=zpk([],[0,-1,-5],2);%原系统开环传递函数 sys=feedback(G,1);%单位负反馈，构成闭环

程序2 tf1=10; a=1;

t=[0:0.1:tf1]; u=a*t;%斜坡输入

y_t=lsim(sys,u,t);

plot(t,y_t,'b-',t,u,'r:') yss=y_t(length(t));

y_deta=t(length(t))-yss

if y_deta<(a/8)%稳态误差小于a/8 disp('yes') else

disp('no') end 程序3

step(sys);%加阶跃响应

程序4

sigmaN=0.2;%要求超调上限 A=log(1/sigmaN);

zeta=sqrt((A^2)/(pi^2+A^2));

程序5

n=[2];d=[1 6 5 0]; rlocus(n,d),hold on grid

z=0.4559;

hs=12*z;

hc=12*sqrt(1-z^2);

plot([0 -hs],[0 hc],'--',[0 -hs],[0 -hc],'--') text(-9.5,1,'Desired performance region')

程序6

G=zpk([],[0,-1,-5],2);%原系统开环传递函数 G2=2.02*G;%引入K/alpha

sys2=feedback(G2,1);%单位负反馈，构成闭环 step(sys2)

30

程序7

Kvc=10;%Kvcopm>8 Kgang=2.02;%K/alpha Kvu=2*Kgang/5;

alpha=Kvc/Kvu%alpha即z/p 程序8

G=zpk([],[0,-1,-5],2);%原系统开环传递函数 p=0.001;z=alpha*p; nc=[1 z];dc=[1,p]; Gc=2.02*tf(nc,dc);

G2=G*Gc;%引入K/alpha

sys2=feedback(G2,1);%单位负反馈，构成闭环

程序9

tf1=10000; a=1;

t=[0:0.1:tf1]; u=a*t;%斜坡输入

y_t=lsim(sys2,u,t);

plot(t,y_t,'b-',t,u,'r:') yss=y_t(length(t));

y_deta=u(length(t))-yss if y_deta<(a/8)%稳态误差小于a/8 disp('yes')%符合要求（1） else

disp('no')%不符合 end

程序10

G=zpk([],[0,-1,-5],2);%原系统开环传递函数 p=0.001;z=alpha*p; nc=[1 z];dc=[1,p]; Gc=2.02*tf(nc,dc);

G2=G*Gc;%引入K/alpha

sys2=feedback(G2,1);%单位负反馈，构成闭环 step(sys2) 程序11

G=zpk([],[0,-1,-5],2);%原系统开环传递函数 p=0.001;z=alpha*p; nc=[1 z];dc=[1,p];

31

Gca=2.02*tf(nc,dc); Gcb=1.90*tf(nc,dc);

Gcc=2.1*tf(nc,dc);%引入K/alpha G2a=G*Gca; G2b=G*Gcb; G2c=G*Gcc;

sys2a=feedback(G2a,1);%单位负反馈，构成闭环 sys2b=feedback(G2b,1); sys2c=feedback(G2c,1); step(sys2a,'b--') hold on

step(sys2b,'r:') step(sys2c,'y-')

程序12 hold off

clg

n=[2];d=[1 6 5 0]; rlocus(n,d),hold on grid

z=0.4559;

hs=12*z;

hc=12*sqrt(1-z^2);

plot([0 -hs],[0 hc],'--',[0 -hs],[0 -hc],'--') text(-9.5,1,'Desired performance region') rlocfind(n,d)

程序13

G=zpk([],[0,-1,-5],2);%原系统开环传递函数 Kvc=10;%Kvcopm>8 y=zeros(200,1);i=0; for K=0.1:0.1:2%K/alpha Kvu=2*K/5;

alpha=Kvc/Kvu;%alpha即z/p p=0.001;z=alpha*p; nc=[1 z];dc=[1,p];

Gc=K*tf(nc,dc);%引入K/alpha G2=G*Gc;

sys2=feedback(G2,1);%单位负反馈，构成闭环 t=[0:0.1:19.9]; i=i+1;

y(:,i)=step(sys2,t);

32

end plot(y)

legend('K=0.1','K=0.2','K=0.3','K=0.4','K=0.5','K=0.6','K=0.7','K=0.8','K=0.9','K=1.0','K=1.1','K=1.2','K=1.3','K=1.4','K=1.5','K=1.6','K=1.7','K=1.8','K=1.9','K=2.0') 程序14

G=zpk([],[0,-1,-5],2);%原系统开环传递函数 Kvc=10;%Kvcopm>8

y=zeros(200,1);i=0;

for K=0.1:0.1:2%K/alpha Kvu=2*K/5;

alpha=Kvc/Kvu;%alpha即z/p p=0.001;z=alpha*p;

nc=[1 z];dc=[1,p];

Gc=K*tf(nc,dc);%引入K/alpha

G2=G*Gc;

sys2=feedback(G2,1);%单位负反馈，构成闭环 t=[0:1:199];%改变时间长度 i=i+1;

y(:,i)=step(sys2,t); end

plot(y)

legend('K=0.1','K=0.2','K=0.3','K=0.4','K=0.5','K=0.6','K=0.7','K=0.8','K=0.9','K=1.0','K=1.1','K=1.2','K=1.3','K=1.4','K=1.5','K=1.6','K=1.7','K=1.8','K=1.9','K=2.0') 程序15

sigma=zeros(20,20);j=0;Kgan=zeros(20,20); for i=1:1:20;

[mp,tf]=max(y(:,i));%求最大幅值点 yss=y(length(t));%求稳态值 tp=t(tf);%峰值时间 j=j+1;

Kgan(j,1)=j/10;%K/alpha

sigma(j,1)=100*(mp-yss)/yss;%超调量 end

plot(Kgan,sigma,[0 2],[20 20]) grid

xlabel('K/alpha') ylabel('Overshoot')

33

程序16

n=find(sigma(:,1)<20);%找满足要求二的K/alpha Kn=Kgan(n,1)

程序17

G=zpk([],[0,-1,-5],2);%原系统开环传递函数 Kvc=10;%Kvcopm>8

y=zeros(200,1);i=0;y_deta=zeros(20,1);Kgan=zeros(20,1); for K=0.1:0.1:2%K/alpha Kvu=2*K/5;

alpha=Kvc/Kvu;%alpha即z/p p=0.001;z=alpha*p; nc=[1 z];dc=[1,p];

Gc=K*tf(nc,dc);%引入K/alpha

G2=G*Gc;

sys2=feedback(G2,1);%单位负反馈，构成闭环 t=[0:1:199]; i=i+1;

Kgan(i,1)=i/10;

u=t;

y(:,i)=lsim(sys2,u,t); yss=y(length(t),i);

y_deta(i,1)=u(length(t))-yss if y_deta(i,1)<(1/8)%稳态误差小于a/8 disp('yes')%符合要求（1） else

disp('no')%不符合 end end

plot(Kgan,y_deta,[0 2],[0.125 0.125])%找到误差小于a/8的直线 grid

xlabel('K/alpha') ylabel('E(s)')

程序18

m=find(y_deta<0.125);%找满足要求一的K/alpha Km=Kgan(m,1) 程序19

Kvc=10;%Kvcopm>8 Kgang=1;%K/alpha Kvu=2*Kgang/5;

34

alpha=Kvc/Kvu;%alpha即z/p

程序20

G=zpk([],[0,-1,-5],2);%原系统开环传递函数 p=0.001;z=alpha*p; nc=[1 z];dc=[1,p]; Gc=1*tf(nc,dc);

G2=G*Gc;%引入K/alpha

sys2=feedback(G2,1);%单位负反馈，构成闭环

程序21 tf1=1000; a=1;

t=[0:0.1:tf1]; u=a*t;%斜坡输入

y_t=lsim(sys2,u,t); plot(t,y_t,'b-',t,u,'r:') grid

yss=y_t(length(t));

y_deta=u(length(t))-yss;

if y_deta<(a/8)%稳态误差小于a/8 disp('yes')%符合要求（1） else

disp('no')%不符合 end

程序22

step(sys2)

程序23

Kvc=10;%Kvcopm>8 Kgang=0.5;%K/alpha Kvu=2*Kgang/5;

alpha=Kvc/Kvu;%alpha即z/p

G=zpk([],[0,-1,-5],2);%原系统开环传递函数 p=0.001;z=alpha*p; nc=[1 z];dc=[1,p]; Gc=Kgang*tf(nc,dc); G2=G*Gc;%引入K/alpha

sys2=feedback(G2,1);%单位负反馈，构成闭环 tf1=1000; a=1;

35

t=[0:0.1:tf1]; u=a*t;%斜坡输入

y_t=lsim(sys2,u,t);

plot(t,y_t,'b-',t,u,'r:') grid

yss=y_t(length(t)); y_deta=u(length(t))-yss; if y_deta<(a/8)%稳态误差小于a/8 disp('yes')%符合要求（1） else

disp('no')%不符合 end 程序24 step(sys2)

程序25

G=zpk([],[0,-1,-5],2);%原系统开环传递函数 bode(G);%原系统伯德图 grid

程序26

G=zpk([],[0,-1,-5],2);%原系统开环传递函数 [mag,phase,w]=bode(G);%原系统伯德图 m=find(phase>=-135); w1=w(max(m))

程序28

G=zpk([],[0,-1,-5],2);%原系统开环传递函数 bode(G);%原系统伯德图 hold on

bode(20*G)%引入增益K grid

legend('K=1','K=20')

程序29

Ko=20;

G=zpk([],[0,-1,-5],2*Ko);%原系统开环传递函数 sys=tf(G);

[mag,phase,w]=bode(sys) figure(1);margin(sys) grid

36

hold on

gama=50;%频域下要求二

pha=gama+5-180;

wgc=spline(phase,w,pha);%找校正后的频率穿越点 na=polyval(sys.num{1},j*wgc); da=polyval(sys.den{1},j*wgc); g=na/da;

g1=abs(g);

h=20*log10(g1);%h的取值等于20lg（beta），才能使幅值等于0 beta=10^(h/20);%求beta T=10/wgc;betat=beta*T;

Gc=tf([T 1],[betat 1])%滞后校正系统函数 bode(G*Gc)

legend('G','G*Gc')

程序30

sys2=feedback(Gc*sys,1);%校正后闭环传递函数 step(sys2)%求阶跃响应 程序31 Ko=20;

G=zpk([],[0,-1,-5],2*Ko);%原系统开环传递函数 sys=tf(G);

[mag,phase,w]=bode(sys) figure(1);margin(sys) grid

hold on

gama=50;%频域下要求二 pha=gama+10-180;

wgc=spline(phase,w,pha);%找校正后的频率穿越点 na=polyval(sys.num{1},j*wgc); da=polyval(sys.den{1},j*wgc); g=na/da; g1=abs(g);

h=20*log10(g1);%h的取值等于20lg（beta），才能使幅值等于0 beta=10^(h/20);%求beta

T=10/wgc;betat=beta*T;

Gc=tf([T 1],[betat 1])%滞后校正系统函数 bode(G*Gc)

legend('G','G*Gc')

37

本文来源：https://www.bwwdw.com/article/w746.html