(1)线性系统结构分析与分解及标准型

实验三十七 线性系统结构分析与分解及标准型

实验类型：验证 难度系数：0.3 实验性质：必做 课内学时：0 课外学时：2 分组人数：2

开课方式：在课外完成在MATLAB平台上完成实验。 实验目的：

掌握线性系统状态空间标准型、解及其模型转换。 实验设备与软件： 1、MATLAB数值分析软件 实验原理：

1、标准型变换、矩阵Jordan型变换阵、特征值

(1)标准型变换 命令格式 csys=canon(sys,’type’)

说明：type制定规范型的形式，包括两种选项：model(模态规范型)、companion(伴随规范型,友矩阵型,能控II型)。

(2)矩阵的Jordan规范型 命令格式 [V J]=Jordan(A) 说明：V特征向量，J是Jordan型

(3)求矩阵特征值和特征向量 命令格式 [V J]=eig(A) cv= eig(A) 说明：V特征向量，J是Jordan型;cv是特征值列向量 2、状态模型的相似变换： 命令格式 sysb=ss2ss(sys,T) 说明：sys是状态空间模型，T是非奇异变换阵的逆阵

? 传递函数模型与状态空间模型之间的相互转换：

命令格式 [A,B,C,D]=tf2ss(num,den)

[num,den]=ss2tf(A,B,C,D,iu)

说明：tf2ss:传递函数?状态空间; ss2tf状态空间?传递函数; iu是第iu个输入有效

? zpk模型与状态空间模型之间的相互转换：

命令格式：[A,B,C,D]=zp2ss(z,p,k)

[z,p,k]=ss2zp(A,B,C,D,iu)

说明：zp2ss :zpk模型?状态空间模型; ss2zp状态空间模型?zpk模型 3、线性定常系统的可控制与可观性及结构分解

关于可控制与可观性及结构分解的理论内容，请看教材学习。

(1)可控性和可观性----一般采用能控性/能观性矩阵类别(适用于离散或连续的情况) 状态可控性和输出可控性子函数如下： function str =pdctrb(A,B) function str=pdctrbo(A,B)%输出可控性 Qc=ctrb(A,B); Co=ctrb(A,B); r=rank(Qc); m=size(C,1);%返回行数 l=length(A); Qyc=[C*Co,D]; if r==l Tm=rank(Qyc);

if m==Tm str=’系统是状态完全可控的!’

else str=’系统不是状态完全可控的!’ end

实际上，rank(Qc)为系统的状态可控性指数，即系统中可控的状态的数目。

str=’系统输出是完全可控的!’ else str=’系统输出不是完全可控的!’ end

实际上，rank(Qyc) 为系统的输出可控性指数，即系统中可控的输出的数目。

状态可观性判别子函数代码如下： function str=pdobsv(A,C) Qo=obsv (A,C); r=rank(Qo); l= size(A,1); if r==l

str=’系统是状态完全可观的!’

else str=’系统不是状态完全可观的!’ end

实际上，rank(Qo)为系统的状态可观性指数，即系统中可观测的状态的数目。 有了上述子函数，在Matlab中可以直接调用这些子函数来判断可控性和可观性。 (2)可控性和可观性Gram矩阵可由下面的函数求得

W=gram(sys,type)----sys是系统的状态空间模型，type可以是’c’或’o’。

通过判断W的正定性判定其可控性或可观性。由于从W的数学表达式上看，Wc和Wo是对称的半正定矩阵，它们分别满足下面的Lyapunov方程(这种方程在稳定性分析中还将提到)

AWc?WcAT??BBT、AWo?WoAT??CTC

所以系统必须稳定才能得到Gram矩阵。

Wc矩阵中的值对应于输入信号对相应状态的贡献：第i个元素越大，则说明输入信号对第i个输入状态的贡献越大。

Wo矩阵中的值对应于每个状态对系统输出的贡献：第i个元素越大，则说明第i状态对系统输出的贡献越大。

(3)结构分解

