线性定常连续系统状态方程的解

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）的分数symsttao%定义符号变量t,taophit=expm(a*t);%求矩阵指数exp(at)if(b==0)b=zeros(size(a,l),l);%重构系数矩阵bendphi=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)拉式反变换后的结果symss%定义符号变量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――THF1求解分量%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));%排序THF1分量求解时变系统转移矩阵的matlab函数transmtx：函数transmtx:解线性时变系统状态迁移矩阵functionphi=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=iforlop=0:naa=a;fori=1:lopif(aa==0)break;endatemp=subs(aa,x,’taoi’);aa=simplify(a*int(atemp,’tao’,a,x));endif(aa==0)break;endatemp=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)z]z反变换后的结果%akbu――输入矩阵[((zi-a)^(-1)*b*u(z)]z反华转换后的结果symsz%定义符号变量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:解线性时变离散系统状态方程的求解functionxk=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)....symsk%定义符号变量kif(bk==0)bk=zeros(size(a,l),l);%重构系数矩阵bendxk=[];forkk=kstart+1:kendaa=eye(size(k));fori=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));fori=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;endb=subs(bk,’k’,kk-1);u=subs(uk,’k’,kk-1);bb=bb+b*u;xk=[xkaa*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为取样周期symstt%定义符号变量ttphit=expm(a*t);%求矩阵指数exp(at)if(b==0)b=zeros(size(a,l),l);%重构系数矩阵bendphitb=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时包含二重积分项，.....symsttphit=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);elsebtao=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)ìanonl求系统x’=ax+bu，y=cx的可控规范i型系数矩阵?ar,bbar,cbar,――变换后的可控规范i型系数矩阵%t――相似变换矩阵n=length(a);co=ctrb(a,b);if(rank(co)~=n),%判断系统可控性error(‘系统不可控！’);。

