基于系统的可控性进行配置极点(系统仿真)

课程设计报告

学生姓名:王强学号：0907240537 学院: 自动化工程学院

班级: 自动095

题目: 专业方向课程设计

基于系统的可控性进行配置极点

指导教师：姜文娟顾大可职称: 助教副教授

2013年 1月 9日

目录

1题目背景与意义 (2)

2 设计题目介绍 (2)

3设计步骤 (2)

3.1状态空间分析 (2)

3.2状态反馈的极点配置问题 (3)

3.3极点配置问题的算法实现 (3)

4理论计算 (4)

5仿真程序 (5)

5.1 主程序流程图 (5)

5.2 极点配置计算程序流程图 (6)

5.3 仿真程序设计 (7)

5.3.1 矩阵及参数输入 (7)

5.3.2 输入矩阵及参数的合法性判断 (7)

5.3.3 系统能控性判断 (8)

5.3.4 对系统进行极点配置 (8)

6仿真及结果分析 (8)

6.1 MATLAB运行及操作 (8)

6.2 仿真结果 (8)

6.3 程序通用性验证 (9)

结论 (9)

参考文献 (9)

附录 (10)

程序清单 (10)

1题目背景与意义

20世纪50年代后期，控制理论由经典控制理论向现代控制理论转变，现代控制理论是在引入状态空间概念的基础上发展起来的。与经典控制理论一样，现代控制系统中仍然主要采用反馈控制结构，但不同的是，经典控制理论中主要采用输出反馈，而现代控制系统中主要采用内部状态反馈。状态反馈可以为系统控制提供更多的信息反馈，从而实现更优的控制。闭环系统极点的分布情况决定于系统的稳定性和动态品质，因此，可以根据对系统动态品质的要求，规定闭环系统的极点所应具备的分布情况，把极点的配置作为系统的动态品质指标。这种把极点配置在某位置的过程称为极点配置。在空间状态法中，一般采用反馈系统状态变量或输出变量的办法，来实现系统的极点配置。

在系统的分析和综合中，所涉及的计算主要为矩阵运算和矩阵变换，MATLAB 为此提供了一个强有力的工具。利用MATLAB 语言编制实现此算法的通用程序。

2 设计题目介绍

本课程设计以自动控制理论、现代控制理论、MATLAB 及应用等知识为基础，基于系统的可控性进行配置极点，目的是使学生在现有的控制理论的基础上，学会用MATLAB 语言编写控制系统极点校正的设计与性能分析的程序，通过上机实习加深对课堂所学知识的理解，掌握实现极点校正的设计的方法。 1 在理论上对连续系统推导出可控的判定方法和计算公式 2 生成可控性矩阵

3 判别可控性矩阵的秩数（画出计算机实现算法的框图、编写程序并调试和运行）

4 实现极点配置

5 以下面的系统为例，进行计算

????

?????

?--=34

0011

121A ，???

??

?????=100B ，[]111

-=c 。

1）判定系统是否可任意配置极点；

2）若指定其极点在 -3+i ，-3-i ，-5处，求状态反馈矩阵K 。 6 秩数计算选主元精度为ep=0.001

7 程序应具有一定的通用性，对不同参数能有兼容性。

3设计步骤

3.1状态空间分析

对于控制系统Cx y Bu Ax x =+=?

,，式中A 为n ×n 系数矩阵；B 为n ×r

输入矩阵；C 为m ×n 输出矩阵。选择控制信号为：μ=﹣KX

图1 状态反馈闭环控制原理图

求解上式，得到Cx y Bv x BK A x =+-=?

