共享关于系统基点配置的 matlab 程序

bluewave

clear all
A=[1 0 0;0 0 1;0 5 0]
B=[1;1;0]
%%*****************************************************************%%
%

%计算变换矩阵－－－－－－－这里只考虑到系数矩阵维数为2到4维的情况，如有需要可添加程序。
%其中系统矩阵已经给出，如有需要，修改系统矩阵即可。
%JJ=eye(n)设所需极点为n维单位阵的特征值,需要时重设其极点即可。
%%*****************************************************************%%
%能控性判断
[m,n]=size(A);%m为行数，n为A的列数
[n,p]=size(B);
Q=[B A*B];
E=A*B;
N=n-1;
for i=2:1:N
   E=A*E;
   Q=[Q E]
end
rank(Q)
if rank(Q)==n
y=1;
f='系统为可控系统';
disp(f);
else
y=0;
f='系统为不可控系统，停止计算';
disp(f);
end

%循环性判断
if y==1
f='判断系统循环性';
disp(f);
D=jordan(A);
[m,n]=size(D)
for j=1:1:n-1
    if (D(j,j)~=D(j+1,j+1))%d当其特征值两两互异时跳出，必为循环阵。
x=0;
        break
elseif (D(j,j)==D(j+1,j+1))&D(j,j+1)==0 %当其特征值相等时，判断是否在一个约当块中。如果不在则必为循环阵。

        x=0;
        break
    else x=1;
    end
end
if x==0
    fprintf('A为循环矩阵')
A=A
K1=zeros(p,n)

else x==1
    fprintf('A为不循环矩阵')

K1=ones(p,n)%根据定理，如果(A,B)能控，对于几乎任意的p*n维常阵，都可使其为循环矩阵，故此处选取此常阵。
     

end
A=A-B*K1
fprintf('A为新构造循环矩阵')
end


%构造能控系统

    e=ones(p,1);%根据定理，几乎任意的p维实数列向量都有(A,Be)为能控，故选取此矩阵。

    b=B*e;
Q1=[b A*b];%判断新构造系统能控性
E1=A*b;
N1=n-1;
for i=2:1:N1
   E1=A*E1;
   Q1=[Q1 E1]
end
rank(Q1)
if rank(Q1)==n
y=1;
f='新构造系统为可控系统';
disp(f);
else
y=0;
f='新构造系统为不可控系统';
disp(f);
end

%计算反馈增益矩阵K
P=poly(A);

[m,n]=size(A);%m为行数，n为A的列数
[n,p]=size(B);
F=[A*B B];
T=A*B;
N=n-1;
for i=2:1:N
   T=A*T;
   F=[T F];

end

%计算变换矩阵－－－－－－－这里只考虑到系数矩阵维数为2到4的情况，如有需要可添加程序。
if n==2
W=[1 0;P(2) 1];
elseif n==3
W=[1 0 0;P(2) 1 0;P(3) P(2) 1];
elseif n==4
W=[1 0 0 0;P(2) 1 0 0;P(3) P(2) 1 0;P(4) P(3) P(2) 1];
elseif n>=5
f='原程序需要增加W的维数';
disp(f);
end
L=F*W %L为变换矩阵
f='L为变换矩阵';
disp(f);
JJ=eye(n);%设所需极点为n维单位阵的特征值
X=poly(JJ);
Z=X(n+1)-P(n+1);
N=n;
for i=n:-1:2
  M=X(i)-P(i);
   Z=[Z M];
end
K=Z*(inv(L))
f='K为状态反馈增益';
disp(f);

[ 本帖最后由 eight 于 2008-2-18 18:30 编辑 ]