a.在matlab中调用ctrbf()函数对系统按能控性分解 [Abar,Bbar,Cbar,T,K]=ctrbf(A,B,C)

[Abar,Bbar,Cbar,T,K]=ctrbf(A,B,C,TOL)

说明：K是可控的状态个数，TOL为误差，这里的变换是这样令的：x=T-1xbar，所以有Abar=TAT-1，Bbar=TB，Cbar=CT-1。显然这与教材中的令法是不同的，且T的选取方式也不同，得到的分解形式也不同，要加以区别。这里的能观性分解形式为

?o?(Ao,Bo,Co)：

?? Cc???AAc??c???A21?0?0??C?，Bc???，Cc???c???BAc??c?b.在matlab中调用ctrbf()函数对系统按能观性分解

[Abar,Bbar,Cbar,T,K]=obsvf(A,B,C)

[Abar,Bbar,Cbar,T,K]= obsvf(A,B,C,TOL)

说明：K是可观的状态个数，TOL为误差，这里的变换是这样令的：x=T-1xbar，所以有Abar=TAT-1，Bbar=TB，Cbar=CT-1。显然这与教材中的令法是不同的，且T的选取方式也不同，得到的分解形式也不同，要加以区别。这里的能观性分解形式为

?o?(Ao,Bo,Co)：

??AAo??o??0?????BA12o?? 0C?，Bo???，Co??o?????BA?o???o?c.按能控、能观性分解----Kalman分解

Kalman分解是先按能控性进行分解，然后对两部分状态再按能观性进行分解，最后得到能控能观部分、能控不能观部分、不能控能观部分、不能控不能观部分。

在Matlab中没有提供对应的函数对系统进行直接的Kalman分解，这里我们给出一个函数实现分解：

function [Gk,T,K]=kalmdec(G) G=ss(G);A=G.a;B=G.b;C=G.c; [Ac,Bc,Cc,Tc,Kc]=ctrbf(A,B,C); nc=rank(ctrb(A,B)); n=length(A); ic=n-nc+1:n;

[Ao1,Bo1,Co1,To1,Ko1]=obsv(Ac(ic,ic),Bc(ic),Cc(ic)); if nc

inc=1:n-nc;

[Ao2,Bo2,Co2,To2,Ko2]=obsvf(Ac(inc,inc),Bc(inc),Cc(inc)); end

[m1,n1]=size(To1); [m2,n2]=size(To2);

To=[To2,zeros(m2,n1);zeros(m1,n2),To1]; T=To*Tc;

n1=rank(obsv(Ac(ic,ic),Cc(ic))); n2= rank(obsv(Ac(inc,inc),Cc(inc)));

K=[zeros(1,n-nc-n2),ones(1,n2),…,2*ones(1,nc-n1),3*ones(1,n1)]; Ak=T*A*inv(T);Bk=T*B;Ck=C*inv(T); Gk=ss(Ak,Bk,Ck,G.d);

此函数分解的形式与教材中的是顺序不一样(变换的令法也不一样)：这里是不能控不能观、不能控能观、能控不能观、能控能观。这一点要注意。 4、定常线性系统的标准型(转换限于SISO系统)

(1)Jordan标准型

化成Jordan标准形是一种并联分解的策略，即将传递函数展开成部分分式和的形式。这种分解的原理在课堂上已讲的很清楚，请复习相关内容。这里就如何用Matlab求取此标准型做说明。

a.互异根的情况下，代码可以采用如下的形式，出可以有采用Jordan指令：

num=[…];den=[…];G=tf [num,den];

[r,p,k]=residue(num,den); A=diag(p);B=ones(length(r),1);

C=r’;%为了得到更准确的结果，可以采用C=(rat(r))’得到有理分式的形式 D=k;

Gss=SS(A,B,C,D) b.有重根的情况

我们直接可以采用Jordan指令。

num=[…]; den=[…]; Gtf=tf(num,den) Gs=ss(G)