,)(

可以看出，如果系统状态完全可控，K 选择适当，对于任意的初始状态，当t 趋于无穷时，都可以使x(t)趋于0。

3.2状态反馈的极点配置问题

线性系统的动态性能，系统稳定性，间域分析中的超调量、过渡时间等指标，要取决于系统的极点位置。点配置的一般方法可以通过换算(如根轨迹法)和经验估计，具体地加以确定。闭环极点组配置到所希望的位置上，价于使综合得到的线性系统的动态性能达到期望的要求。以证明，定常系统可通过状态变量反馈来任意配置其全部极点的充要条件是该系统为完全能控的，系统是完全能控的，特征值为Ki ，定的闭环系统的n 个所期望的极点是( K*1 , K*2 , , , K*n ) ，有K*i = Ki( A- BK) i = 1, 2, , , n 。

因此，于一个完全能控的线性系统的极点配置问题，际上转化为求解状态反馈增益矩阵K 。

3.3极点配置问题的算法实现

单输入系统极点配置的计算步骤如下:

(1)检查系统的能控性。如果系统能控，则按以下步骤计算状态反馈矩阵。

(2)计算矩阵A 的特征多项式，即0111]det[a s a s a s A sI n n n ++++=--- 。 (3)计算由希望的极点}{2

1

n λλλ

所决定的状态反馈系统特征多项式

*

1*1

1

*21*

)())(()(a

s a s

a

s s s s s a n n n

n ++++=---=-- λλλ。

(4)计算-

k

][11

*

1

1

*0

0*---

---=n n a a

a a

a a k 。 (5)确定变换矩阵P

???????

?????????=---00

1

001011][132

121

1

1

n n n c a a a a a a B A AB B

p (6)计算状态反馈矩阵K

1

1--

=c P k k

注意，这里的变换矩阵1

1-c P 就是将不是能控标准形的状态方程变成能控标准形的变换矩阵。当系统方程为能控标准形的情况，1

1-c P =I ，这时状态反馈矩阵K=-

k 。

在MATLAB 中用rank 指令易的矩阵的秩，用place()函数易得极点配置，但

在本次设计中将不采用这两条语句。这样就需要我们自己设计算法来实现课题需要。

4理论计算

由题目可知如下条件：

????

?????

?--=34

0011

121

A ，???

??

?????=100B ，[]111

-=c ，P=[-3+i ，-3-i ，-5]。

按步骤进行计算： 1) 系统能控性矩阵：Q c =][1

B A AB

B

n -

=????

?

????

?---93

1100

410

经计算矩阵Q c 的秩为3，所以系统能控。

2) 计算矩阵A 的特征多项式：Δ(s)==-]det[A sI s^3 - 5*s^2 + 5*s – 1 3) 希望极点所决定的状态反馈系统特征多项式为：

=

-∏

=3

1

)(i si s s^3 +11*s^2 + 40*s +50

4) 计算-k ：-

k =[a*0-a 0 a*1-a 1 a*2-a 2]=[51 35 16] 5)

计算得变换矩阵11-c P =????

????

?

?-----13

2011

010 6) 计算状态反馈矩阵K

1

1--

=c P k k =[-67 -134 16]

5仿真程序

5.1 主程序流程图

图2主程序流程图

5.2 极点配置计算程序流程图

图3极点配置计算程序流程图

5.3 仿真程序设计

5.3.1 矩阵及参数输入

在MATLAB中，从键盘输入数据需要用到input指令，对一个系统进行极点配置，需要输入系统矩阵、输入矩阵、输出矩阵、指定极点。所以程序可写成：

A=input('请输入系统矩阵A = ');

B=input('请输入输入矩阵 B = ');

C=input('请输入输出矩阵C = ');

P=input('请输入指定极点P= ');

5.3.2 输入矩阵及参数的合法性判断

对系统进行极点配置，所输入的矩阵及指定极点都是有要求的，其中系

统矩阵A必须是方阵，输入矩阵B和系统矩阵A的行数必须相同，指定闭

环极点P的数量必须和A的行数相同，且P的实部应小于0，以及P中的

复数应共轭。这些判断的程序如下：

[nx,na]=size(A);

if na~=nx,

error('A 必须是方阵！')

elseif nx==0,

return

end

[n,m]=size(B);

if nx~=n,

error('B和A的行数必须相同！')

end

if length(P)~=nx,

error('P的数量必须和A的行数相同！')

end

j=0;

for i=1:n,

if(real(P(i))>0),

j=1;

end

end

if j~=0,

error ('P的实部应小于0！')

else

if ~isequal(sort(P(imag(P)>0)),sort(conj(P(imag(P)<0))))

error ('P中复数应共轭！')

end

end

???

?

??

??????????--0001001011132121

n n a a a a a a 5.3.3 系统能控性判断

