自动控制原理MATLAB第四五作业.docx

1、自动控制原理MATLAB第四五作业基于MATLAB 的状态反馈极点配置函数设计 本文用MATLAB设计了一个用于控制系统状态反馈极点配置的placen()函数。该函数调用格式为：k=placen(A,B,P)。式中：输入参量A 为系统的状态矩阵；B为输入矩阵；P为指定的系统闭环极点。返回参量k 为反馈增益矩阵即状态反馈矩阵。其功能和MATLAB自带的place()函数基本相同。极点配置算法如下：（1）判断系统能控性（2）计算A的特征多项式（3）计算由所决定的希望的特征多项式（4）计算 （5）计算变换阵（6）计算反馈增益向量 程序清单如下：function K=placen(A,B,P)nx,n

2、a=size(A); %判断A矩阵是否合法if na=nx, error(A 必须是方阵！)elseif nx=0, K=; returnendn,m=size(B); %判断B 矩阵是否合法if nx=n, error(B和A的行数必须相同！)endS=ctrb(A,B); %S为能控性矩阵if rank(S)=n, %判断系统是否可控 error(系统不可控！)endif length(P)=nx, %判P（系统指定闭环极点）向量是否合法 error(P的数量必须和A的行数相同！)endj=0;for i=1:n, if(real(P(i)0), j=1; endendif j=0, wa

3、rning(P的实部应小于0！)elseif isequal(sort(P(imag(P)0),sort(conj(P(imag(P)0) warning(P中复数应共轭！)endsyms s;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=;for i=1:n, for j=1:n, if(jj, A1(i,j)=a

4、(i); else A1(i,j)=a(j); end elseif(j=n-i+1), A1(i,j)=1; else A1(i,j)=0; end endendPc1=S*A1; %构建Pc1K=k*inv(Pc1); %现举例运行如下：例一：（课本P224页习题6.15）：输入： A=-2 1;0 -1;B=0;1;P=-3 -3;K=placen(A,B,P)输出： K =1 3二维系统结果正确！例二：（课本P204页例题6.14）：输入： A=0 10 0;0 0 10;0 -6 -32;B=0;0;450;P=-7+7i -7-7i -100;K=placen(A,B,P)输出：

5、K =0.2178 0.3196 0.1822三维系统结果正确！实际设计状态反馈向量K时，不必通过繁琐的线性变换途径。在检查系统完全能控之后，只需令A-BK的特征多项式等于要求的特征多项式，由两多项式同次幂系数相等，便可解出K中的各元素，但需要已知A、B和P。以例二数据为例，算法如下：输入： A=0 10 0;0 0 10;0 -6 -32;B=0;0;450;P=-7+7i -7-7i -100;syms s k1 k2 k3;I=eye(3);k=k1 k2 k3;as=det(s*I-A+B*k)ao=poly(P)输出： as =s3+32*s2+450*s2*k3+60*s+4500

6、*s*k2+45000*k1ao =1 114 1498 9800由对应系数相等可得：k1 k2 k3=solve(45000*k1=9800,4500*k2+60=1498,450*k3+32=114);K=k1 k2 k3输出：K = 49/225, 719/2250, 41/225结果正确!本人目前还不能够利用此算法设计出一个类似于place()函数功能的函数，还有待进一步提高。P225 6.23设系统状态空间描述为， 试利用MATLAB的CAD程序设计状态反馈向量使闭环极点在，-2，并画出系统波特图、根轨迹图，求闭环系统零、极点与静态增益，并求闭环系统的单位阶跃响应。分析：此题意在对比

7、极点配置前后系统外部性能的变化。系统： 引入状态反馈后的闭环系统为： 程序清单如下：%本程序能对任意能控系统进行极点配置，并画出配置前后系统波特图、根轨迹图，求出配置前后闭环系统零、极点与静态增益，并绘出配置前后闭环系统的单位阶跃响应曲线。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);rlo

8、cus(num,den); %绘制根轨迹title(极点配置前根轨迹);grid; %求零极点和静态增益Z P K=tf2zp(num,den)figure(3);step(num,den);grid; %绘单位阶跃响应曲线title(极点配置前单位阶跃响应);P=input(请输入期望的系统闭环极点向量 P =);K=place(A,B,P); %place函数求出状态反馈向量，若A,B,Pdisp(状态反馈向量 K =);disp(K); %不合法或系统不可控，place()函数将给出提示A=A-B*K; %求极点配置后的系统矩阵num,den=ss2tf(A,B,C,D,1);figur

9、e(4);bode(num,den);grid; %绘制波特图title(极点配置后波特图);figure(5);rlocus(num,den); %绘制根轨迹title(极点配置后根轨迹);grid;Z P K=tf2zp(num,den) %求零极点和静态增益figure(6);step(num,den); %绘单位阶跃响应曲线title(极点配置后单位阶跃响应);grid;运行过程及结果：输入：请输入系统矩阵 A = -0.4 0 -0.01;1 0 0;-1.4 9.8 -0.02请输入输入矩阵 B = 6.3;0;9.8请输入输出矩阵 C = 0 0 1请输入前向反馈矩阵 D = 0

10、请输入期望的系统闭环极点向量 P =-1+i -1-i -2输出：状态反馈向量 K = 0.4706 1.0000 0.0627%配置前Z =0.2500 + 2.4975i 0.2500 - 2.4975iP = -0.6565 0.1183 + 0.3678i 0.1183 - 0.3678iK =9.8000%配置后Z =0.2500 + 2.4975i 0.2500 - 2.4975iP =-2.0000 -1.0000 + 1.0000i -1.0000 - 1.0000iK =9.8000%由上可知：状态反馈可以改变系统特征值，也就是改变了系统传递函数的极点，但却不能改变系统传递函

11、数的零点，也没有改变静态增益。由上可知：若系统完全能控，状态反馈极点配置可以任意配置系统极点，从而改变系统的动态特性。P225 6.24利用MATLAB的CAD程序对上题设计带全阶观测器状态反馈控制系统（设观测器极点为，-4）；求反馈向量L，求观测器的数学模型，并画出系统的波特图，根轨迹图，求系统增益欲量和相位欲量，求闭环系统零、极点与静态增益，球闭环系统的单位阶跃响应。设状态初态为，观测器初态为，比较与的动态误差情况。分析：系统： 带观测器的状态反馈控制系统为：误差方程为： 程序清单如下：%本程序能对任意能控能观系统进行带全阶观测器极点配置，并画出配置后系统波特图、根轨迹图，求出配后闭环系统

12、零、极点与静态增益，绘出配置后闭环系统的单位阶跃响应曲线，并求出估计误差衰减方程。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(A

13、HC);na=rank(A);A=A-B*K B*K;zeros(na) A-L*C; %nb1,nb2=size(B);B=B;zeros(na,nb2); %C=C zeros(nb2,na); %num,den=ss2tf(A,B,C,D,1);figure(1);bode(num,den);grid; %绘制波特图title(波特图);figure(2);rlocus(num,den);grid; %绘制根轨迹title(根轨迹);Z P K=tf2zp(num,den) %求闭环零极点figure(3);step(num,den); %求闭环单位阶跃响应title(单位阶跃响应);g

14、rid;c1=input(请输入状态初态：); %输入状态初态并判断是否输入正确nc11,nc12=size(c1);if nc11=nb1, error(输入行数不对)elseif nc12=nb2, error(输入列数不对)endc2=input(请输入观测器初态：); %输入观测器初态并判断是否输入正确nc21,nc22=size(c2);if nc21=nb1, error(输入行数不对)elseif nc22=nb2, error(输入列数不对)endc=c1-c2;I=eye(na);syms s;x=ilaplace(inv(s*I-AHC)*c %求出估计误差运行过程及结果：

15、输入：请输入系统矩阵 A = -0.4 0 -0.01;1 0 0;-1.4 9.8 -0.02请输入输入矩阵 B = 6.3;0;9.8请输入输出矩阵 C = 0 0 1请输入前向反馈矩阵 D = 0请输入期望的系统闭环极点向量 CP =-1+i -1-i -2请输入观测器极点向量 OP =-3+3i -3-3i -4输出：反馈向量 L = 5.4664 4.6762 9.5800状态观测器矩阵 AHC = -0.4000 0 -5.4764 1.0000 0 -4.6762 -1.4000 9.8000 -9.6000Z = -4.0000 -3.0000 + 3.0000i -3.000

16、0 - 3.0000i 0.2500 + 2.4975i 0.2500 - 2.4975iP = -3.0000 + 3.0000i -3.0000 - 3.0000i -4.0000 -2.0000 -1.0000 + 1.0000i -1.0000 - 1.0000iK = 9.8000%将上面结果与6.23结果相比较，可知：全阶观测器增加了N个系统极点，但同时也增加了N个相同的零点，所以复合系统的零极点并没有改变，仍然是状态反馈极点配置的极点。同时静态增益也没有改变。 绘出系统波特图后，即可得到：增益裕量：6.8相位裕量：143+180=323 5.96+180=185.96 -60.4

17、+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/

18、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) hold on%绘制 plot(t,266814/45325*exp(-3*t).*cos(3*t)-240888/45325*exp(-3*t).*sin(3*t)-85514/45325*exp(-4*t),b) hold on%绘制 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秒后消除。