第三章 系统的分析——状态方程的解§3-1线性连续定常齐次方程求解一、齐次方程和状态转移矩阵的定义１、齐次方程状态方程的齐次方程部分反映系统自由运动的状况（即没有输入作用的状况），设系统的状态方程的齐次部分为：)()(t Ax t x= 线性定常连续系统：Ax x= 初始条件：00x x t ==２、状态转移矩阵的定义齐次状态方程Ax x = 有两种常见解法：（1）幂级数法；（2）拉氏变换法。

其解为)0()(x e t x At ⋅=。

其中At e 称为状态转移矩阵（或矩阵指数函数、矩阵指数），记为：At e t =)(φ。

若初始条件为)(0t x ，则状态转移矩阵记为：)(00)(t t A e t t -=-Φ 对于线性时变系统，状态转移矩阵写为),(0t t φ，它是时刻t ，t 0的函数。

但它一般不能写成指数形式。

（１）幂级数法——直接求解设Ax x= 的解是t 的向量幂级数 +++++=k k t b t b t b b t x 2210)(式中 ，，，，，k b b b b 210都是n 维向量，是待定系数。

则当0=t 时， 000b x x t ===为了求其余各系数，将)(t x 求导，并代入)()(t Ax t x = ，得：+++++=-1232132)(k k t kb t b t b b t x)(2210 +++++=k k t b t b t b b A上式对于所有的t 都成立，故而有：⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧======00323021201!1!31312121b A k b b A Ab b b A Ab b Ab b K K且有：00x b =故以上系数完全确定，所以有：+++++=k k t b t b t b b t x 2210)(+++++=k k t b A k t b A t Ab b 020200!1!21)0()!1!21(22x t A k t A At I kk +++++=定义（矩阵指数或矩阵函数）：∑∞==+++++=022!1!1!21K kk k k AttA k t A k t A At I e则)0()(x e t x At⋅=。

实验二 线性连续定常系统的运动分析一、实验目的1．掌握线性连续定常系统的状态转移矩阵的求法，学会用MATLAB 求解状态转移矩阵。

2．掌握线性连续定常系统的状态方程的求解方法，学会用MATLAB 求解线性连续定常系统的时间响应，并绘制相应的状态响应曲线和输出响应曲线。

二、实验原理1．线性连续定常系统状态转移矩阵的计算设线性连续定常系统的状态空间表达式为'=+⎧⎨=+⎩x Ax Buy Cx Du ，则其状态转移矩阵为()t t e =A Φ从时间角度看，状态转移矩阵使状态向量随着时间的推移不断地作坐标变换，不断地在状态空间中作转移，故称为状态转移矩阵。

对于线性连续定常系统，其状态转移矩阵与其矩阵指数函数相同，可利用直接求解法、拉氏变换法、标准型法和待定系数法等方法对其进行求解。

（1）直接求解法220111()!2!!tk k k kk t e t t t t k k ∞====+++++∑A A I A A A Φ（2）拉氏变换法()11()t t e L s --⎡⎤==-⎣⎦A I A Φ（3）标准型法对系统矩阵A 进行线性非奇异变换，将其变换为对角线矩阵或约旦矩阵1-=A P AP ，其中P 为非奇异变换阵。

状态转移矩阵为1()t t t e e -==A A P P Φ，其中1-=A P AP若A 的特征值12,,,n λλλ两两互异，则A 为对角线矩阵，此时1110()0n t tt t e t e e e λλ--⎡⎤⎢⎥===⎢⎥⎢⎥⎣⎦A A P P P P Φ 若A 有n 重特征值i λ，则A 为约旦矩阵，此时1111(1)!()0i i i i i ttt n t t tte te t e n t e e te e λλλλλ---⎡⎤⎢⎥-⎢⎥===⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦A A Q Q Q Q Φ （4）待定系数法根据凯莱－哈密顿（以下简称C-H ）定理，线性连续定常系统的状态转移矩阵为110110()()()()()n tj n j n j t ea t a t a t a t ---====+++∑A A I A A Φ其中，011(),(),,()n a t a t a t -为t 的标量函数，可按A 的特征值确定。

连续系统的状态变量方程求解连续系统的状态变量方程求解通常采用数值方法，例如龙格-库塔法（Runge-Kutta）等。

在这个过程中，需要将连续系统的状态方程离散化，即将连续时间步长的微分方程转化为离散时间步长的离散方程。

求解离散方程可采用递推的方式，根据系统的初始条件和上一时刻的状态变量值，计算出当前时刻的状态变量值。

以下是一个求解连续系统状态变量方程的步骤：1. 确定连续系统的状态变量方程。

例如，给定线性定常系统dx/dt = Ax + Bu，其中x为状态变量，A和B为系统矩阵。

2. 离散化。

将状态变量方程转化为离散方程。

常见的离散化方法有前项差分变换、后项差分变换和Tustin变换。

具体变换方法取决于系统的特性以及所需的数值稳定性和精度。

例如，使用Tustin变换将连续系统离散化，得到离散状态方程x[k+1] = A*x[k] + B*u[k]。

3. 初始化。

给定初始条件，如x[0] 和u[0]，初始化状态变量值。

4. 数值求解。

使用数值方法（如龙格-库塔法）递推计算离散方程，得到一系列状态变量值x[1], x[2], ...，以及对应的输出值y[1], y[2], ...。

5. 分析结果。

根据求解得到的状态变量值和输出值，分析系统的性能，如稳定性、收敛速度等。

在MATLAB中，可以使用ode45等函数求解连续系统的状态变量方程。

以下是一个简单的示例：```MATLAB定义系统矩阵A、B和输入信号uA = [1 0; -1 1];B = [0 1];u = [1; 0.5];定义初始条件x0 = [1; 2];设置求解参数tspan = [0, 10];options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);求解状态变量方程[x, u] = ode45(@(t, x) A*x + B*u, tspan, x0, options);绘制状态变量曲线figure;plot(t, x(:, 1), 'b', 'LineWidth', 2);hold on;plot(t, x(:, 2), 'r', 'LineWidth', 2);xlabel('Time');ylabel('State Variables');legend('x1', 'x2');```这个示例中，我们使用ode45函数求解了一个线性定常系统在给定输入信号下的状态变量演化。

１．状态空间表达式n 阶DuCx y Bu Ax x+=+= 1:⨯r u 1:⨯m y n n A ⨯:r n B ⨯:n m C ⨯:rm D ⨯:A 称为系统矩阵，描述系统内部状态之间的联系；Ｂ为输入（或控制）矩阵，表示输入对每个状态变量的作用情况；C 输出矩阵，表示输出与每个状态变量间的组成关系，Ｄ直接传递矩阵，表示输入对输出的直接传递关系。

２．状态空间描述的特点①考虑了“输入－状态－输出”这一过程，它揭示了问题的本质，即输入引起了状态的变化，而状态决定了输出。

②状态方程和输出方程都是运动方程。

③状态变量个数等于系统包含的独立贮能元件的个数，n 阶系统有n 个状态变量可以选择。

④状态变量的选择不唯一。

⑤从便于控制系统的构成来说，把状态变量选为可测量或可观察的量更为合适。

⑥建立状态空间描述的步骤：a 选择状态变量；b 列写微分方程并化为状态变量的一阶微分方程组；c 将一阶微分方程组化为向量矩阵形式，即为状态空间描述。

⑦状态空间分析法是时域内的一种矩阵运算方法，特别适合于用计算机计算。

３．模拟结构图（积分器加法器比例器）已知状态空间描述，绘制模拟结构图的步骤：积分器的数目应等于状态变量数，将他们画在适当的位置，每个积分器的输出表示相应的某个状态变量，然后根据状态空间表达式画出相应的加法器和比例器，最后用箭头将这些元件连接起来。

４．状态空间表达式的建立1由系统框图建立状态空间表达式：a 将各个环节（放大、积分、惯性等）变成相应的模拟结构图；b 每个积分器的输出选作i x ，输入则为i x；c 由模拟图写出状态方程和输出方程。

2由系统的机理出发建立状态空间表达式：如电路系统。

通常选电容上的电压和电感上的电流作为状态变量。

利用KVL 和KCL 列微分方程，整理。

③由描述系统的输入输出动态方程式（微分方程）或传递函数，建立系统的状态空间表达式，即实现问题。

实现是非唯一的。

方法：微分方程→系统函数→模拟结构图→状态空间表达式。