自动控制原理MATLAB第四五作业.docx
《自动控制原理MATLAB第四五作业.docx》由会员分享，可在线阅读，更多相关《自动控制原理MATLAB第四五作业.docx（16页珍藏版）》请在冰豆网上搜索。
自动控制原理MATLAB第四五作业
基于MATLAB的状态反馈极点配置函数设计
本文用MATLAB设计了一个用于控制系统状态反馈极点配置的placen（）函数。
该函数调用格式为：
k=placen（A,B,P）。
式中：
输入参量A为系统的状态矩阵；B为输入矩阵；P为指定的系统闭环极点。
返回参量k为反馈增益矩阵即状态反馈矩阵。
其功能和MATLAB自带的place（）函数基本相同。
极点配置算法如下：
（1）判断系统能控性
（2）计算A的特征多项式
（3）计算由
所决定的希望的特征多项式
（4）计算
（5）计算变换阵
（6）计算反馈增益向量
程序清单如下：
functionK=placen（A,B,P）
[nx,na]=size（A）;%判断A矩阵是否合法
ifna~=nx,
error（'A必须是方阵！
'）
elseifnx==0,
K=[];
return
end
[n,m]=size（B）;%判断B矩阵是否合法
ifnx~=n,
error（'B和A的行数必须相同！
'）
end
S=ctrb（A,B）;%S为能控性矩阵
ifrank（S）~=n,%判断系统是否可控
error（'系统不可控！
'）
end
iflength（P）~=nx,%判P（系统指定闭环极点）向量是否合法
error（'P的数量必须和A的行数相同！
'）
end
j=0;
fori=1:
n,
if（real（P（i））>0）,
j=1;
end
end
ifj~=0,
warning（'P的实部应小于0！
'）
elseif~isequal（sort（P（imag（P）>0））,sort（conj（P（imag（P）<0））））
warning（'P中复数应共轭！
'）
end
symss;
I=eye（n）;%创建单位阵
Det=det（s*I-A）;%计算A的特征多项式
aq=fliplr（sym2poly（Det））;%取A的特征多项式系数，并反转
ah=fliplr（poly（P））;%取由P构成的特征多项式的系数，并反转
k=ah-aq;
k=k（1:
n）;%计算
a=aq（2:
n）;%以下循环构建矩阵
A1=[];
fori=1:
n,
forj=1:
n,
if（j<=n-i）,
ifi>j,
A1（i,j）=a（i）;
else
A1（i,j）=a（j）;
end
elseif（j==n-i+1）,
A1（i,j）=1;
else
A1（i,j）=0;
end
end
end
Pc1=S*A1;%构建Pc1
K=k*inv（Pc1）;%
现举例运行如下：
例一：
（课本P224页习题6.15）：
输入：
A=[-21;0-1];
B=[0;1];
P=[-3-3];
K=placen（A,B,P）
输出：
K=13
二维系统结果正确！
例二：
（课本P204页例题6.14）：
输入：
A=[0100;0010;0-6-32];
B=[0;0;450];
P=[-7+7i-7-7i-100];
K=placen（A,B,P）
输出：
K=0.21780.31960.1822
三维系统结果正确！
实际设计状态反馈向量K时，不必通过繁琐的线性变换途径。
在检查系统完全能控之后，只需令A-BK的特征多项式等于要求的特征多项式，由两多项式同次幂系数相等，便可解出K中的各元素，但需要已知A、B和P。
以例二数据为例，算法如下：
输入：
A=[0100;0010;0-6-32];
B=[0;0;450];
P=[-7+7i-7-7i-100];
symssk1k2k3;
I=eye（3）;
k=[k1k2k3];
as=det（s*I-A+B*k）
ao=poly（P）
输出：
as=s^3+32*s^2+450*s^2*k3+60*s+4500*s*k2+45000*k1
ao=111414989800
由对应系数相等可得：
[k1k2k3]=solve（'45000*k1=9800','4500*k2+60=1498','450*k3+32=114'）;
K=[k1k2k3]
输出：
K=[49/225,719/2250,41/225]
结果正确!
本人目前还不能够利用此算法设计出一个类似于place（）函数功能的函数，还有待进一步提高。
P2256.23
设系统状态空间描述为
，
试利用MATLAB的CAD程序设计状态反馈向量使闭环极点在
，-2，并画出系统波特图、根轨迹图，求闭环系统零、极点与静态增益，并求闭环系统的单位阶跃响应。
分析：
此题意在对比极点配置前后系统外部性能的变化。
系统：
引入状态反馈后的闭环系统为：
程序清单如下：
%本程序能对任意能控系统进行极点配置，并画出配置前后系统波特图、根轨迹图，求出配置前后闭环系统零、极点与静态增益，并绘出配置前后闭环系统的单位阶跃响应曲线。
A=input（'请输入系统矩阵A='）;
B=input（'请输入输入矩阵B='）;
C=input（'请输入输出矩阵C='）;
D=input（'请输入前向反馈矩阵D='）;
[num,den]=ss2tf（A,B,C,D,1）;
figure
（1）;
bode（num,den）;grid;%绘制波特图
title（'极点配置前波特图'）;
figure
（2）;
rlocus（num,den）;%绘制根轨迹
title（'极点配置前根轨迹'）;grid;%求零极点和静态增益
[ZPK]=tf2zp（num,den）
figure（3）;
step（num,den）;grid;%绘单位阶跃响应曲线
title（'极点配置前单位阶跃响应'）;
P=input（'请输入期望的系统闭环极点向量P='）;
K=place（A,B,P）;%place函数求出状态反馈向量，若A,B,P
disp（'状态反馈向量K='）;disp（K）;%不合法或系统不可控，place（）函数将给出提示
A=A-B*K;%求极点配置后的系统矩阵
[num,den]=ss2tf（A,B,C,D,1）;
figure（4）;
bode（num,den）;grid;%绘制波特图
title（'极点配置后波特图'）;
figure（5）;
rlocus（num,den）;%绘制根轨迹
title（'极点配置后根轨迹'）;grid;
[ZPK]=tf2zp（num,den）%求零极点和静态增益
figure（6）;
step（num,den）;%绘单位阶跃响应曲线
title（'极点配置后单位阶跃响应'）;grid;
运行过程及结果：
输入：
请输入系统矩阵A=[-0.40-0.01;100;-1.49.8-0.02]
请输入输入矩阵B=[6.3;0;9.8]
请输入输出矩阵C=[001]
请输入前向反馈矩阵D=0
请输入期望的系统闭环极点向量P=[-1+i-1-i-2]
输出：
状态反馈向量K=
0.47061.00000.0627
%配置前
Z=0.2500+2.4975i
0.2500-2.4975i
P=-0.6565
0.1183+0.3678i
0.1183-0.3678i
K=9.8000
%配置后
Z=0.2500+2.4975i
0.2500-2.4975i
P=-2.0000
-1.0000+1.0000i
-1.0000-1.0000i
K=9.8000
%由上可知：
状态反馈可以改变系统特征值，也就是改变了系统传递函数的极点，但却不能改变系统传递函数的零点，也没有改变静态增益。
由上可知：
若系统完全能控，状态反馈极点配置可以任意配置系统极点，从而改变系统的动态特性。
P2256.24
利用MATLAB的CAD程序对上题设计带全阶观测器状态反馈控制系统（设观测器极点为
，-4）；求反馈向量L，求观测器的数学模型，并画出系统的波特图，根轨迹图，求系统增益欲量和相位欲量，求闭环系统零、极点与静态增益，球闭环系统的单位阶跃响应。
设状态初态为
，观测器初态为
，比较
与
的动态误差情况。
分析：
系统：
带观测器的状态反馈控制系统为：
误差方程为：
程序清单如下：
%本程序能对任意能控能观系统进行带全阶观测器极点配置，并画出配置后系统波特图、根轨迹图，求出配后闭环系统零、极点与静态增益，绘出配置后闭环系统的单位阶跃响应曲线，并求出估计误差衰减方程。
A=input（'请输入系统矩阵A='）;
B=input（'请输入输入矩阵B='）;
C=input（'请输入输出矩阵C='）;
D=input（'请输入前向反馈矩阵D='）;
CP=input（'请输入期望的系统闭环极点向量CP='）;
OP=input（'请输入观测器极点向量OP='）;
K=place（A,B,CP）;%求出反馈向量
L=acker（A',C',OP）;
L=L';
disp（'反馈向量L='）;disp（L）;
AHC=A-L*C;%求观测器数学模型
disp（'状态观测器矩阵AHC='）;disp（AHC）;
na=rank（A）;
A=[A-B*KB*K;zeros（na）A-L*C];%
[nb1,nb2]=size（B）;
B=[B;zeros（na,nb2）];%
C=[Czeros（nb2,na）];%
[num,den]=ss2tf（A,B,C,D,1）;
figure
（1）;
bode（num,den）;grid;%绘制波特图
title（'波特图'）;
figure
（2）;
rlocus（num,den）;grid;%绘制根轨迹
title（'根轨迹'）;
[ZPK]=tf2zp（num,den）%求闭环零极点
figure（3）;
step（num,den）;%求闭环单位阶跃响应
title（'单位阶跃响应'）;grid;
c1=input（'请输入状态初态：
'）;%输入状态初态并判断是否输入正确
[nc11,nc12]=size（c1）;
ifnc11~=nb1,
error（'输入行数不对'）
elseifnc12~=nb2,
error（'输入列数不对'）
end
c2=input（'请输入观测器初态：
'）;%输入观测器初态并判断是否输入正确
[nc21,nc22]=size（c2）;
ifnc21~=nb1,
error（'输入行数不对'）
elseifnc22~=nb2,
error（'输入列数不对'）
end
c=c1-c2;
I=eye（na）;
symss;
x=ilaplace（inv（s*I-AHC））*c%求出估计误差
运行过程及结果：
输入：
请输入系统矩阵A=[-0.40-0.01;100;-1.49.8-0.02]
请输入输入矩阵B=[6.3;0;9.8]
请输入输出矩阵C=[001]
请输入前向反馈矩阵D=0
请输入期望的系统闭环极点向量CP=[-1+i-1-i-2]
请输入观测器极点向量OP=[-3+3i-3-3i-4]
输出：
反馈向量L=5.4664
4.6762
9.5800
状态观测器矩阵AHC=
-0.40000-5.4764
1.00000-4.6762
-1.40009.8000-9.6000
Z=-4.0000
-3.0000+3.0000i
-3.0000-3.0000i
0.2500+2.4975i
0.2500-2.4975i
P=-3.0000+3.0000i
-3.0000-3.0000i
-4.0000
-2.0000
-1.0000+1.0000i
-1.0000-1.0000i
K=9.8000
%将上面结果与6.23结果相比较，可知：
全阶观测器增加了N个系统极点，但同时也增加了N个相同的零点，所以复合系统的零极点并没有改变，仍然是状态反馈极点配置的极点。
同时静态增益也没有改变。
绘出系统波特图后，即可得到：
增益裕量：
6.8
相位裕量：
143+180=323
5.96+180=185.96
-60.4+180=119.6
将波特图、根轨迹图和单位阶跃响应曲线同6.23相比较，基本相同（仅此处多了3个对消的零极点）。
可知：
加入观测器不影响闭环系统外部特性
输入：
请输入状态初态：
[1;2;3]
请输入观测器初态：
[-1;-2;-3]
输出：
%估计误差衰减方程
x=182556/32375*exp（-3*t）*cos（3*t）-337752/32375*exp（-3*t）*sin（3*t）-117806/32375*exp（-4*t）
266814/45325*exp（-3*t）*cos（3*t）-240888/45325*exp（-3*t）*sin（3*t）-85514/45325*exp（-4*t）
-299/125*exp（-4*t）+1049/125*exp（-3*t）*cos（3*t）-233/125*exp（-3*t）*sin（3*t）
由上已经得到
，可通过如下命令绘制误差估计衰减曲线：
>>t=0:
0.1:
5;
>>figure（4）;
%绘制
>>plot（t,182556/32375*exp（-3*t）.*cos（3*t）-337752/32375*exp（-3*t）.*sin（3*t）-117806/32375*exp（-4*t）,'r'）
>>holdon
%绘制
>>plot（t,266814/45325*exp（-3*t）.*cos（3*t）-240888/45325*exp（-3*t）.*sin（3*t）-85514/45325*exp（-4*t）,'b'）
>>holdon
%绘制
>>plot（t,-299/125*exp（-4*t）+1049/125*exp（-3*t）.*cos（3*t）-233/125*exp（-3*t）.*sin（3*t）,'g'）
>>grid;
>>title（'估计误差衰减曲线'）;
、
、
图形分别如下所示：
由以上三估计误差衰减曲线可以看出：
误差大约在1.5秒后消除。
升级会员 