routh.m

华东理工大学自动化系《控制系统分析、设计和应用》教材配套程序实例,有关说明： 1. 安装MATLAB软件。安装时

function [row,flag]=routh(den,lamda,strf)
%flag=0,1,2 means system stable,marginally stable and unstable
% lamda used for draw the s=-lamda line 
%den=[1 0 -3 2];
%den=[1 10 35 50 24];
%den=[1 2 1 2];
%den=[1 10 36 60 59 50 24];

% Copyright 2000-2001 ECUST. 
%   $Revision: 1.1 $  $Date: 2001/08/03 16:48:42 $

if nargin==1,lamda=0;end;
if nargin==3,strs=strf;else strs='s';end;
syms s system1;system1=poly2sym(den,'s');den=sym2poly(subs(system1,(s-lamda)));
system=poly2sym(den,'s');no=nargout;
if no==0,disp(['系统特征多项式：',char(system)]);end; dl=length(den);
row1=den(1:2:dl);row2=den(2:2:dl); l1=length(row1);l2=length(row2);
row=[row1;row2,zeros(1,l1-l2)];flag=0;msg=[];
if no==0,disp('Routh阵列如下：');end;
if row(2,1)<eps,row(2,1)=eps;disp('第一元素出现零');end;
for i=1:dl-2,a(i)=row1(1)/row2(1);
    for j=1:(l2-(l2==l1)),b(j)=row(i,j+1)-a(i)*row(i+1,j+1);end;
    if sum(abs(b))==0,
        r2(1:2:2*length(row2)-1)=row2(1:length(row2));
        r2(2:2:2*length(row2)-2)=0; bp=polyder(r2);b=bp(1:2:length(bp));
        msg=['第',num2str(i+2),'行的所有元素为零;'];flag=1;
    elseif abs(b(1))<eps,b(1)=0.00001;msg=[msg,'第一元素出现零;'];
    end;
    row=[row;b,zeros(1,l1-length(b))];row1=row2;row2=b;
end;
[r,u]=size(row);
for i=1:r;bbc=[];
    for j=1:u,
        fi=deblank(num2str(row(i,j)));
        bbc=[bbc,blanks(12-length(fi)),num2str(row(i,j))];end;
    if no==0,disp([strs,'^',int2str(r-i)',' | ',bbc]);end;
end;
disp(['劳斯稳定性判别结果：',msg]);
bf=0;for i=1:r-1,
    if row(i,1)*row(i+1,1)<0,bf=bf+1;end;end;
if bf~=0 ,disp(['系统有',num2str(bf),'个正实部根，系统不稳定']);flag=2;
elseif flag==1,disp('系统临界稳定');else,disp('系统稳定');end;