线性系统理论的matlab应用.

MATLAB的应用曾水平控制系统的分析方法早期的控制系统分析过程复杂而耗时，如想得到一个系统的冲激响应曲线，首先需要编写一个求解微分方程的子程序，然后将已经获得的系统模型输入计算机，通过计算机的运算获得冲激响应的响应数据，然后再编写一个绘图程序，将数据绘制成可供工程分析的响应曲线。MATLAB控制系统工具箱和SIMULINK辅助环境的出现，给控制系统分析带来了福音。控制系统的分析包括系统的稳定性分析、时域分析、频域分析及根轨迹分析。本课件是关于状态空间的时域分析一利用MATLAB进行模型转换1-1传递函数与状态空间表达式之间的转换1）.连续系统状态空间表达式MATLAB是当今世界上最优秀的科技应用软件之一，它以强大的科学计算能力和可视化功能，简单易用的编程语言以及开放式的编程环境等一些显著的优点，使得它在当今许许多多科学技术领域中成为计算机辅助分析和设计、算法研究和应用开发的基本工具和首选平台。在本书中，用它作为系统分析和设计的软件平台，更显示出独特的优势。本节利用MATLAB实现数学模型的转换。可以用ss命令来建立状态空间模型。对于连续系统，其格式为sys=ss(A,B,C,D)，其中A，B，C，D为描述线性连续系统的矩阵。当sys1是一个用传递函数表示的线性定常系统时，可以用命令sys=ss(sys1)，将其转换成为状态空间形式。也可以用命令sys=ss(sys1,’min’)计算出系统sys的最小实现。例1控制系统微分方程为uuuuyyyyy2424724503510)4(求其状态空间表达式。解可以先将其转换成传递函数2450351024247)()()(23423ssssssssusysG输入下列命令语句执行结果为这个结果表示，该系统的状态空间表达式为uyu]0[1875.0375.04375.01000102000040000161875.07813.0188.210xxx注意，在输入命令中，sys=ss(G)也可以改用[A,B,C,D]=tf2ss(num,den)，在本例中其作用和sys=ss(G)近似，也可以计算出矩阵A、B、C、D。2）.离散系统的状态空间表达式（不要求）离散系统的状态空间表达式为)()()()()()1(kdukCxkykHukGxkx和连续系统状态空间表达式的输入方法相类似，如果要输入离散系统的状态空间表达式，首先需要输入矩阵G、H、C、d，然后输入语句，即可将其输入到MATLAB的workspace中，并且用变量名来表示这个离散系统，其中T为采样时间。如果Gyu表示一个以脉冲传递函数描述的离散系统，也可以用ss(Gyu)命令，将脉冲传递函数模型转换成状态空间表达式。),,,,(TdCHGsssys例2假设某离散系统的脉冲传递函数为47.022.298.323.389.038.057.031.0)(23423zzzzzzzzGyu采样周期为，将其输入到MATLAB的workspace中，并且绘制零、极点分布图。并且将该离散系统脉冲传递函数模型转换成状态空间表达式。sT1.0解输入下列语句语句执行的结果为再输入语句，绘制出零、极点分布图如下在执行完上述语句后，Gyu已经存在于MATLAB的workspace中，这时再执行语句执行结果为结果表示，离散系统的状态空间表达式为)(0001)(05.00000200002235.0555.099.123.3)1(kukkxx)(445.0095.0285.0310)(k.kyx1-2求传递函数矩阵在已知线性定常系统中的A、B、C和D矩阵之后，则该系统的传递函数矩阵可以按下式求出DBAsICsusysG1][)()()(例3已知系统状态方程为21212111013210uuxxxx2121321100003121112uuxxyyy输入以下语句解其中inv()函数是求矩阵的逆矩阵，而simple()函数是对符号运算结果进行简化。执行结果如下这表示ssssAsI213231][211))1/(1()1/(3)2/(1)2/(2)1/(13))1/(3()(sssssssG1-3.线性变换1）化为对角矩阵函数eig()可以计算出矩阵A的特征值以及将A阵转换成对角阵的线性变换矩阵。其语句格式为[Q,D]=eig(A)，则D为对角阵并且对角线上各元素为矩阵A的特征值，满足，因为即：。DAQQ11PQDPAP1例4线性控制系统的状态方程为u1006116100010xx试作线性变换,要求变换后系统矩阵A为对角阵。xQxPx1解先求出系统矩阵的特征值，Q阵可以选择为由特征值构成的范德蒙特矩阵。输入语句可以求出A阵的特征值为－1、－2和－3。因此9413211111PQ输入以下语句执行结果如下由以上计算数据可得系统经过线性变换后的方程为u5.015.0300020001xx也可以输入语句运行结果为再计算线性变换矩阵P，并且验证结果如下可见，两种线性变换虽然不同，却都可以将A阵转换为对角阵2）化为约当矩阵在MATLAB中用函数命令jordan()来求矩阵的约当标准形。其命令格式为：[Q,J]=jordan(A)。输入参量A是系数矩阵，输出参量J是矩阵A的约当标准形矩阵，而就是线性变换矩阵，满足。1QP11PAPQAQJ例5将化为标准形矩阵。452100010A解首先输入语句运行结果为可见，不满秩，即矩阵A的特征值中有重特征值,并且A的独立特征向量的个数小于n。因此输入语句语句执行结果为计算结果表明，矩阵A的约当阵为。100110002J我们验证如下执行结果为所计算出的结果表明，满足1PAPJ二用MATLAB求解系统方程2.1线性齐次状态方程的解使用MATLAB可以方便地求出状态方程的解。我们通过例子来说明。例已知线性系统齐次状态方程为xx3210初始条件01)0(x求系统状态方程的解。解用以下MATLAB程序计算齐次状态方程的解，其中collect()函数的作用是合并同类项，而ilaplace()函数的作用是求取拉普拉斯逆变换，函数det()的作用是求方阵的行列式。程序执行结果这表示ttttttttt2222e2ee2e2eeee2)(ttttt22e2e2ee2)(x2.2线性非齐次状态方程的解通过以下例子说明。例2已知系统状态方程为u103210xx01)0(x)(1)(ttu解用以下MATLAB程序求系统方程的解。其中，语句phi=subs(phi0,’t’,(t-tao))表示将符号变量phi0中的自变量t用(t-tao)代换就构成了符号变量phi，而语句x2=int(F,tao,0,t)表示符号变量F对tao在0到t的积分区间上求积分，运算结果返回到x2。程序执行结果为ttttt22eee5.0e5.0)(x这表示2.3连续系统状态方程的离散化（不要求）在MATLAB中，函数c2d（）的功能就是将连续时间的系统模型转换成离散时间的系统模型。其调用格式为：sysd=c2d(sysc,T,method)。其中，输入参量sysc为连续时间的系统模型；T为采样周期（秒）；method用来指定离散化采用的方法。‘zoh’——采用零阶保持器；‘foh’——采用一阶保持器；‘tustin’——采用双线性逼近方法；‘prewarm’——采用改进的tustin方法；‘matched’——采用SISO系统的零极点匹配方法；当method为缺省时（即：调用格式为sysd=c2d(sysc,T)时），默认的方法是采用零阶保持器。例2-3某线性连续系统的状态方程为BuAxxDuCxy其中6116100010A201201B112011C0000D采用零阶保持器将其离散化，设采样周期为0.1秒。求离散化的状态方程模型。解输入以下语句，其中D=zeros(2)表示，将D赋值为2×2维的全零矩阵。语句执行的结果为计算结果表示系统离散化后的状态方程为)(1936.01164.00902.01959.00047.01099.0)(5112.08366.04429.00738.09541.00246.00041.00984.09991.0)1(kkkuxx)(0000)(112011)(kkkuxy三MATLAB的应用3-1判断线性系统的能控性和能观测性用MATLAB可以很方便地求出线性控制系统的能控性矩阵和能观测性矩阵，并且求出它们的秩。从而判断系统的能控性和能观测性。函数ctrb()和obsv()分别计算系统的能控性矩阵和能观测性矩阵。格式为：Qc=ctrb(A,B),Qo=obsv(A,C)。例判断下面的线性系统是否能控？是否能观测？CxBAxxyu其中103021101A201201B010001C解先分别计算系统的能控性矩阵和能观测性矩阵。然后，再用rank()函数计算这两个矩阵的秩。输入以下语句这些语句的执行结果为从计算结果可以看出，系统能控性矩阵和能观测性矩阵的秩都是3，为满秩，因此该系统是能控的，也是能观测的。注：当系统的模型用sys=ss(A,B,C,D)输入以后，也就是当系统模型用状态空间的形式表示时，我们也可以用Qc=ctrb(sys)，Qo=obsv(sys)的形式求出该系统的能控性矩阵和能观测性矩阵。3-2线性系统按能控性或者能观测性分解在用MATLAB进行结构分解时，不能控（不能观）的系统，其结构分解的系统方程形式有点不同。)ctrbf(][CB,A,KT,,C,B,An)rank(CQ当系统能控性矩阵的秩时，我们可以使用函数命令ctrbf()可以对线性系统进行能控性分解。其调用格式为。其中，T为相似变换矩阵。CCAAAA210CBB0CCCCC输出为一个向量，sum(K)可以求出能控的状态分量的个数。类似地，当系统能观测性矩阵的秩时，我们可以使用函数命令obsvf()可以对线性系统进行能观测性分解。其调用格式为。n)rank(OQ)obsvf(][CB,A,KT,,C,B,A其中，T为相似变换矩阵。O12OAAAA0输出为一个向量，sum(K)可以求出能观测的状态分量的个数。OOBBBOCC0例系统方程为CxBAxxyu其中6101101600A013B100C试按能控性进行结构分解。解输入下列语句语句执行结果为从输出的向量可以看出有两个状态分量是能控的。可以验证，输入语句TTATA得到的结果为可见，A1=Abar，所得到的结果是正确的。3-3线性系统转换成能控标准形和能观标准形下面通过两个例子来说明将系统变换成能控标准形和能观标准形的方法。例系统方程为CxBAxxyu其中