chenjob22ding.m

用MATLAB实现用C代码可控性和客观性的判断和转移矩阵的求解等各种功能。

A=input('please enter the matrix A= ');
B=input('please enter the matrix B= ');
C=input('please enter the matrix C= ');
D=input('please enter the matrix D= ');
U=input('please enter the  u= ');
T=input('please enter sampling time T= ');
Th=input('please enter the simulation time Th= ');
x=input('please enter the initial value of x0= ');
y=input('please enter the initial value of y0= ');
E=input('please enter the set error E= ');
[m,n]=size(A);
[p,q]=size(B);
[s,r]=size(C);
M= expm(A*T);
disp('the state transition matrix ')
M
a=0;
z=zeros(n);
for i=1:1:n
    for j=1:1:n
a=a+abs(A(i,j)); %a为A的范数
end
end
for i=0:1:1
z=z+(A*T)^i*T/prod(1:(i+1));
end
m1=det(z);
L=1;
while(1)
M1=zeros(m,n);
 for i=0:1:L
 M1=M1+((A*T)^i*T/(prod(1:i)*(i+1)));
end
if (m1>det(M1))   %求最小值
    m1=det(M1);
end
E=1e-7;   %确定误差范围
e=a*T/(L+2);%计算e来确定是否符合精度,当小于1时前提成立
F=(a*T)^(L+1)/prod(1:(L+1))*(1/(1-e));
if ((e<1)&&(F<=E*m1))
M1=M1*B;
M1
break
else 
    L=L+1;
end
end
Uc=B;
Vo=C;
k=rank(A); %求系统矩阵的维数
disp('the dimension of system is')
k
if(k==n)%判断系统是否满秩
    disp('the matrix of system is full rank')
end
for i=1:n-1 
    U1=A^i*B;
    V1=C*A^i;
    Uc=[Uc, U1];%可控性判断矩阵
    Vo=[Vo;V1];%可观性判断矩阵
end
if(n==rank(Uc))
    if(p==rank(Vo))
        disp('the system is controllable and observable')
    else
        disp('the system is controllable and unobservable')
    end
else if(p==rank(Vo))
        disp('the system is uncontrollable and observable')
    else
        disp('the system is uncontrollable and unobservable')
    end
end
X=zeros(m,(Th/T+1));
Y=zeros(s,(Th/T+1));
X(:,1)=x;
Y(:,1)=y;
t=0:T:Th;
for i=2:Th/T
    for j=1:m
        X(j,i)=M(j,:)*X(:,i-1)+M1(j,:)*U;   %第j个状态变量X(j）在iT时刻的值        
    end
    for j=1:s
        Y(j,i)=C(j,:)*X(:,i)+D(j,:)*U;       %第j个输出在iT时刻的值
    end
end
for i=1:m           %前三行窗口画出状态变量X（t）
    subplot(6,2,i)
    plot(t,X(i,:));
    grid on;    
end
for j=1:s
    subplot(6,2,6+j)   %第后三行窗口画出Y(t)
    plot(t,Y(j,:));
    grid on;
end