自动控制系统仿真

曹晓平　刘渊

【摘要】 衡量控制系统调节品质优劣的性能指标可以理解为对每一类自动控制系统被控量变化全过程提出的共同基本要求，即稳定性、快速性和准确性。其中，稳定性是首要的，那么如何准确判断控制系统是否稳定就是十分重要的问题，通过设计程序算法和控制系统仿真，能够很好地解决这个问题。

【关键词】 系统稳定性 劳斯判据 仿真 算法 性能

引言

本论文通过阐述如何判断自动控制系统品质优劣的问题入手，研究了判断稳定性如今最前沿的方法，利用软件编程的策略，这样可以通过仿真的方法准确判断任一闭环控制系统是否稳定，是将理论结合实际工作的典范。

一、设计内容及意义

系统特征方程如下，试用劳思判据和霍尔维茨判据判断系统稳定性.s^5+3s^4+10s^3+24s^ 2+32s+48=0

自动控制系统的共同基本要求可以归结为稳、准、快三点。稳定是控制系统的重要性能，也是系统能够正常运行的首要条件。一个稳定的控制系统，其被控量偏离期望值的初始偏差经过一个过渡过程时间应逐渐减小并趋于零。而不稳定的系统无法实现预定的控制任务。因此，稳定性判断对自动控制系统是非常紧要的，而自动控制理论的基本任务之一就是如何分析系统的稳定性并提出保证系统稳定的措施。根据题目要求，运用控制原理相关知识，分析所给系统的稳定性，并结合控制系统仿真和matlab知识，编写出计算机实现程序，判断稳定性。

二、稳定性分析

平衡状态稳定性概念是由俄国学者李雅普诺夫1892年首先提出，根据该稳定性理论，线性系统稳定性定义为：

线性控制系统在初始扰动影响下，其动态过程随时间推移逐渐衰减并趋于零（或原平衡工作点），系统最后可以达到平衡状态，则称该系统渐进稳定，简称稳定；若在初始扰动影响下，其动态过程随时间推移而发散，系统被控量失控，则称系统不稳定；若在初始扰动影响下，其动态过程随时间的推移虽不能回到原平衡点，但能够保持在原工作点附近的某一有限区域内运动，表现为等幅振荡形式则称系统临界稳定。

线性系统的稳定性取决于系统自身的结构和参数，与外界因素条件无关。线性系统稳定的充分必要条件是：闭环系统特征方程的所有根全部都具有负实部；或者说，闭环传递函数的极点全部都严格位于左半s平面。

根据稳定的充分必要条件判别系统的稳定性，必然要求准确求出系统的全部特征根。对于高阶系统来说，求根的工作量十分庞大，那么就希望使用一种间接判断系统特征根是否全部严格位于s左半平面的代替简便方法。

三、算法及对象选定

3.1 劳斯判据

将各项系数，按下面的格式排成劳斯表

s^n a0 a2 a4 a6 …

s^（n-1） a1 a3 a5 a7 …

s^（n-2） b1 b2 b3 b4

s^（n-3） c1 c2 c3

.

.

s^2 d1 d2 d3

s^1 e1 e2

s^0 f1

表中：

这样可求得n+1行系数

劳斯稳定判据是根据所列劳斯表第一列系数符号的变化，去判别特征方程式根在S平面上的具体分布，过程如下：

①如果劳斯表中第一列的系数均为正值，则其特征方程式的根都在S的左半平面，相应的系统是稳定的。

②如果劳斯表中第一列系数的符号有变化，其变化的次数等于该特征方程式的根在S的右半平面上的个数，也代表方程的正实根个数，相应的系统为不稳定。

③如果劳斯表中某一列出现全零行时，需要用上一行的系数构造一个辅助方程F（s）=0，并将辅助方程对变量s求导，用所得导数方程的系数取代全零行的元，便可按劳斯稳定判据的要求继续运算下去，直到得出完整的劳斯计算表。

④如果劳斯表中某一行的第一列為零，而其余各项不为零，或不全为零时，以一个很小的正数来代替为零的这项，据此算出其余的各项，完成劳斯表的排列。

