基于MATLAB对控制课程的辅助应用分析.docx
《基于MATLAB对控制课程的辅助应用分析.docx》由会员分享，可在线阅读，更多相关《基于MATLAB对控制课程的辅助应用分析.docx（20页珍藏版）》请在冰豆网上搜索。
基于MATLAB对控制课程的辅助应用分析
基于MATLAB对控制课程的辅助应用分析
姓名XXX班级08自动化学号200805070029
申请得分：
95
摘要：
早期的控制系统设计可以由纸笔等工具容易地计算出来，如Ziegler与Nichols于1942年提出的PID经验公式就可以十分容易地设计出来。
随着控制理论的迅速发展，光利用纸笔以及计算器等简单的运算工具难以达到预期的效果，加之在计算机领域取得了迅速的发展，于是很自然地出现了控制系统的计算机辅助设计（computer-aidedcontrolsystemdesign,CACSD）方法。
随着MATLAB语言出现以来，就深受控制领域学生和研究者的欢迎，已经成为控制界最流行、最有影响的通用计算机语言，MATLAB作为主要程序设计语言来介绍控制系统计算机辅助设计的算法，可以使学者将主要精力集中在控制系统理论和方法上，而不是将主要精力花费在没有太大价值的底层重复性机械性劳动上，这样可以对控制系统计算机辅助设计技术有较好的整体了解，避免“只见树木，不见森林”的认识偏差，提高控制器设计的效率和可靠性。
子曰:
“工欲善其事，必先利其器”。
跟踪国际最先进的CACSD软件环境及发展，以当前国际上最流行的CACSD软件环境MATLAB为基本出发点来系统地介绍控制系统计算机辅助设计技术及软件实现，从而大大提高CACSD算法研究与实际应用的效率和可靠性。
关键词：
Matlab经典控制现代控制数学模型传递函数根轨迹稳定性劳斯判据奈氏图伯德图能控性能观性状态空间模型
正文：
控制理论的主要内容主要分为两大部分：
经典控制理论和现代控制理论。
但是将这两者总结道一起又可以详尽的分为以下6个分类：
1．数学模型
经典（时域法）
现代（频域法）
理论基础
建立在以
1.常微分方程稳定性理论
2.Fourier变换
为基础的根轨迹和奈奎斯特判据理论之上
1.常微分方程稳定性理论
2.状态空间分析
3.泛函分析、微分几何等现代数学分支
数学模型
传递函数
（研究系统外部特性,属于外部描述,不完全描述。
）
状态空间表达式
（深入系统内部,是内部描述,完全描述。
）
适用对象
仅适用于:
单输入单输出
线性
定常
集总参数
可推广至:
多输入多输出
非线性
时变
分布参数
传递函数表示方法
MATLAB语言表达
Num=[b1,b2,…，bm];
Den=[a1,a2,…,an];
G=tf（num,den）;
G=ss（A,B,C,D）
2.稳定性分析及时域分析
2.1对控制系统的性能的要求，主要是稳定性、暂态性能和稳态性能几个方面。
系通过分析是系统设计的基础，而大部分的系统设计方法都是在系统稳定性的基础上发展起来的。
常见的线性定常系统的时域分析方法有劳斯判据、赫尔维兹稳定判据等。
线性定常连续系统稳定的充要条件是系统的全部特征根或闭环极点都具有负实部，或者说位于复平面左半部。
劳斯判据不仅能够判别系统是否稳定，而且能够确定有多少正实部根，也能够具体确定对称于远点的特征根。
线性系统的动态性能取决于系统的闭环极点和零点的分布。
虽然可以通过计算机直接求解得到闭环极点，但不能看到系统闭环极点随着系统参数变化的情况。
但是通过根轨迹法，能够根据系统的开环零、极点分布，用图解的方法画出系统闭环极点随着系统参数变化的轨迹。
虽然用手工精确绘制系统更轨迹是非常困难的，但是Matlab中专门提供了绘制根轨迹有关的函数。
[r,k]=rlocus（num,den）和[r,k]=rlocus（num,den,k）的功能是绘制根轨迹图。
[r,k]=rlocus（num,den）是绘制
部分的根轨迹。
如果要以给定的参数范围绘制根轨迹，则执行命令[k,poles]=rlofind（num,den,p）。
[k,poles]=rlofind（num,den,）和[k,poles]=rlofind（num,den,p）的功能是确定根轨迹上poles处的根轨迹放大系数的值。
例：
键入：
n=[1,3];
d=[1,13,54,82,60,0];
rlocus（n,d）
可得图1所示根轨迹图：
图1
3.频域分析
当系统是高阶系统时，系统微分方程的求解时很困难的；另外，系统的时间响应没有明确反映出系统响应与系统结构、参数之间的关系，一旦系统很难满足要求，就很难确定如何去调整系统的结构和参数。
频域分析法克服了时域的不足。
根据系统的频率特性，可以直观地分析系统的稳定性，并且很容易的二和系统的结构、参数联系起来，因此可以根据系统的频率特性选择系统的结构和参数，使之满足控制要求。
实验中常用的频率分析方法是伯德图和奈氏图。
因此，可以使用Matlab绘制系统的博得图、奈氏图，并确定系统的相位裕度和幅值裕度。
用Matlab很容易精确的绘制奈氏图和伯德图，但是只能绘制其中的一部分，不能绘制出
的全貌，所以不能用来分析系统的稳定性。
3.1用Matlab绘制奈氏图
绘制奈氏图的Matlab命令是nyquist（num,den）,若要指定频率ω时，可用函数nyquist（num,den，ω）。
还有以下两种形式：
[re,im,w]=nyquist（num,den）
[re,im,w]=nyquist（num,den,w）
通过以上两种形式的调用，可计算频率特性的实部和虚部，但不能产生奈氏图，还需调用plot（re,im）函数才能得到。
例如：
开环传递函数为
用Matalb键入：
G=tf（75*[1,2],conv（[1,10],[1,3,2,5]））,nyquist（G）
得到如图2所示：
图2
3.2用Matlab绘制伯德图
绘制伯德图可用命令bode（num,den）.如果需要给出频率ω的频率范围，可调用之林w=logspace（a,b,n）,频率ω的采样点可在一定范围内产生n个十进制对数分点的等距离点。
如果需要制定幅值mag和phase范围，这执行命令[mag,phase,w]=bode（num,den），Matlab在频率响应范围内能够自动选取频率ω值绘图。
例：
用Matlab绘制伯德图
键入命令：
G=tf（2000*[1.5],conv（[1,2,0],[1,4,100]））,bode（G）
输出结果如下：
Transferfunction:
3000
-----------------------------
s^4+6s^3+108s^2+200s
图3
若键入G=tf（2000*[1.5],conv（[1,2,0],[1,4,100]））;[kg,r]=margin（G）
可以得到G的幅值域度kg=0.8296和幅值穿越频率r=-14.4062
4.结构问题
关于结构问题主要涉及现代控制理论，它包含了系统的能控能观性、以及对应的能控规范型和能观规范型、系统的能控分解和能观分解、系统的最小实现。
线性系统的能控性和能观性事基于状态方程的控制理论的基础.关于线性系统空间结构涉及的两个主要问题就是系统的状态空间模型的结构性分解以及传递函数阵与能控性/能观性的关系。
而由系统的传递函数建立状态空间模型这类问题称为系统实现问题,而求得的状态空间模型称为相应的传递函数的一个实现。
两种系统实现方法--能控/能观规范形实现。
4.1线性系统的可控性判定
基于MATLAB的判定方法：
rank（T）
构造可控性判定矩阵:
例：
离散状态方程的可控性
MATLAB求解
判定矩阵
判定矩阵构造方法
4.1.1线性系统的能观性判定
4.1.1.1判定矩阵4.1.1.2MATLAB求解
4.1.2可控标准型和能观标准型
使用下面两个转换函数时线输入如下函数sscanform（）
functionGs=sscanform（G,type）
switchtype
case'ctrl'
G=tf（G）;Gs=[];
G.num{1}=G.num{1}/G.den{1}
（1）;%传递函数归一化
G.den{1}=G.den{1}/G.den{1}
（1）;d=G.num{1}
（1）;
G1=G;G1.ioDelay=0;G1=G1-d;
num=G1.num{1};den=G1.den{1};n=length（G.den{1}）-1;
A=[zeros（n-1,1）eye（n-1）;-den（end:
-1:
2）];
B=[zeros（n-1,1）;1];C=num（end:
-1:
2）;D=d;
Gs=ss（A,B,C,D,'Ts',G.Ts,'ioDelay',G.ioDelay）;
case'obsv'
Gc=sscanform（G,'ctrl'）;
Gs=ss（Gc.a',Gc.c',Gc.b',Gc.d','Ts',G.Ts,'ioDelay',G.ioDelay）;
otherwise
error（'Onlyoptions''ctrl''and''obsv''areapplicable.'）
end
可控标准型Gs=sscanform（G,'obsv'）
例如：
对
求解可控标准型
键入命令：
num=[602810];den=[20648];
G=tf（num,den）;Gs=sscanform（G,'obsv'）
标准型：
可观测标准型
4.1.2.1Kalman规范分解
子空间示意图4
图4
4.1.2.2最小实现
例：
多变量模型
是否能最小实现？
MATLAB求解
A=[-6,-1.5,2,4,9.5;-6,-2.5,2,5,12.5;-5,0.25,-0.5,3.5,9.75;
-1,0.5,0,-1,1.5;-2,-1,1,2,3];
B=[6,4;5,5;3,4;0,2;3,1];D=zeros
（2）;
C=[2,0.75,-0.5,-1.5,-2.75;0,-1.25,1.5,1.5,2.25];
G=ss（A,B,C,D）;G1=minreal（G）
可得如下：
a=
x1x2x3
x1-2.169-1.9670.1868
x20.1554-0.1811-0.2789
x30.14031.613-1.65
b=
u1u2
x1-7.907-5.5
x23.2282.506
x32.4615.046
c=
x1x2x3
y1-0.8326-0.1501-0.04028
y2-0.38420.27640.4348
d=
u1u2
y100
y200
5．系统设计
系统设计主要是基于输入输出模型和基于状态空间模型的计。
由此可以使用Matlab设计出相应的算法直接设计，并进行整个系统的仿真分析。
5.1基于输入输出模型的设计方法
5.1.1.超前滞后校正器设计方法
串联超前滞后校正器，如图5
图5
超前滞后校正器α>1;β>1：
超前滞后校正器零极点分配图与伯德图，如图6所示
图6
5.1.2超前滞后校正器的设计方法
基于剪切频率和相位裕度的设计方法
式中
，
为校正器的增益
根据超前滞后校正器的涉及规则，可以写出相应的Matlab语言的超前滞后校正器的设计函数leadladc（），内容如下
functionGc=leadlagc（G,Wc,Gam_c,Kv,key）
G=tf（G）;[Gai,Pha]=bode（G,Wc）;
Phi_c=sin（（Gam_c-Pha-180）*pi/180）;
den=G.den{1};a=den（length（den）:
-1:
1）;
ii=find（abs（a）<=0）;num=G.num{1};G_n=num（end）;
iflength（ii）>0,a=a（ii
（1）+1）;else,a=a
（1）;end;
alpha=sqrt（（1-Phi_c）/（1+Phi_c））;Zc=alpha*Wc;Pc=Wc/alpha;
Kc=sqrt（（Wc*Wc+Pc*Pc）/（Wc*Wc+Zc*Zc））/Gai;K1=G_n*Kc*alpha/a;
ifnargin==4,key=1;
ifPhi_c<0,key=2;else,ifK1end
switchkey
case1,Gc=tf（[1Zc]*Kc,[1Pc]）;
case2
Kc=1/Gai;K1=G_n*Kc/a;Gc=tf（[10.1*Wc],[1K1*Gcn
（2）/Kv]）;
case3
Zc2=Wc*0.1;Pc2=K1*Zc2/Kv;Gcn=Kc*conv（[1Zc],[1,Zc2]）;
Gcd=conv（[1Pc],[1,Pc2]）;Gc=tf（Gcn,Gcd）;
end
例：
，
，设计超前滞后校正器。
键入命令：
s=tf（'s'）;G=4*（s+1）*（s+0.5）/s/（s+0.1）/（s+2）/（s+10）/（s+20）
wc=20;f1=figure;f2=figure;%打开两个图形窗口
forgam=20:
10:
90
Gc=leadlagc（G,wc,gam,1000,3）;
figure（f1）;step（feedback（G*Gc,1）,1）;holdon
figure（f2）;bode（Gc*G）;holdon;
end
Gc1=zpk（leadlagc（G,20,60,1000,3））%设计超前滞后校正器
Gc2=zpk（leadlagc（G,20,60,1000,1））%设计超前校正器
figure;
step（feedback（G*Gc1,1）,'-',feedback（G*Gc2,1）,':
',1）%绘制闭环响应
结果如下：
Transferfunction:
4s^2+6s+2
-------------------------------------------
s^5+32.1s^4+263.2s^3+426s^2+40s
其中
为超前滞后校正器，
为超前校正器，如图7，图8，图9所示。
如果要知道相位裕度为60，不同剪切频率
对应的单位阶跃响应和伯德图，键入如下命令即可，如图10、11、12、13所示。
图7图8
图9
图10图11
图12图13
gam=60;f1=figure;f2=figure;%打开两个图形窗口
forwc=5:
5:
30
Gc=leadlagc（G,wc,gam,1000,3）;[a,b,c,d]=margin（Gc*G）;
figure（f1）;step（feedback（G*Gc,1）,3）;holdon
figure（f2）;bode（Gc*G）;holdon;
end
5.2基于状态空间模型的设计方法
使用Matlab辅助设计状态空间模型主要是对状态空间模型控制器的设计。
5.2.1状态反馈控制结构框图
如图14所示
图14
将u（t）=v（t）-Kx（t）代入开环系统的状态方程模型，则在状态反馈矩阵K下，系统的闭环状态方程模型可以写成
如果系统（A,B）完全可控，则选择合适的K矩阵，可以将闭环系统矩A-BK的特征值配置到任意地方。
对于开环状态空间模型（A，B，C，D）来说，在状态反馈向量K下，闭环系统的状态方程可以写成（A-BK，B，C-DK，D）。
就极点配置来说，使用以下算法可以实现，算法如下：
Ackermann算法：
K=acker（A,B,p）;
鲁棒极点配置算法：
K=place（A，B，p）;
place（）函数不适用于含有多重期望极点的问题
acker（）函数可以求解配置多重极点的问题
例：
，
键入命令：
A=[0100;00-10;0001;0050];
B=[01;0-1;00;00];
p=[-0.1;-0.2;-0.5+0.2i;-0.5-0.2i];
K=place（A,B,p）
rank（ctrb（A,B））
5.2.2观测器设计及基于观测调节器设计
如图15所示
图15
原系统的（A，C）为完全可观测，则状态观测器数学模型的状态空间表示为：
方程解析解为:
稳定，
输入.m函数，如下
function[xh,x,t]=simobsv（G,L）
[y,t,x]=step（G）;G=ss（G）;A=G.a;B=G.b;C=G.c;D=G.d;
[y1,xh1]=step（（A-L*C）,（B-L*D）,C,D,1,t）;
[y2,xh2]=lsim（（A-L*C）,L,C,D,y,t）;xh=xh1+xh2;
G为对象的状态方程的状态对象模型，L为观测器向量。
例：
。
键入命令：
P=[-10;-10;-10;-10];L=acker（A',C',P）';L'%设计新观测器
[xh,x,t]=simobsv（ss（A,B,C,D）,L）;
plot（t,x,t,xh,':
'）;set（gca,'XLim',[0,30],'YLim',[-0.5,4]）
figure
结果如图16所示：
图16
如果极点位于10
键入命令：
P=[-10;-10;-10;-10];L=acker（A',C',P）';L'%设计新观测器
[xh,x,t]=simobsv（ss（A,B,C,D）,L）;
plot（t,x,t,xh,':
'）;set（gca,'XLim',[0,30],'YLim',[-0.5,4]）
可得结果：
ans=
-421.1634260.725533.2946-20.8091
带有观测器的状态反馈控制结构图
如图17所示：
图17
求取控制器和反馈环节的Matlab函数
function[Gc,H]=obsvsf（G,K,L）
H=ss（G.a-L*G.c,L,K,0）;
Gc=ss（G.a-G.b*K-L*G.c+L*G.d*K,G.b,-K,1）;
如果参考输入信号r（t）=0,控制结构Gc（s）化简为
c=reg（G,K,L）
则有Gc=reg（G,K,L）
例：
，选择加权矩阵Q=diag（0.01,0.01,2,3）,R=1,设计LQ最有控制器。
键入命令：
A=[0,2,0,0;0,-0.1,8,0;0,0,-10,16;0,0,0,-20];
B=[0;0;0;0.3953];C=[0.09882,0.1976,0,0];D=0;
Q=diag（[0.01,0.01,2,3]）;R=1;%输入加权矩阵
K=lqr（A,B,Q,R）,step（ss（A-B*K,B,C,D））%设计{\rmLQ}最优控制器
结果如下：
K=
0.10000.94290.76630.6387
如图18所示：
设观测器极点均位于-5，用极点配置发设计极点观测器：
键入命令：
P=[-5;-5;-5;-5];G=ss（A,B,C,D）;L=acker（A',C',P）';%设计观测器
[Gc,H]=obsvsf（G,K,L）;%设计控制器
step（ss（A-B*K,B,C,D）,feedback（G*Gc,H））%比较基于观测器的控制器与状态反馈
结果如图19所示：
图18
图19
基于观测器下控制器闭环系统的最小实现模型：
键入命令：
zpk（minreal（feedback（G*Gc,H）））%最小实现模型
zpk（minreal（ss（A-B*K,B,C,D）））%这个结果和上式忽略两个一阶零点几乎一致
结果如下：
Zero/pole/gain:
9.9982（s+1）
----------------------------------------------
（s+20.01）（s+10.01）（s^2+0.3341s+0.05052）
Zero/pole/gain:
9.9982（s+1）
----------------------------------------------
（s+20.01）（s+10.01）（s^2+0.3341s+0.05052）
参考文献：
[1]薛定宇控制系统计算机辅助设计——MATLAB语言与应用（第二版）.北京：
清华大学出版社
[2]王万良自动控制原理.北京：
高等教育出版社
[3]刘豹，唐万生现代控制理论（第三版）.北京：
机械工业出版社
升级会员 