Matlab程序解现代控制理论与工程中的状态方程_免费下载.doc

用矩阵指数法解状态方程的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,taoPhit=expm(A*t); %求矩阵指数exp(At)if (B=0) B=zeros(size(A,l),l); %重构系数矩阵Bendphi=sub(Phit,t,t-tao); %求expA(t-tao)PhitBu=int(phi*B*ut,tao,0,t); %求expA(t-tao)*B*u(tao)在0t区间的积分用拉氏变换法解状态方程的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,taoAA=s*eye(size(A)-A; %求sI-AinvAA=inv(AA); %求（sI-A）矩阵的逆intAAtAA=ilaplace(intAA) ; %求intAA的拉氏反变换sI_A=simplify; %简化拉式反变换的结果if (B=0) B=zeros(size(A,l),l); %重构系数矩阵BendtAB=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)endutao=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=Ifor 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); %修正Phiend求解线性定常离散系统状态方程的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)zZ反变换后的结果%AkBu输出矩阵(zI-A)(-1)*B*u(z)Z反变换后的结果syms z %定义符号变量zAA=z*eye(size(A)-A; %求zI-AinvAA=inv(AA); %求（zI-A）矩阵的逆intAAtAA=iztrans(intAA*z) ; %求intAA*z的Z反变换Ak=simple(tAA); %简化Z反变换的结果if (B=0) B=zeros(size(A,l),l); %重构系数矩阵BendtAB=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 %定义符号变量kif (Bk=0) Bk=zeros(size(A,l),l); %重构系数矩阵Bend 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;endAAB=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 TPhit=expm(A*t); %求矩阵指数exp(At)if (B=0) B=zeros(size(A,l),l); %重构系数矩阵BendPhitB=int(Phit*B,t,0,T); %求exp(At)*B在0T区间的积分 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 符号变量，指明矩阵AB中的时变参数，通常为时间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)endPhitB=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)%ccanonl 求系统X=AX+Bu，y=Cx的可控规范I型系数矩阵%Abar,Bbar,Cbar,变换后的可控规范I型系数矩阵%T 相似变换矩阵n=length(A);Co=ctrb(A,B);if (rank(Co)=n), %判断系统可控性 error(系统不可控！);endRs=sym(polymtx(A); %计算R矩阵并转变为符号矩阵Rs形式As=sym(A); %讲矩阵A转变为符号矩阵AsBs=sym(B); %讲矩阵B转变为符号矩阵BsTs=Bs; for i=1:n-1, Ts=Asi*Bs Ts;end Ts=Ts*Rs; %计算相似变换符号矩阵TsAbar=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|=sn+a(n-1)*s(n-1)+.+a(1)*s+a(0)%的各项系数d=size(A);if (length(d)=2), %判断系统可控性 error(错误：非二维矩阵！);endif (d(1)=d(2), %判断系统可控性 error(错误：A非方针！);endn=d(1);As=sym(A);p=sym2poly(poly(As);R=;for i=1:n R=R p(1:n); p=0 p(1:n); endR=numeric(R);线性定常系统可控分解函数cdecomp：函数cdecomp: 线性定常系统的可控性分解function Abar,Bbar,Cbar,P=cdecomp(A,B,C)%cdecomp 可控性分解%Abar,Bbar,Cbar,P=cdecomp(A,B,C)%若系统不完全可控，则存在相似变换矩阵P，使得% -1 -1 %Abar=P*A*P, Bbar=P*B, Cbar=C*P%其中%Abar=|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); %求可控性矩阵Mn=numeric(rank(Ms);if (nnA), P=; %系统不可控，计算变换矩阵P i=1; while numeric (rank(P)n, %取M的n个线性无关列矢量至P矩阵 P=P Ms(:,i); if (numeric(rank(P)size(P,2), P (:,size(P,2)= ; end i=i+1;endE=sym(eye(size(