经过手算，对于特征方程 s^5+3s^4+10s^3+24s^2+32s+48=0

列出劳斯表： s^5 1 10 32

s^4 3 24 48

s^3 2 16

s^2 0 48

s^2 e（-10） 48 辅助方程

s^1 16

s^0 48

根据劳斯判据，此系统是临界稳定的。

3.2霍尔维茨判据

系统特征方程s^5+3s^4+10s^3+24s^2+32s+48=0

系统稳定的充要条件：

特征方程的全部系数都为正，且主行列式及对角线上的子行列式都大于零。

3.3 用求特征方程根的方法，判断实部是否小于零来判断稳定性，存在实部大于零的根，则系统不稳定。

3.4也可构成单位负反馈系统，绘制开环幅相曲线，利用奈奎斯特判据来判断系统的稳定性。

四、仿真算法选取与性能分析

根据题目要求，利用劳斯判据和霍而维茨判据判断所给系统的稳定性，这两种算法直接并且实现较为简单。

构造系统的routh表

格式： [rtab，info]=routh（den）

说明：其中den是系统的分母多项式向量，rtab是构造的routh表矩阵，info为字符串型变量，返回有关信息。系统不稳定极点的数目等于所产生的routh表中第一列元素的符号变化次数。

function [rtab，info]=routh（den）

info=[]；

vec1=den（1：2：length（den））； nrT=length（vec1）；

vec2=den（2：2：length（den））；

rtab=[vec1； vec2， zeros（1，nrT-length（vec2））]；

for k=1：length（den）-2

alpha（k）=vec1（1）/vec2（1）；

if mod（length（den），2）==0，n=length（vec1）-1；

else n=length（vec2）；

end

for i=1：n

a3（i）=rtab（k，i+1）-alpha（k）*rtab（k+1，i+1）；

end

if sum（abs（a3））==0

a3=polyder（vec2）；

info=[info，'All elements in row '，...

int2str（k+2） ' are zeros；']；

elseif abs（a3（1））

a3（1）=1e-6；

info=[info，'Replaced first element；']；

end

构造Hurwitz矩阵。

格式： [H，Hz_det]=hurwitz（den）

说明：H为构造的Hurwitz矩阵，Hz_det为各阶主子式的行列式值，den为系统的分母多项式D（S）。

D（S）=a0Sn+ a1Sn-1 + a2Sn-2+ ……+ an-1S1 + an

H 矩阵

function [H，Hz_det]=hurwitz（den）

n=length（den）-1；

for i=1：n

i1=floor（i/2）；

if i==i1*2

hsub1=den（1：2：n+1）； i1=i1-1；

else hsub1=den（2：2：n+1）；

end

l1=length（hsub1）；

H（i，：）=[zeros（1，i1），hsub1，zeros（1，n-i1-l1）]；

end

[nr，nc]=size（H）；

for i=1：nr

Hz_det（i，1）=det（（H（1：i，1：i）））；

end

經过分析计算，仿真编程序实现，知道这个系统是临界稳定的，劳斯表中第四行的第一列为零，而其余各项不为零，或不全为零时，以一个很小的正数来代替为零的这项，据此算出其余的各项，完成劳斯表的排列。而运用霍尔维茨判据，得到行列式为零，也可以判断系统是临界稳定的。

五、结论

经过程序运行与调试，判断出系统是临界稳定的，系统性能并不理想，若在初始扰动影响下，其动态过程随时间的推移虽不能回到原平衡点，但可以保持在原工作点附近的某一有限区域内运动，则称系统临界稳定，在外加一定扰动下，会偏离平衡状态而变为不稳定系统，对控制性能带来不利影响，是我们应该避免的状态，在实际控制系统中应该使系统阶跃相应为衰减振荡的，才能满足一般自动系统的最基本要求，稳定性。

我们设计的这个程序针对不同的情况，可以判断出系统是稳定的，不稳定或者是临界稳定的，适应性很强，而不是简单的只是判断是稳定还是不稳定，这也是它的优越性所在。