判断系统能控性，需要先求系统能控性矩阵，在MATLAB 中，使用指令ctrb(A,B)便可方便得到。之后需要确定系统能控性矩阵是否满秩，如果满秩则能控，否则不可控。由于不使用rank 指令求秩，则需要对系统能控性矩阵进行初等行变换，在MATLAB 中可以使用指令rref ，经两次化成最简形后，求矩阵的行列式，如果为0则说明不是满秩，进而得到系统不可控，反之，则系统可控，可以进行极点配置。以上过程实现的程序如下： S=ctrb(A,B); Q=rref(rref(S)')'; if det(Q)==0.000

disp ('系统不可控！不能任意配置极点！') else

disp ( '系统可控，理论上能任意配置极点。') end

5.3.4 对系统进行极点配置

在MATLAB 中用指令place(A,B,P)可以直接进行极点配置，但在本次课程设计中本着探究的目的需要自行设计算法进行实现极点配置。所以这段的设计中，我参照极点配置的常规步骤，按流程用MATLAB 程序语言实现极点配置的过程。计算A 的特征多项式时，直接使用det(s*I-A)便可以实现；在构建P c 1前，需要用for 循环及条件指令，构造一个矩阵，形如：

最后计算K ，直接按照公式1

1--

=c P k k ，在MATLAB 中将之前求得的部分相乘即得。

6仿真及结果分析

6.1 MATLAB 运行及操作

打开MATLAB 软件，在菜单栏里点击File 找到new 之后弹出并点击Script 或者使用快捷键Ctrl+N 新建一个.m 文件，然后按照之前的设计方案将语句输入，之后在工具栏点击Debug 找到Save file and Run ，就可以运行所编的程序。

6.2 仿真结果

在m 文件运行后，按照提示输入如下：

请输入系统矩阵 A = [1,2,-1;1,1,0;0,-4,3] 请输入输入矩阵 B = [0;0;1]

请输入输出矩阵C = [1,-1,1]

请输入指定极点P= [-3+i,-3-i,-5]

之后按回车键，软件开始运算。结果如下：

系统可控，理论上能任意配置极点。

K =

[ -57.0, -74.0, 16.0]

该结果与前边手工计算的结果一样，初步说明该程序具有可行性。

6.3 程序通用性验证

再次运行该程序，输入与课题给出的一组不同的数据，如下：

请输入系统矩阵A = [-6,-11,-6;1,0,0;0,1,0]

请输入输入矩阵B = [1;0;0]

请输入输出矩阵C = [0,0,10]

请输入指定极点P= [-10,-11,-12]

运行结果如下：

系统可控，理论上能任意配置极点。

K =

[ 27.0, 351.0, 1314.0]

经手工计算及使用place指令求得的极点配置结果验证一致，该程序具有一定的通用性。

结论

通过此次的课程设计，使我对MATLAB有了更加深入的了解，在设计过程中，遇到不会的，通过与老师交流及大量参考资料的查找，强化了对MATLAB 及现代控制理论的认识。并且在过程中，我明白了，实现一个同样的结果可以有多种方法，通过抛开MATLAB中的简便函数使用，更加锻炼我的分析能力及编程思维，在大量语句的编写过程，细心是很重要的。出现错误之后，积极寻找可行的解决办法，拓展了对算法的理解。

参考文献

[1]孙增晰.计算机辅助设计.北京：清华大学出版社，1995

[2]郑大钟.线性系统理论.北京：清华大学出版社，2004

[3]姚全珠.软件技术基础.北京：电子工业出版社，2006

[4]王小丹.基于MA TLAB的系统分析与设计.西安：西安电子科技大学出版社，2007

[5] 王正林.MA TLAB/Simulink与控制系统仿真（第2版）.电子工业出版社，1998.48-54.

[6]赵景波.MA TLAB控制系统仿真与设计.机械工业出版社，2010-11-1.78-90.

附录

程序清单

A=input('请输入系统矩阵A = ');

B=input('请输入输入矩阵 B = ');

C=input('请输入输出矩阵C = ');

P=input('请输入指定极点P= ');

[nx,na]=size(A); %判断A矩阵是否合法

if na~=nx,

error('A 必须是方阵！')

elseif nx==0,

return

end

[n,m]=size(B); %判断B 矩阵是否合法

if nx~=n,

error('B和A的行数必须相同！')

end

S=ctrb(A,B); % S为能控性矩阵]

B n

[1B

AB

A

Q=rref(rref(S)')'; %初等行变换化行最简形

if det(Q)==0.000 %判断系统是否可控

disp ('系统不可控！不能任意配置极点！')

else

disp ( '系统可控，理论上能任意配置极点。')

end

if length(P)~=nx, %判P（系统指定闭环极点）是否合法error('P的数量必须和A的行数相同！')

end

j=0;

for i=1:n,

if(real(P(i))>0),

j=1;

end

end

if j~=0,

error ('P的实部应小于0！')

else

if ~isequal(sort(P(imag(P)>0)),sort(conj(P(imag(P)<0))))

error ('P中复数应共轭！')

end

end

syms s; %定义符号变量

???????

????????

?--00

1001011132121

n n a a a a a a 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);

%计算-

k

a=aq(2:n);

%以下循环构建矩阵:

A1=[]; for i=1:n,

for j=1:n,

if(j<=n-i), if i>j,

A1(i,j)=a(i); else

A1(i,j)=a(j); end else if(j==n-i+1),

A1(i,j)=1; else

A1(i,j)=0; end

end end end

Pc1=S*A1; %构建Pc1 digits(6)

K=vpa(k*inv(Pc1))

%1

1--

=c P k k