毕业设计（论文）-基于Mathematica函数StateSpaceModel的系统分析方法.doc

陕西理工学院毕业设计论文基于Mathematica函数StateSpaceModel的系统分析方法与程序设计（陕西理工学院物理与电信工程学院电子信息科学与技术专业1103班，陕西汉中723000）指导老师：摘要研究调用Mathematica内部函数StateSpaceModel求解高阶连续时不变系统响应思路与方法步骤。手工和调用Mathematica内部函数NDSolve求解高阶连续LTI系统初值问题都十分困难。我们试图调用Mathematica 内部函数StateSpaceModel及其相关函数，采用建模方法来求连续LTI系统初值问题的解析解或者数值解。研究发现，无论系统阶数高低，用这一方法总可以求出数值解，但是仅当系统阶数n较低时(n4)，才能用这一方法求出系统的解析解。我们编写了采用StateSpaceModel函数思想的高阶连续LTI系统数值解求解实用程序，可以快速求得一般高阶连续LTI系统初值问题的数值解，如果系统阶数较低，也可以求得系统的解析解。关键词 Mathematica;StateSpaceModel;程序设计窗体顶端Systems analysis and programming Mathematica functions based StateSpaceModelLi Lingxi(Grade11,Class3,Major Electronic Information Science and Technology School of Physics and Telecommunication Engineering, Shannxi University of Technology,Hanzhong,723000,Shannxi)Tutor :Long Shuming窗体底端窗体顶端窗体顶端Abstract Research Mathematica call an internal function StateSpaceModel Solving the Higher Order Thinking continuous time-invariant system response and method steps. Manual and call the internal function DSolve Mathematica Solving Initial Value Problems of Higher Order continuous LTI systems are very difficult. We tried to call the internal function StateSpaceModel Mathematica and related functions, using modeling methods to continuous LTI systems analytical solutions or numerical solution of initial value problem. The study found that regardless of the level of order of the system, the total can be obtained by this method numerical solution, but only when the system is low order n (n 4), to the system of analytical solution obtained by this method. We prepared a high value uses StateSpaceModel continuous LTI system function thought to Resolving utility, general numerical solution of initial value problem of higher-order continuous LTI system can be quickly obtained, if the lower order of the system, the system can also be obtained analytical solution.窗体底端Keywords Mathematica, State Space Model, program design窗体底端窗体底端目录1 手工求解连续LTI系统响应的方法和存在的困难12 用Mathematica软件内部函数求连续LTI系统响应的方法43 Matlab的Simulink求解连续系统53.1适用范围53.2连续时间系统建模与仿真分析53.2.1结构图数学模型63.2.2微分方程数学模型83.3解的特点104 Mathematica的StateSpaceModel方法求解连续LTI 系统104.1 创建StateSpaceModel的几种方法104.2 方法步骤105 连续LTI系统初值问题的状态空间模型描述116 连续LTI系统StateSpaceModel解法的程序设计126.1 建立描述系统的微分方程及初值条件126.2 创建状态空间函数模型求零输入响应126.3 创建状态空间函数模型求零状态响应127 StateSpaceModel求解LTI系统的实际应用13结语16附录A18附录B18附录C18附录D18求解系统是科学研究和生产实践当中要解决的很现实的问题，很多的系统设计最后都简化为用一个微分方程来描述。如果我们求出的系统的解不符合用户要求，则必须通过修改系统的结构和修改系统里面元件的参数，重新搭建系统直到设计出来的系统符合用户的需求，而且性价比最好。在这个过程中，我们总是要不断地求解微分方程。但是高阶系统求解，无论是数值解还是解析解手工解都很困难。特别是解析解，高阶系统的解析解手工、编程序都解不出来，但是高阶系统的数值解总是可以求出。用手工计算高阶系统的数值解不可能实现，因为采样点一般都是成千上万的数据。求解系统，科学研究和工程应用当中比较多的要用到数值解。数值解的求解必须编程来计算，因为系统的数据量很大，尽管解法简单，手工解还是太过费时费力。因为我们更多的看中的是系统的数值解，所以数值解我们必须编程来做1，比如可以用NDSolve，Matlab及Simulink里面建模的方法。如今，Mathematica的高版本里面有传递函数、状态空间函数的一些方法，用他们来求解高阶系统的数值解非常容易和简洁2。所以，我们选了这个课题来研究，探索使用StateSpaceModel求解连续LTI系统的响应的方法以及程序设计。1 手工求解连续LTI系统响应的方法和存在的困难手工求LTI连续系统响应(解)，有时域解法、变换域解法。我们仅讨论连续系统，用手工求解系统响应，就是求微分方程初值问题的解。求解可以在时域进行，也可以在S域进行。对于高阶系统，手工求解的难点在于，求一元高次代数方程的根和求多元一次代数方程组的解十分困难3。所以手工只能解低阶系统。下面我们讨论n阶连续LTI时域求解的方法步骤，并展示手工求解的困难所在4。在科学研究和各类电子产品研发调试过程中，我们总要反复求解n阶LTI连续系统初值问题 （1-1）求解这一问题的步骤是：(1)求系统的零输入yzi(t)，它由方程 (1-2)决定。要解这一问题，需要先求一元n次代数方程 (1-3)的n个根，设其由n个单根l1， l2，ln，则得到 (1-4)为了确定其中的n个叠加常数，需要将(1-2)中的初值条件代入(1-4)，得到下面的n元一次代数方程组 (1-5)如果n2，手工求解(1-3)和(1-5)比较困难或者是很费时间。那么，必须考虑在Mathematica,Matlab或Maple环境中编程或调用相关的内部函数来求解才是明智之举。(2)求问题(1-1)的单位冲激响应h(t)求h(t)分两步：第 2 页 共 19 页求h(t)的第一步，先求下面方程决定的h1(t) (1-6)问题(1-6)求解过程与(1-2)过程完全相同，但定解条件不同。分析给出 (1-7)上式中的是单位阶跃信号函数，即自变量大于等于0处，函数值等于1否则函数值为0。(1-7)中的迭加系数满足n元一次代数方程组 (1-8)求h(t)的第二步：解出h1(t)后，由(1-1)方程右端的组合系数分别乘h1(t)的各阶导函数叠加构造h(t)，即 (1-9)(3)由h(t)和输入信号f(t)求系统的零状态响应问题(1-1)的零状态响应为 (1-10)(4)求系统的完全响应 (1-11)例1 求2阶LTI连续系统初值问题解：(1)yzi(t)由方程y(t)+2y(t)+5y(t)=0， y(0_)=5， y (0_)=1给出(2)单位冲激响应h1(t)由方程第 3 页 共 19 页(3)求零状态响应yzs(t)上面的时域卷积，手工计算过于复杂，可采样Laplace变换，各自的像乘积后取逆变换计算yzs(t)，即(4)求系统的全响应这个例子显示，即便是二阶LTI连续系统初值问题求解析解，手工计算也是很难度的。例2 求解三阶LTI连续系统初值问题解：(1) 求零输入响应yzi(t)因初值条件全0，所以yzi(t)=0(2) 求单位冲激响应h(t)由于特征方程手工无法快速求解，调用Mathematica 内部函数Solvel3+2.7l2+5.123l+12.78=0解得l1=-2.61192， l2=-0.0440388-2.21156i， l3=-0.0440388+2.21156i将两个复数特征根重组，得(3) 求零状态响应yzs(t) (4) 全响应y(t)对于这些例子，虽然我们求得了解析解，但我们不得不借助于安装了理工软件包的计算机系统，手工计算一元三次代数方程的根以及简化以s为自变量的有理分式都很困难。由以上求解步骤和实例可见，求解高阶LTI系统解析解必须解决两个初等代数问题:(1) 求一元n次代数方程(1-3)的根，n2时，这个问题手工很难快速求解。(2) 求两个n元一次代数方程组(1-5)和(1-8)，n较大时，手工求解很困难。由此可见，对于n2的高阶LTI系统初值问题， 手工求解析解很困难。2 用Mathematica软件内部函数求连续LTI系统响应的方法Mathematica软件中有多个内部函数可用来求较低阶LTI连续系统初值问题的解析解和高阶LTI连续系统的数值解。调用NDSolve函数，求解连续LTI系统yt+7yt+12yt=200Cos10tExp-0.25t，y0=0，y0=-5数值解需要的程序如下：Cleary;eq=yt+7yt+12yt=200Cos10tExp-0.25t，y0=0，y0=-5;sol=NDSolveeq，y，t，0，10;yt_=yt/.sol1;Plotyt，t，0，10，PlotRange-All，AxesLabel-t，y(t)程序运行后画出系统响应波形曲线见图1。图 1二阶连续LTI系统初值问题响应的衰减震荡型波形曲线这个例子说明NDSolve函数求连续LTI系统数值响应非常方便。我们也可以调用Mathematica状态空间建模思想求解连续系统响应的相关函数TransferfunctionModel， StataSpaceModel，OutputResponse等函数来求解系统y(t)+7yt+12yt=200Cos10tExp-0.25t，y0=0，y(0)=0，y(0)=-5的解析解，需要的程序如下：Cleary，ss，sol;tf=TransferfunctionModel1/(s2+7s+12)，s;ss=StateSpaceModeltfft_=200Cos10tExp-0.25t;yt_=(OutputResponsess，0，-5，ft，t/Chop/Expand)1Plotyt，t，0，10，PlotRange-All，AxesLabel-t，y(t)运行后给出y(t)=11.5753 E(-4. t)-10.1133 E(-3. t)-1.46204 E(-0.25 t) Cos10. t+1.05959 E(-0.25 t)* Sin10. t;其波形曲线见图1。如果要求系统的数值解，将程序修改为Cleary，ss，sol;tf=TransferfunctionModel1/(s2+7s+12)，s;ss=StateSpaceModeltfft_=200Cos10tExp-0.25t;yt_=(OutputResponsess，0，-5， ft， t， 0， 10 /Chop/Expand)1Plotyt，t，0，10，PlotRange-All，AxesLabel-t，y(t)其波形曲线见图1。3 Matlab的Simulink求解连续系统3.1适用范围用Matlab的Simulink求解连续系统的方法主要是用来求线性时不变系统的解。但是因为Simulink方法求解出的是数值解，所以在使用该方法求解连续系统的过程中，系统结构必须要是给定的，而且模型中各个部件的参数值都必须是给定的，不能含有未知数。参数给定后微分方程、系统函数等才能确定，才能使用该方法求解系统。所以，对于系统结果未知或者含有未知参数的系统，该方法不适用。3.2连续时间系统建模与仿真分析通常，线性连续时间系统的数学模型主要有：包含传递函数的结构图和包含状态空间表达式在内的微分方程。相应地，根据数学模型的不同，线性连续时间系统的模型可采用传递函数模块、积分模块或状态方程模块来构建。3.2.1结构图数学模型一个三阶控制系统结构如图2：c(t)r(t)（_+ 0.5图2三阶控制系统结构图运行Matlab2012b，点击home-Simulink Library或在提示符后运行Simulink，进入Simulink环境，点击Simulink library Browser上的file，再选择new-Model进入模型搭建编辑环境。分别从Simulink library中的Sources， Math Operations， Continuous， Sinks各目录中找出单位阶跃信源，加法器，增益器，示波器，输出等模块，并拖入模型文件编辑窗口。设置模块属性。最后点击编辑窗口上的file-Save， 给文件命名xxx.slx并保存模型文件。如图3 所示的系统零输入响应模型。建立系统的Simulink模型：图3系统的Simulink零输入响应模型设置模块属性：设置模块Sum 的模块属性：图4模块Sum 的模块属性图设置模块Zero-Pole 的模块属性：图5模块Zero-Pole 的模块属性图设置模块Transfer fcn的模块属性：图6模块Transfer fcn 的模块属性图设置模块Gain的模块属性：图7模块Gain 的模块属性图仿真结果：图8系统响应波形图3.2.2微分方程数学模型用微分方程描述的数学模型，可利用积分模型直接构建Simulink模型。例子：求解二阶微分方程 x(t)+0.2x(t)+0.4x (t)=0.2u (t)，其中是脉冲信号。根据微分方程构建Simulink模型，如图9：图9系统的Simulink零输入响应模型设置模块属性：设置模块Gain的模块属性：图10模块Gain 的模块属性图设置模块Add的模块属性：图11模块Add 的模块属性图设置模块Gain1的模块属性：图12模块Gain1 的模块属性图设置模块Gain2的模块属性：图13模块Gain2 的模块属性图运行结果：图14系统响应波形图不难看出，Matlab的函数运算及函数绘图功能的强大，使得计算效率大为提高。3.3解的特点Simulink的求解方法所计算出来的系统所有解都为数值解，解出的是给定时间点上面的采样值，是没有办法得到解析函数的形式。Simulink在对连续系统进行求解时，其核心是对系统微分或者偏微分方程进行求解。因此，使用Simulink对连续系统进行求解仿真时所得到的结果均为近似解，只要此近似解在一定的误差范围内便可。该误差是与其系统的采样间隔有关系，采样间隔大的话误差就大，采样间隔小了，误差就比较小。我们可以利用Simulink求解出来的数值解的折线图来得到系统随时间变化的趋势，但是如果我们要得到系统的解随着时间变化的关系或者规律的话，就要进一步进行数据的拟合。4 Mathematica的StateSpaceModel方法求解连续LTI 系统4.1 创建StateSpaceModel的几种方法(1)StateSpaceModela，b，c，d表示状态矩阵为a、输入矩阵为b、输出矩阵为c和转移矩阵为d的状态空间模型。(2)StateSpaceModela，b，c，d，e表示具有描述器矩阵e的描述器状态空间模型。(3)StateSpaceModelsys给出对应于系统模型sys的状态空间模型。(4)StateSpaceModeleqns，x1，x10，.， u1，u10，.，g1，给出通过输出为y、自变量为的常微分或者差分方eqns关于点的泰勒线性化而取得的状态空间模型。窗体顶端4.2 方法步骤求解步骤：(1)用拉普拉斯变换写出系统的象域代数方程；(2)对信号源取拉普拉斯变换；(3)解代数方程，得到信号源与要研究对象的关系式；(4)写出系统函数；(5)求解零输入响应；(6)求解零状态响应；(7)得到完全响应；(8)画出各响应部分的图形。以电路系统为例：(1)首先，根据描述及问题，构造时域电路，画出时域电路图。其次，根据时域电路图映射到S域，画出S域电路图。要求对于线性时不变二端元件R、L、C，若规定其端电压u（t）与电流i（t）为关联参考方向，其相应的象函数分别为U(s)和I（s），那么由拉普拉斯变换的线性性质及微分性质、积分性质可得到它们的S域模型。由时域电路图画s域电路图，只需把时域电路中的电容C改为容抗1/(sC)，电感L改为感抗sL，输入信号f（t）改为象函数f（s）。再次，有序化电路，标清节点5。(2)根据s域电路列写s域节点电压（代数）方程。根据基尔霍夫电压定律、基尔霍夫电流定律及节点电压法，得到自电导乘以节点电压减去互电导乘以“它”节点电压等于流入节点的电流源电流代数和。(3)求解代数方程。(4)根据上一步列出的代数方程，解出输出电压y（s）。系统函数即：系统零状态响应的象函数与激励的象函数f（s）之比称为系统函数，用H（s）表示，即(5)求解零输入响应。(6)求解零状态响应。(7)得到完全响应。(8)画出各响应部分的图形。5 连续LTI系统初值问题的状态空间模型描述对于不同的系统，我们所知道的初值条件是不同的，连续LTI系统的状态空间模型描述方法一共有三种。第一种，微分方程和初使条件已知的情况。由微分方程写出系统函数，对应的也就已知传递函数，根据传递函数模型我们可以得出状态空间模型。求解系统y(t)+7yt+12yt=200Cos10tExp-0.25t，y0=0，y(0)=0，y(0)=-5的解析解，由微分方程写出系统函数H(s)=，将系统函数代入传递函数模型，即tf=TransferfunctionModel1/(s2+7s+12)，s再将传递函数模型代入状态空间模型，即ss=StateSpaceModeltf第二种，像电路系统等系统，微分方程一般写不出，我们可以直接写成系统函数，即系统函数已知的情况6。把系统映射到S域里，直接导出系统函数，再建立传递函数模型，然后把传递函数模型代入状态空间模型中，就建立了状态空间模型。电路系统如图15所示，图15时域电路图我们根据系统电路图15可以画出映射到S域里的象域电路图16。 图16象域电路图需要的程序如下：ClearU2，U3，U4，Us，y;R1=10;R2=100;R3=50;c1=2*10-2;c2=5 10-2;L1=0.5;L2=0.8;eq=(1/R1+c1 s+1/(L1 s)U2-Us/R1-(1/(L1 s)U3=0，(1/(L1 s)+1/R2+1/R3+c2 s)U3-U2/(L1 s)-U4/R3=0，(1/R3+1/(L2 s)U4-U3/R3=0;sol=Solveeq，U2，U3，U4;U2，U3，U4=U2，U3，U4/.sol1;Hs_=Chop(U3-U4)/Us)/Simplify，10-8;tf=TransferfunctionModelHs，sss=StateSpaceModeltf第三种，直接知道描述系统函数的四个矩阵A、B、C、D。把四个矩阵直接写进StateSpaceModel里面，就可以得到状态空间模型7。我们已知一个系统的四个矩阵，如：则可以写出如下程序得到状态空间模型：StateSpaceModel0，1-12，-7，0，1，1，0，06 连续LTI系统StateSpaceModel解法的程序设计6.1 建立描述系统的微分方程及初值条件y(t)+a*y(t)+b*y(t)=c*f(t)+d*f(t)， y(0)=p， y(0)=q给出系统S域系统函数h(s)=(c+d*s)/(s2+a*s+b)6.2 创建状态空间函数模型求零输入响应 Hs_=(d*s+c)/(s2+a*s+b);tfun=TransferfunctionModelHs，s; ss=StateSpaceModeltfun; eqzi=StateResponsess，p，q， 0，t yzit_=eqzi1; 如果方程过于复杂，可以用NDSolve求数值解，再作数据拟合或频谱分析找出解析解8。6.3 创建状态空间函数模型求零状态响应 Hs_=(d*s+c)/(s2+a*s+b); tfun=TransferfunctionModelHs，s; ss=StateSpaceModeltfun; eqzs= OutputResponsess，ft，t yzst_=eqzs1; yt=yzi t+yzst; Plotyt，t，0，tf7 StateSpaceModel求解LTI系统的实际应用在一个实际电路系统中，如图15所示，输入信号us(t)=，电阻R1、R2、R3分别为10、100、50，电容c1、c2为0.02f、0.05f，电感L1、L2为500mh、800mh，求R3两端的输出。图15时域电路图这是一个比较复杂的四阶系统，为了便于分析求解我们要将电感以及电容等效为电阻，画出这个电路的S域等效电路图，如下图16所示：图16象域电路图根据上面的象域电路图，我们可以根据节点电压法9列出如下表达式：这是我们自己的分析过程，计算交给计算机完成。使用Mathematica软件，编程如下：ClearU2，U3，U4，Us，y;t2=4;R1=10;R2=100;R3=50;c1=2*10-2;c2=5 10-2;L1=0.5;L2=0.8;ust_=10Exp-tUnitStept;eq=(1/R1+c1 s+1/(L1 s)U2-Us/R1-(1/(L1 s)U3=0，(1/(L1 s)+1/R2+1/R3+c2 s)U3-U2/(L1 s)-U4/R3=0，(1/R3+1/(L2 s)U4-U3/R3=0;sol=Solveeq，U2，U3，U4;U2，U3，U4=U2，U3，U4/.sol1;Hs_=Chop(U3-U4)/Us)/Simplify，10-8;tf=TransferfunctionModelHs，sss=StateSpaceModeltfyt_=StateResponsess，ust，t，0，t2;(*Printy(t)=，yt*)Plotyt1，t，0，t2，PlotRange-All，AxesLabel-t，y1(t)Plotyt2，t，0，t2，PlotRange-All，AxesLabel-t，y2(t)Plotyt3，t，0，t2，PlotRange-All，AxesLabel-t，y3(t)Plotyt4，t，0，t2，PlotRange-All，AxesLabel-t，y4(t)yzit_=OutputResponsess，0.0001，-0.001，0.02，-0.002，0，t，0，t2;(*Printyzi(t)=，yzit*)Plotyzit，t，0，t2，PlotRange-All，AxesLabel-t，yzi(t)yzst_=OutputResponsess，0，0，0，0，ust，t，0，t21/Chop/Expand;(*Printyzs(t)=，yzst*)Plotyzst，t，0，t2，PlotRange-All，AxesLabel-t，yzs(t)y0t_=yzit+yzst;(*Printy0t_=，y0t_*)Ploty0t，t，0，t2，PlotRange-All，AxesLabel-t，y0(t)程序运行后，传递函数模型为：状态空间函数模型为：系统响应波形曲线见图1723。图17状态变量1波形曲线图图18状态变量2波形曲线图图19状态变量3波形曲线图图20状态变量4波形曲线图图21系统零输入响应波形图22系统零状态响应波形图23系统全响应波形结语本次毕业设计，我如期完成了任务书中的内容，设计出了可以求解高阶连续LTI系统数值解的通用程序。该程序通过调用Mathematica软件中的StateSpaceModel函数构造系统的状态空间模型来求解连续LTI系统，并对系统进行分析。解决了人工计算求解复杂LTI系统数值解的棘手问题，节省了大量的人力物力，为研究连续LTI系统提供了简便快捷的方法。致谢毕业设计的完成，也就意味着自己给大学四年交了一份试卷，至于合格不合格也就在于论文的内容精彩与否。本次毕业设计从三月份选题开始，一切都进行得有条不紊，按照学校的安排以及龙老师的安排，我们小组的同学都积极有效的完成了一系列工作。挑选课题，下载任务书，完成开题报告，接着就到了最重要的环节，论文初稿的完成。我大约用了一个月的时间来做。首先最重要的就是温故所学知识，包括信号与系统，数学物理方法中的相关内容，比如系统、信号的基本概念，连续系统响应的解法，连续系统的响应等。学习了基本知识，还要回顾软件的使用，我所需要使用的软件有Matlab、Mathematica，虽然我是基于Mathematica的分析，但是Matlab这一重要软件的学习也是不可小觑的。在做毕业设计的过程中，除了自己的努力之外，我还得到了同学和龙老师的鼎力相助。他们都很热心的帮我查资料，给我讲解部分知识的难点，对于这次毕设，给了我很大的帮助，我很感谢他们和龙老师。特别是龙老师，经常抽出自己的休息时间给我讲解我不懂的知识，每次我都受益匪浅。还有图书馆的管理员，他们热心的帮我查找资料。现在，毕业设计已经基本完成，在这一过程中，我经历了很多，也学到了很多，凡是都不容易，我们要刻苦努力、乐观向上，这样才会收获自己想要的结果。参考文献1李汉龙 ，缪淑贤 ，韩婷.Mathematica基础及其在数学建模中的应用M 国防工业出版社，2013：.2孔祥强.Mathematica软件在数学教学中的应用探索J计算机应用与软件，2012，29（6）:141-144.3王福贵，侯建文.基于Mathematica的线性方程组的可视化求解J计算机与现代化，2013，5:149-154.4章美月.数学软件Mathematica及其应用M中国矿业大学出版社，2010:34-355熊青玲，潘林峰，程衍富.Mathematica软件在示波器实验中的应用J实验科学与技术，2013，11（2）：37-39.6杜海莲，王占峰，刘波粒，付超.Matlab在电路分析实验教学中的应用J江苏技术师范学院学报，2011，17(10):89-947龙飞，刘娇.Matlab在LTI连续系统响应分析中的应用J四川理工学院学报，2006，19（3）：101-104.8Tin Win Mon， Myo Myint Aung.JWorld Academy of Science， Engineering and Technology，2008，2:113-118.9Romke Bontekoe，MBontekoe Research， Amsterdam romkebontekoe.nl第 20 页 共 19 页附录ACleary;eq=yt+7yt+12yt=200Cos10tExp-0.25t，y0=0，y0=-5;sol=NDSolveeq，y，t，0，10;yt_=yt/.sol1;Plotyt，t，0，10，PlotR