[V J]=jordan(Gs.a) %V是特征向量，由特征向量组成非奇异变换矩阵, Gss=ss2ss(Gs,inv(V)) %Jordan型系统 (2)能控标准型-限于SISO系统

若系统能控，则可转换成能控标准I型和II型。 转换成能控标准II型代码：

A=[…];B=[…];C=[…];D=[…];Gs=ss(A,B,C,D);

T=ctrb(Gs.a,Gs.b) % x =Tz，T=[B,AB,A2B,…, A n-1B]， 能控标准II型 Abar=inv(T)*A*T;Bar=inv(T)*B;Cbar=C*T,Dar=D; Gss=ss(Abar,Bbar,Cbar,Dbar) 转换成能控标准I型代码：

A=[…];B=[…];C=[…];D=[…];Gs=ss(A,B,C,D); Tt=ctrb(Gs.a,Gs.b);Ttt=fliplr(Tt); cp=poly(Gs.a); n= length(Gs.a); Tea=eye(n) for i=2:n for j=1:(n-1) if i>j

Tea(i,j)=cp(i-(j-1)); end end end T=Ttt*Tea;

Abar=inv(T)*A*T;Bar=inv(T)*B;Cbar=C*T,Dar=D; Gss=ss(Abar,Bbar,Cbar,Dbar) (3)能观标准型-限于SISO系统

若系统能观，则可转换成能观标准I型和II型。 转换成能观标准I型代码:

A=[…];B=[…];C=[…];D=[…];Gs=ss(A,B,C,D);

x=Tz,T=V Tinv=obsv (Gs.a,Gs.b);T=inv(Tinv);

Abar=inv(T)*A*T;Bar=inv(T)*B;Cbar=C*T,Dar=D; Gss=ss(Abar,Bbar,Cbar,Dbar) 转换成能观标准II型代码

A=[…];B=[…];C=[…];D=[…];Gs=ss(A,B,C,D); Tt=obsv (Gs.a,Gs.c);Ttt==flipud(Tt); cp=poly(Gs.a); n= length(Gs.a); Tea=eye(n) for i=2:n for j=1:(n-1) if i>j

Tea(i,j)=cp(i-(j-1)); end end end Tea=Tea’; T= Tea *Ttt;

Abar=inv(T)*A*T;Bar=inv(T)*B;Cbar=C*T,Dar=D; Gss=ss(Abar,Bbar,Cbar,Dbar) 实验内容： 1、已知线性系统

??6?0.6250.75??1??x??0?u???8x00????

??20??0??0??y??1?0.250.0625?x?u(1) 判断其状态可控性、可观性和传递函数的关系，并加以说明分析。 (2) 对系统分别按能控性分解、能观性分解以及能控能观性分解。 2、在Matlab中建立并运行如下的.m代码，回答下面的问题。

num=[1 2 3];

den=conv([1 6 25],[1 12 35]); G=tf(num,den) Gs=ss(G)

[V J]=jordan(Gs.a) %求特征向量和Gs.a的Jordan标准型 Gss=ss2ss(Gs,inv(V)) %Jordan型系统 Gsm=canon(Gs,'model') %模态型系统 Gsf=canon(Gs,'companion') %能控标准II型系统

(1)给出无分号行的运行结果，并比较几个状态方程。 (2)在什么情况下,cannon得到的是对角型系统?请举例说明。

Gss=ss2ss(Gs,inv(V)) Gsm=canon(Gs,'model') Gsf=canon(Gs,'companion')

运行Untitled6.m

>> Untitled6 G =

s^2 + 2 s + 3 ------------------------------------

s^4 + 18 s^3 + 132 s^2 + 510 s + 875

Continuous-time transfer function. Gs = a =

x1 x2 x3 x4 x1 -18 -8.25 -3.984 -3.418 x2 16 0 0 0 x3 0 8 0 0 x4 0 0 2 0

b =

u1 x1 0.25 x2 0 x3 0 x4 0

c =

x1 x2 x3 x4 y1 0 0.25 0.0625 0.04688

