Matlab程序解现代控制理论与工程中的状态方程.doc
《Matlab程序解现代控制理论与工程中的状态方程.doc》由会员分享,可在线阅读,更多相关《Matlab程序解现代控制理论与工程中的状态方程.doc(8页珍藏版)》请在装配图网上搜索。
用矩阵指数法解状态方程的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) %ccanonl 求系统X’=AX+Bu,y=Cx的可控规范I型系数矩阵 %Abar,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) %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)); %求可控性矩阵M n=numeric(rank(Ms)); if (n- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- Matlab 程序 现代 控制 理论 工程 中的 状态方程

链接地址:https://www.zhuangpeitu.com/p-9415583.html