Matlab程序解现代控制理论与工程中的状态方程

用矩阵指数法解状态方程的MATLAB函数vslove1：

函数vslove1:求解线性定常连续系统状态方程的解 function [Phit,PhitBu]=vsolves1(A,B,ut)

%vsolves1求线性连续系统状态方程X’=AX+Bu的解 %[Phit,phitBu]=vsolves1(A,B,ut) %A,B系数矩阵

%ut控制输入，必须为时域信号的符号表达式，符号变量为t %Phit——输出Phi(t)

%PhitBu——输出phi(t-tao)*B*u(tao)在区间（0，t）的积分

syms t tao %定义符号变量t,tao Phit=expm(A*t); %求矩阵指数exp(At) if (B==0)

B=zeros(size(A,l),l); %重构系数矩阵B end

phi=sub(Phit,’t’,’t-tao’); %求exp[A(t-tao)]

PhitBu=int(phi*B*ut,’tao’,’0’,’t’); %求exp[A(t-tao)]*B*u(tao)在0~t区间的积分

用拉氏变换法解状态方程的MATLAB函数vslove2：

函数vslove2:求解线性定常连续系统状态方程的解 function [sl_A,sl_ABu]=vsolves1(A,B,us)

%vsolves2求线性连续系统状态方程X’=AX+Bu的解 %[sl_A,sl_ABu]=vsolves1(A,B,ut) %A,B系数矩阵

%us控制输入，必须为拉氏变换后的符号表达式，符号变量为s %sl_A——输出矩阵(sl-A)^(-1)拉式反变换的结果

%sl_ABu——输出(sl-A)^(-1)*B*u(s)拉式反变换后的结果

syms s %定义符号变量t,tao AA=s*eye(size(A))-A; %求sI-A

invAA=inv(AA); %求（sI-A）矩阵的逆intAA tAA=ilaplace(intAA) ; %求intAA的拉氏反变换 sI_A=simplify; %简化拉式反变换的结果 if (B==0)

B=zeros(size(A,l),l); %重构系数矩阵B end

tAB=ilaplace(intAA*B*us) ; %求intAA*B*us的拉氏反变换 sI_ABu=simplify(tAB); %化简拉式反变换的结果

求解时变系统状态方程的MATLAB函数tslove：

函数tslove:求解线性时变连续系统状态方程的解 function [Phi,PhiBu]=tsolves(A,B,u,x,a,n) %tsolves求时变系统状态方程 %[Phi,phiBu]=vsolves1(A,B,u,x,a,n) %A,B时变系数矩阵

%Phi——状态转移矩阵计算结果 %PhiBu——受控解分量

%u——控制输入向量，时域形式

%x——符号变量，指明矩阵A中的时变参数，通常为时间t %a——积分下限

%n——时变状态转移矩阵中计算重积分的最大项数，n=0时无重积分项 % n=1时包含二重积分项，.....

Phi=transmtx(A,x,a,n); %计算状态转移矩阵 Phitao=subs(Phi,x,’tao’); %求Phi(tao)

if (B==0)

Btao=zeros(size(A,l),l); %求B(tao) end

utao=subs(u,x,’tao’); %求u(tao)

PhiBu=simple(int(Phitao*Btao*utao,’tao’,a,x)); %计算受控分量

求解时变系统转移矩阵的MATLAB函数transmtx：

函数transmtx:求解线性时变系统状态转移矩阵 function Phi=transmtx(A,x,a,n)

%transmtx计算时变系统状态转移矩阵 %Phi=transmtx(A,x,a,n)

%Phi——状态转移矩阵计算结果 %A 时变系数矩阵

%x——符号变量，指明矩阵A中的时变参数，通常为时间t %a——积分下限

%n——时变状态转移矩阵中计算重积分的最大项数，n=0时无重积分项 % n=1时包含二重积分项，.....

Phi=eye(size(A)); %初始化Phi=I for lop=0:n AA=A; for i=1:lop if (AA==0) break; end

Atemp=subs(AA,x,’taoi’);

AA=simplify(A*int(Atemp,’tao’,a,x));

end

if (AA==0) break; end

Atemp=subs(AA,x,’taoi’);

AA=simplify(A*int(Atemp,’tao’,a,x)); %计算重积分 Phi=simplify(Phi+AA); %修正Phi end

求解线性定常离散系统状态方程的MATLAB函数disolve：

函数disolve:求解线性定常离散系统状态方程的解 function [Ak,AkBu]=disolve(A,B,uz)

%disolve 求线性离散系统状态方程x(k+1)=Ax(k)+Bu(k)的解 %[Ak,AkBu]=disolve(A,B,uz) %A,B 系数矩阵

