c5

MTALAB判断可控可观性一．已知系统状态方程为{dXAXBu已知其中A为1×1~5×5YCXDu①编写(T),m(T)数值计算函数②代入实际的A、B、C、D，并判断维数及是否满秩，计算X（t）,Y(t)。当A满秩时判断能控能观性。二．分析1．题意要求（1）对于任意的维数在1×1~5×5的系统编写一个通用的(T),m(T)计算公式；2）维数的计算，满秩的判断；3）能控性和能观性的判断；4）系统响应和输入响应的计算。2．解题思路对于1×1~5×5的系统编写一个通用的(T),m(T)，所以系统矩阵是由键盘任意输入，(T)即计算状态转移矩阵eAT，将eAT用无穷级数来示，然后截取其中的前若干项，若级数在i=L处截断，则可以写成:LiiiiATATATRei!Mi0iL1i!上式中的第一项为级数的近似值，而第二项为余数项。要求：|rmax|E|mmin|，其中E=10-d（d为正整数），其中和对应为R与M的元素，rmax为R中最大元素，mmin为M中绝对值最小的元素。因为R仍是一个无穷项的和。下面来估计。令||R||为矩阵R的范数，根据矩阵的范数定义：nrmax||R|||rij|i,j1||R||||AiTi||A||iTi||A||L1TL1||A||T||A||2T2)||i!(L1)!(1(L3)(LiL1i!iL1L22)||A||L1TL1||A||T||A||2T2)(L(1L2(L2)21)!1令：||A||TL2上式为：||A||L1TL12)||R||(L1)!(1如果1，则：||R||||A||L1TL1(1)(L1)!1nL1TL11所以，如果rmax||A||(|rij|||R||)Emmini,j1(L1)!1则满足要求。计算过程如下：1）选择初始L2）计算矩阵M及mmin3）求（4）判断是否满足精度要求，如满足，则M代替eAT5）如不满足，令L=L+1,重新计算。6）同时Tm(T)(T)Bd,令Tt,0TTTAitiAiti1m(T)(t)BdteAtBdtdtBB00i00i!i0(i1)i!即可以同时计算。系统的维数等于系统矩阵的秩。可以通过使用求秩函数rank()求出。系统的可控性可以通过状态可控性矩阵Uc的秩进行判断，如果Uc的秩等于系统中状态变量X的个数则系统完全可控，否则不完全可控。Uc的计算方法如下：Uc[B,A*BA(n1)*B]系统的可观性可以通过状态可观矩阵Vo的秩进行判断，如果Vo的秩等于系统中输入变量Y的个数，则系统完全可观，否则不完全可观。Vo的计算方法如下：Vo[C,C*AC*A(n1)]T求X(t)建了一个矩阵X，X的行数为系统中状态变量的个数，列数为采样2时间的点数，所以X(i,j)为第i个变量在第j个时刻的值。求Y(t)也用上述类似的方法。最后在一个图中把各个变量随时间的变化和输出随时间的变化画出。因为是系统是1×1~5×5，为了适用系统的任意性，将一个图形窗口分割成6×2部分，前三行的6个窗口用于绘制各个系统变量X（i）随时间的变化图，后三行的6个窗口用于绘制输入变量Y(i)随时间的变化图。三．实现与结果分析1.求状态转移矩阵（1）代码如下A=input('pleaseenterthematrixA=');B=input('pleaseenterthematrixB=');T=input('pleaseenterthesamplingtimeT=');E=input('pleaseentertheseterrorE=');L=1;a=0;[m,n]=size(A);[p,q]=size(B);z1=zeros(n);z=zeros(n);fori=1:1:nforj=1:1:na=a+abs(A(i,j));%B为A的范数endendfori=0:1:1z=z+(A*T)^i/prod(1:i);z1=z1+(A*T)^i*T/prod(1:(i+1));endm1=det(z1);m0=det(z);while(1)M=zeros(m,n);fori=0:1:LM=M+(A*T)^i/prod(1:i);endif(m0>det(M))%求最小值m0=det(M);ende=a*T/(L+2);%计算e来确定是否符合精度,当小于1时前提成立3F=(a*T)^(L+1)/prod(1:(L+1))*(1/(1-e));if((e<1)&&(F<=E*m))Ldisp('Theapproximationofstatetransitionmatrixis')MbreakelseL=L+1;endendL=1;while(1)M1=zeros(m,n);fori=0:1:LM1=M1+((A*T)^i*T/prod(1:i)*(i+1));endif(m1>det(M1))%求最小值m1=det(M1);ende=a*T/(L+2);%计算e来确定是否符合精度,当小于1时前提成立F=(a*T)^(L+1)/prod(1:(L+1))*(1/(1-e));if((e<1)&&(F<=E*m1))LM1=M1*B;M1breakelseL=L+1;endend在命令窗口根据提示输入任意系统矩阵A（例如：A=[1032;0123;9342;0021]），B（例如：B=[10;01;22;32]），仿真步长T（例如：取T=0.1），误差范围（例如：E=1e-7），结果如图2.1中的两个图所示。运行结果如下：pleaseenterthematrixA=[1032;0123;9342;0021];pleaseenterthematrixB=[10;01;22;32]pleaseenterthesamplingtimeT=0.1pleaseentertheseterrorE=1e-7L=184TheapproximationofstatetransitionmatrixisM=1.28320.05930.43710.27760.12481.14680.31250.37101.23530.41181.76180.44520.11380.03790.27451.1414L=23M1=0.59350.36590.44940.45911.26000.94240.57150.41225图2.1求转移矩阵仿真图（2）通过调用MATLAB库函数求状态转移矩阵A=[1032;0123;9342;0021]A=1032012393420021T=0.1T=0.1000expm(A*T)ans=1.28320.05930.43710.27760.12481.14680.31250.37101.23530.41181.76180.445260.11380.03790.27451.1414结果显示如下：图2.2调用库函数求转移矩阵图3）对比通过M近似代替状态转移矩阵的值和通过调用MATLAB库函数求状态转移矩阵比较可以看出在取合适的仿真步长和设定合适的误差范围内，M的值和用MATLAB库函数求的值相同2.求秩X(t),Y(t)，维数及满秩判断、可控性和可观性判断（1）代码如下：A=input('pleaseenterthematrixA=');B=input('pleaseenterthematrixB=');C=input('pleaseenterthematrixC=');D=input('pleaseenterthematrixD=');U=input('pleaseentertheu=');T=input('pleaseentersamplingtimeT=');Th=input('pleaseenterthesimulationtimeTh=');x=input('pleaseentertheinitialvalueofx0=');y=input('pleaseentertheinitialvalueofy0=');E=input('pleaseentertheseterrorE=');[m,n]=size(A);[p,q]=size(B);7[s,r]=size(C);M=expm(A*T);disp('thestatetransitionmatrix')Ma=0;z=zeros(n);fori=1:1:nforj=1:1:na=a+abs(A(i,j));%a为A的范数endendfori=0:1:1z=z+(A*T)^i*T/prod(1:(i+1));endm1=det(z);L=1;while(1)M1=zeros(m,n);fori=0:1:LM1=M1+((A*T)^i*T/(prod(1:i)*(i+1)));endif(m1>det(M1))%求最小值m1=det(M1);endE=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;M1breakelseL=L+1;endendUc=B;Vo=C;k=rank(A);%求系统矩阵的维数disp('thedimensionofsystemis')kif(k==n)%判断系统是否满秩disp('thematrixofsystemisfullrank')endfori=1:n-18U1=A^i*B;V1=C*A^i;Uc=[Uc,U1];%可控性判断矩阵Vo=[Vo;V1];%可观性判断矩阵endif(n==rank(Uc))if(p==rank(Vo))disp('thesystemiscontrollableandobservable')elsedisp('thesystemiscontrollableandunobservable')endelseif(p==rank(Vo))disp('thesystemisuncontrollableandobservable')elsedisp('thesystemisuncontrollableandunobservable')endendX=zeros(m,(Th/T+1));Y=zeros(s,(Th/T+1));X(:,1)=x;Y(:,1)=y;t=0:T:Th;fori=2:Th/Tforj=1:mX(j,i)=M(j,:)*X(:,i-1)+M1(j,:)*U;%第j个状态变量X(j）在iT时刻的值endforj=1:sY(j,i)=C(j,:)*X(:,i)+D(j,:)*U;%第j个输出在iT时刻的值endendfori=1:m%前三行窗口画出状态变量X（t）subplot(6,2,i)plot(t,X(i,:));gridon;endforj=1:ssubplot(6,2,6+j)%第后三行窗口画出Y(t)plot(t,Y(j,:));gridon;end例如：激励u为0，系统矩阵A=[-0.300;2.9-0.62-2.3;02.30]，B=[100]'，9C=[110;1-31]，D=[01]'，初始条件为X（0）=[111]’，结果如图2.3和2.4运行结果如下：pleaseenterthematrixA=[-0.300;2.9-0.62-2.3;02.30]pleaseenterthematrixB=[100]'pleaseenterthematrixC=[110;1-31]pleaseenterthematrixD=[01]'pleaseentertheu=[0]pleaseentersamplingtimeT=0.1pleaseenterthesimulationtimeTh=10pleaseentertheinitialvalueofx0=[111]'pleaseentertheinitialvalueofy0=[2;-1]pleaseentertheseterrorE=1e-7thestatetransitionmatrixM=0.9704000.27450.9146-0.22110.03220.22110.9742M1=0.09850.01400.0011thedimensionofsystemisk=3thematrixofsystemisfullrankthesystemiscontrollableandobservable10图2.3系统可控可观判断图11图2.4X(t)和Y(t)仿真图2）下面通过Simulink建模仿真①建了一个Simulink框图模型保存名为job2s;②将输入模块，输入模块，分离模块，状态空间表达式模块分别复制到job2s文件中，并连接上；③在状态空间模块的参数对话框，将A,B,C,D分别填入参数对话框中的A,B,C,D四个矩阵输入参数栏，在初始条件栏中直接填入初始向量参数x(0)的值。图2.5仿真系统的框图模型12④在命令窗口输入[t,x,y]=sim('job2s',10)plot(t,y(:,1),':b',t,y(:,2),'-r');>>grid图2.6系统的输入曲线3）用MATLABLK库函数Uc=ctrb(A,B);Vo=obsv(A,C);求系统的可控和可观性代码可以在命令窗口输入执行：A=[-0.300;2.9-0.62-2.3;02.30];B=[100]';C=[110;1-31];D=[01]';k=rank(A);disp('therankofsystemis')kUc=ctrb(A,B);Vo=obsv(A,C);if(rank(Uc)==n)if(rank(Vo)==p)disp('thesystemiscontrolableandobservable')13elsedisp('thesystemiscontrolableandunobservable')endelseif(rank(Vo)==p)disp('thesystemisobservableisuncontrolable')elsedisp('thesystemisuncontrolableandunobservable')endend图2.7用可控可观库函数求解命令图3）结果比较通过比较用C代码实现的y（t）和用Simulink模块实现的y（t）进行比较可以看出用C代码实现的y（t）和Simulink模块实现的y（t）实现的接近。用C代码实现的可控可观性判断与用库函数实现的一致。用C代码实现也适用与多输入多输出系统，但在输入各个矩阵和初始值时要注意维数：A的行数和B的行数，C的列数，初值X(0)的列数要一致。B的列数，D的列数和U的行数要一致，D的行数和C的行数和Y(0)的列数要一致。多输入多输出例如下：pleaseenterthematrixA=[1032;2315;9843;2371]pleaseenterthematrixB=[13;45;86;12]pleaseenterthematrixC=[1234;5438]pleaseenterthematrixD=[01;23]14pleaseentertheu=[12]'pleaseentersamplingtimeT=0.1pleaseenterthesimulationtimeTh=10pleaseentertheinitialvalueofx0=[0000]'pleaseentertheinitialvalueofy0=[00]'pleaseentertheseterrorE=1e-7thestatetransitionmatrixM=1.35000.23670.52850.33900.47011.59180.46160.73141.47841.36662.01340.85070.73640.80931.10261.4584M1=0.32840.53060.66940.82821.42191.35760.63490.7396thedimensionofsystemisk=4thematrixofsystemisfullrankthesystemiscontrollableandobservable>>15图2.8多输入多输出实现图1617