d = u1 y1 0

Continuous-time state-space model. V =

-0.4883 -1.3398 0.4570 - 0.1719i 0.4570 + 0.1719i

1.5625 3.0625 -0.4375 + 1.5000i -0.4375 - 1.5000i -2.5000 -3.5000 -1.5000 - 2.0000i -1.5000 + 2.0000i

1.0000 1.0000 1.0000 1.0000 J =

-5.0000 0 0 0 0 -7.0000 0 0 0 0 -3.0000 - 4.0000i 0 0 0 0 -3.0000 + 4.0000i

Gss = a =

x1 x2 x3 x1 -5+1.04e-15i -2.13e-14+2.19e-15i 2.43e-14+2.27e-15i x2 -1.78e-15-4.62e-16i -7-1.06e-15i -5.7e-15+2.51e-15i x3 -1.38e-15+1.78e-15i 4.01e-15+2.22e-15i -3-4i x4 -1.73e-15+4.44e-16i 7.54e-15+4.44e-16i -4.88e-15+1.11e-15i

x4 x1 2.01e-14-4.83e-15i x2 -4.96e-15-2.82e-15i x3 0+1.78e-15i x4 -3+4i

b =

u1 x1 1.6-2.22e-16i x2 -1+1.67e-16i x3 -0.3-0.1i x4 -0.3+0.1i

c =

x1 x2 x3 x4 y1 0.281-1.77e-16i 0.594-2.78e-16i -0.156+0.25i -0.156-0.25i

d = u1 y1 0

Continuous-time state-space model.

Gsm =

a =

x1 x2 x3 x4 x1 -3 4 0 0 x2 -4 -3 0 0 x3 0 0 -7 0 x4 0 0 0 -5

b =

u1 x1 0.9492 x2 -0.1931 x3 4.942 x4 5.042

c =

x1 x2 x3 x4 y1 0.1699 0.09054 -0.1201 0.08925

d = u1 y1 0

Continuous-time state-space model.

Gsf = a =

x1 x2 x3 x4 x1 0 0 0 -875 x2 1 0 0 -510 x3 0 1 0 -132 x4 0 0 1 -18

b = u1 x1 1 x2 0

x3 0 x4 0

c =

x1 x2 x3 x4 y1 0 1 -16 159 d = u1 y1 0

Continuous-time state-space model.

(2)当系统的 | sI - A |= 0 得到的特征多项式的解无重根时，则cannon得到的是对角型系统。 例如：

>> A=[0 1 0;0 0 1;-6 -11 -6]; B=[0 ;0; 1]; C=[1,0,0]; D=[0]; E=eig(A)

sys=ss(A,B,C,D);

[Gt,P]=canon(sys,'model')

E =

-1.0000 -2.0000 -3.0000 Gt =

a =

x1 x2 x3 x1 -3 0 0 x2 0 -2 0 x3 0 0 -1

b =

u1 x1 -3.881 x2 -4.899 x3 1.436

c =

x1 x2 x3 y1 -0.1288 0.2041 0.3482

d = u1 y1 0

Continuous-time state-space model. P =

-7.7621 -11.6431 -3.8810 -14.6969 -19.5959 -4.8990

8.6168 7.1807 1.4361 (3)

能控标准II型 程序：

function Gss2 = NK2(A,B,C,D) d6;

A=Gs.a;B=Gs.b;C=Gs.c;D=Gs.d; Gs2=ss(A,B,C,D); T=ctrb(Gs.a,Gs.b); Abar=inv(T)*A*T; Bbar=inv(T)*B; Cbar=C*T; Dbar=D;

Gss2=ss(Abar,Bbar,Cbar,Dbar) end

结果： >> NK2 Gss2 = a =

x1 x2 x3 x4 x1 0 0 0 -875 x2 1 0 0 -510 x3 0 1 0 -132 x4 0 0 1 -18

b = u1 x1 1 x2 0 x3 0 x4 0 c =

x1 x2 x3 x4 y1 0 1 -16 159