%uz 控制输入，必须为Z变换后的符号表达式，符号变量为z %Ak——输出矩阵[((zI-A)^(-1)z]Z反变换后的结果

%AkBu——输出矩阵[((zI-A)^(-1)*B*u(z)]Z反变换后的结果

syms z %定义符号变量z AA=z*eye(size(A))-A; %求zI-A

invAA=inv(AA); %求（zI-A）矩阵的逆intAA tAA=iztrans(intAA*z) ; %求intAA*z的Z反变换 Ak=simple(tAA); %简化Z反变换的结果 if (B==0)

B=zeros(size(A,l),l); %重构系数矩阵B end

tAB=iztrans(intAA*B*uz) ; %求intAA*B*uz的Z反变换 AkBu=simple(tAB); %化简Z反变换的结果

求解线性时变离散系统状态方程的MATLAB函数tdsolve：

函数tdsolve:求解线性时变离散系统状态方程的解 function xk=tsolve(Ak,Bk,uk,x0,kstart,kend)

%tdsolve 求线性时变离散系统状态方程x(k+1)=A(k)x(k)+B(k)u(k)的解 %xk=tsolve(Ak,Bk,uk,x0,kstart,kend) %Ak,Bk 系数矩阵

%uk 控制输入，必须为时域符号表达式，符号变量为k %x0 初始状态 %kstart——初始时刻 %kend——终止时刻

%xk——输出结果，矩阵每一列分别对应x(k0+1),x(k0+2)....

syms k %定义符号变量k if (Bk==0)

Bk=zeros(size(A,l),l); %重构系数矩阵B end

xk=[] ;

for kk=kstart+1: kend AA=eye(size(k));

for i=kstart : kk-1 %计算A(k-1)A(k-2)....A(k0+1)A(k0) A=subs(Ak,’k’,i); AA=A*AA; end

AAB=eye(size(Ak)); BB=zeros(size(Bk));

for i=kk-1 :-1:kk+1 %计算A(k-1)A(k-2)....A(j+1)B(j)u(j)的累加和 A=subs(Ak,’k’,i); AAB=AAB*A;

B=subs(Bk,’k’,kk-1+i+kstart); u=subs(uk,’k’,kk-1+i+kstart); BB=BB+AAB*B*u; end

B=subs(Bk,’k’,kk-1); u=subs(uk,’k’,kk-1); BB=BB+B*u;

xk=[xk AA*x0+BB]; %计算x(k) end

连续系统状态方程离散化的MATLAB符号函数sc2d：

函数 sc2d: 线性连续系统状态方程的离散化 function [Ak,Bk]=sc2d(A,B)

%sc2d 离散化线性连续系统状态方程X’=AX+Bu %sysd=sc2d(A,B)

%A,B ——连续系统的系数矩阵 %Ak,Bk ——离散系统系数符号矩阵

% 离散状态方程为：x(k+1)=Ak*x(k)+Bk*u(k) % Ak,Bk中变量T为采样周期

syms t T %定义符号变量t T Phit=expm(A*t); %求矩阵指数exp(At) if (B==0)

B=zeros(size(A,l),l); %重构系数矩阵B end

PhitB=int(Phit*B,’t’,0,’T’); %求exp(At)*B在0~T区间的积分 Ak=simple(subs(Phit,’t’,’T’)); Bk=simple(PhitB);

线性时变系统离散化的MATLAB函数tc2d：

函数 tc2d: 线性时变系统的离散化 function [Ak,Bk]=tc2d(A,B,x,n) %tc2d 线性时变系统的离散化 %[Ak,Bk]=tc2d(A,B,x,n)

%A,B ——连续系统的系数矩阵

%Ak ——离散化后的系数矩阵A(kt) %Bk ——离散化后的系数矩阵B(kt)

%x ——符号变量，指明矩阵A\\B中的时变参数，通常为时间t

%n ——时变状态转移矩阵中计算重积分的最大项数，n=0时无重积分项， % n=1时包含二重积分项，.....

syms t T

Phit=transmtx(A,x,k*T,n); %计算时变系统的状态转移矩阵 Ak=simplify(subs(Phi,x,(k+1)*T)); %计算离散化后的系数矩阵A(kT)

Phitao=subs(Phi,x,’tao’); %求Phi(tao) if (B==0)

Btao=zeros(size(A,l),l); else

Btao=subs(B,x,’tao’); %求B(tao) end

PhitB=simple(int(Phitao*Btao,’tao,k*T,x’)); %计算受控分量

Bk=simplify(subs(PhiB,x,(k+1)*T)); %计算离散化后的系数矩阵B(kT)

定常系统可控规范I型变换函数ccanonl：

函数ccanonl: 求线性定常系统的可控规范I型形式 function [Abar,Bbar,Cbar,T]=ccanonl(A,B,C)

ìanonl 求系统X’=AX+Bu，y=Cx的可控规范I型系数矩阵 ?ar,Bbar,Cbar,——变换后的可控规范I型系数矩阵 %T ——相似变换矩阵

n=length(A); Co=ctrb(A,B);

if (rank(Co)~=n), %判断系统可控性 error(‘系统不可控！’);

end

Rs=sym(polymtx(A)); %计算R矩阵并转变为符号矩阵Rs形式 As=sym(A); %讲矩阵A转变为符号矩阵As Bs=sym(B); %讲矩阵B转变为符号矩阵Bs Ts=Bs; for i=1:n-1,

Ts=[As^i*Bs Ts]; end

Ts=Ts*Rs; %计算相似变换符号矩阵Ts

Abar=numeric(inv(Ts)*As*Ts); %实现矩阵A的相似变换并转变为数值形式 Bbar=numeric(inv(Ts)*Bs); %实现矩阵B的相似变换并转变为数值形式 Cbar=numeric(inv(C)*Ts); %实现矩阵C的相似变换并转变为数值形式 T=numeric(Tc); %相似变换矩阵T转变为数值形式

求系统相似变换的R矩阵函数polymtx：

函数polymtx: 求系统相似矩阵变换的R矩阵 function R=polymtx(A)

%R=polymtx(A)——A须为方针

% [1 0 0 0] % [a(n-1) 1 0 0] % R= [a(n-2) a(n-1) 0 0] % [... ... ... ...] % [a(2) a(3) 1 0] % [a(1) a(2) a(n-1) 1]

%其中a(i)(i=1,...n-1)为矩阵A的特征多项式 %|s*I-A|=s^n+a(n-1)*s^(n-1)+...+a(1)*s+a(0) %的各项系数

d=size(A);

if (length(d)~=2), %判断系统可控性 error(‘错误：非二维矩阵！’); end

if (d(1)~=d(2)), %判断系统可控性 error(‘错误：A非方针！’); end n=d(1); As=sym(A);

p=sym2poly(poly(As)); R=[]; for i=1:n

R=[R p(1:n)’];

p=[0 p(1:n)’]; end

R=numeric(R);

线性定常系统可控分解函数cdecomp：

函数cdecomp: 线性定常系统的可控性分解 function [Abar,Bbar,Cbar,P]=cdecomp(A,B,C) íecomp 可控性分解

%[Abar,Bbar,Cbar,P]=cdecomp(A,B,C)

%若系统不完全可控，则存在相似变换矩阵P，使得 % -1 -1

?ar=P*A*P, Bbar=P*B, Cbar=C*P %其中

?ar=|Ac A12| ,Bbar=|Bc |,Cbar=|Cc Cuc| |0 Auc| |0|

%(Ac,Bc)构成系统的可控子空间

As=sym(A); %转变为符号矩阵求解 Bs=sym(B); Cs=sym(C); nA=size(As,l);

Ms=sym(ctrb(As,Bs)); %求可控性矩阵M

n=numeric(rank(Ms)); if (n

P=[]; %系统不可控，计算变换矩阵P i=1;

while numeric (rank(P))

if (numeric(rank(P))

E=sym(eye(size(A))); i=1;

while numeric (rank(P))

if (numeric(rank(P))

end else

P=eye(size(A)); %若系统可控，取P=1 end

Abar=numeric(inv(P)*A*P); %转变为数值矩阵输出 Bbar=numeric(inv(P)*B); Cbar=numeric(inv(C*P)); P=numeric(P);

线性定常系统可观分解函数odecomp：

函数odecomp: 线性定常系统的可观性分解 function [Abar,Bbar,Cbar,P]=odecomp(A,B,C) %odecomp可控性分解

%[Abar,Bbar,Cbar,P]=odecomp(A,B,C)

%若系统不完全可观，则存在相似变换矩阵P，使得 % -1 -1

?ar=P*A*P, Bbar=P*B, Cbar=C*P %其中

?ar=|Ao 0| ,Bbar=|Bo |,Cbar=|Co 0| |A21 Ano| |Bno| %(Ao,Bo)构成系统的可观子空间

%根据对偶原理，应用可控性分解函数cdecomp实现可观性分解

[abar,bbar,cbar,P]=cdecomp(A’,B’,C’); Abar=abar’; Bbar=cbar’; Cbar=bbar’;

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