d = u1 y1 0

Continuous-time state-space model.

能控标准I型： 程序：

function Gss1 = NK1(A,B,C,D) d6;

A=Gs.a;B=Gs.b;C=Gs.c;D=Gs.d; Gs1=ss(A,B,C,D); Tt=ctrb(A,B); Ttt=fliplr(Tt); cp=poly(A); n=length(A); Tea=eye(n); for i=2:n for j=1:(n-1) if i>j

Tea(i,j)=cp(i-(j-1)); end end end T=Ttt*Tea; Abar=inv(T)*A*T; Bbar=inv(T)*B; Cbar=C*T; Dbar=D;

Gss1=ss(Abar,Bbar,Cbar,Dbar) end

结果： >> NK1

Gss1 =

a =

x1 x2 x3 x4 x1 2.842e-14 1 x2 -2.842e-14 0 x3 -1.023e-12 0 x4 -875 -510

b = u1 x1 0 x2 0 x3 0 x4 1

c =

x1 x2 x3 x4 y1 3 2 1 0

d = u1 y1 0

Continuous-time state-space model.

能观I型： 程序：

function Gss3 = NG1(A,B,C,D) d6;

A=Gs.a;B=Gs.b;C=Gs.c;D=Gs.d; Gs3=ss(A,B,C,D);

Tinv=obsv(Gs.a,Gs.c); T=inv(Tinv);

Abar=inv(T)*A*T; Bbar=inv(T)*B; Cbar=C*T;

0 0 1 0 0 1 -132 -18 Dbar=D;

Gss3=ss(Abar,Bbar,Cbar,Dbar) end

结果： >> NG1 Gss3 =

a =

x1 x2 x3 x4 x1 -1.176e-14 1 -4.677e-16 -3.431e-17 x2 1.99e-13 7.105e-14 1 7.772e-16 x3 -3.183e-12 -9.095e-13 -1.137e-13 1 x4 -875 -510 -132 -18

b =

u1 x1 6.896e-17 x2 1 x3 -16 x4 159

c =

x1 x2 x3 x4 y1 1 -3.331e-16 -5.551e-17 -3.469e-18

d = u1 y1 0

Continuous-time state-space model.

能观II型： 程序：

function Gss4 = NG2(A,B,C,D) d6;

A=Gs.a;B=Gs.b;C=Gs.c;D=Gs.d; Gss4=ss(A,B,C,D); Tt=obsv(A,C); Ttt=fliplr(Tt); cp=poly(A);

n=length(A); Tea=eye(n); for i=2:n for j=1:(n-1) if i>j

Tea(i,j)=cp(i-(j-1)); end end end Tea=Tea'; T=Ttt*Tea; Abar=inv(T)*A*T; Bbar=inv(T)*B; Cbar=C*T; Dbar=D;

Gss4=ss(Abar,Bbar,Cbar,Dbar) end

结果： >> NG2 Gss4 =

a =

x1 x2 x3 x4 x1 -1.877e+06 -3.572e+07 -2.86e+08 -1.279e+09 x2 1.699e+05 3.233e+06 2.588e+07 1.158e+08 x3 -9914 -1.887e+05 -1.511e+06 -6.757e+06 x4 227.2 4324 3.462e+04 1.548e+05

b =

u1 x1 676 x2 -61.17 x3 3.57 x4 -0.08175

c =

x1 x2 x3 x4 y1 9.399 178.9 1432 6405

d =

u1 y1 0

Continuous-time state-space model.

分析:能控型I型与能观II型对偶；能控II型与能观I型对偶。

实验总结

本实验主要以MATLAB为工具对线性系统的能关能控进行判断和分析，包括对系统进行能控能观性分解和对模型进行交换。研究了能控和能观标准型的关系。通过这次实验，我熟悉了MATLAB对线性系统的基本操作，重新复习了对M文件的编写方法，对系统的能控能管性有了更深的了解，收货很大。

本文来源：https://www.bwwdw.com/article/609v.html