连续系统TransferFuctionModel分析方法程序设计终稿

连续系统TransferFuctionModel分析方法程序设计摘要采用Mathematica软件的传递函数模型分析方法，讨论低阶连续系统的分析方法，并编写了通用的低阶连续系统解析解求解程序，以实例展示程序的用法和程序分析的高效率。关键词连续系统;传输函数模型;拉普拉斯变换;MathematicaContinuousanalysisofprogramdesignsystemTransferFuctionModelAbstract:ThetransferfunctionmodelusingMathematicasoftware,analysis,discussionandanalysismethodsofcontinuouslow-levelsystem,andthepreparationofacommonsystemofcontinuouslow-levelanalyticalsolutionsolver,withexamplesshowingtheuseofprogramanalysisandprogramefficiency.Keywords:continuoussystem;transferfunctionmodel;Laplacetransform;Mathematica目录引言.11用Simulink求解连续系统的思路.11.1用Matlab的Simulink求解连续系统方法.11.2用Matlab的Simulink求解连续系统的步骤.11.3Simulink方法分析求解连续系统的特点.11.4适用的范围.22Mathematica的TransferFunctionModel分析方法求解连续LTI系统的思路.22.1用Mathematica求解LTI连续系统的方法有两种.22.2用Mathematica的TransferFunctionModel分析方法求解连续LTI系统的思路.22.3TransferFunctionModel分析方法的特点.23系统0输入响应计算和0状态响应的计算思路.33.1低阶系统计算思路.33.2高阶系统计算思路.34用Mathematica的TransferFunctionModel分析方法求解连续LTI系统程序设计思路.44.1求解连续LTI系统的设计思路.44.2二阶系统的求解程序.45程序应用实例.55.1程序应用实例一.55.2程序应用实例二.76结语.9致谢.10参考文献.11附录1.12附录2.12第0页共12页引言系统是由若干个相互联系、相互作用的单元组成的具有一定功能的有机整体。LTI系统的分析具有非常重要的意义，一方面是由于在实际应用中经常遇到LTI系统，并且有一些非线性系统或时变系统在限定范围与指定条件下都遵从线性时不变特性的规律；另一方面是LTI系统的分析方法已经形成了严密的、完整的体系，日趋完善和成熟1,2。但是所有的连续系统要进行分析首先要建立微分方程2，所有的离散系统要进行分析需要先建立差分方程，无论是连续系统的微分方程或微分方程组还是离散系统的差分方程或差分方程组要进行求解，都会遇到很多困难。特别是当系统比较复杂时，微分方程的阶数和差分方程的阶数比较高或者微分方程的个数和差分方程的个数比较多的时候，这个方程时域里面求解，是非常具有挑战性的。但是我们可以在Mathematica中用传递函数3或传递函数模型这种方法来求解系统。这种方法比较方便,因为不是在时域里求解。连续系统是在复频域里面进行分析4，离散系统是在Z域里面进行分析，那么在复频域里面或在Z域里面描述系统的方程将转化成简单的代数方程或代数方程组，转换后不会在出现微分方程和差分方程，这就是把人工需要解决的复杂问题转交给了计算机系统来做。计算机在后台来解方程，这样就会很方便的完成任务。所有的电子研发都需要研究系统，都需要求解系统，然而设计的系统合理不合理，对给定的信号有没有加工处理的功能，那就需要我们设计出系统来，给定参数后把解求出来，看一看，是否符合用户的要求，不符合用户的要求就需要改变系统的结构或在系统结构不变的情况下，改变系统的参数的取值，然后再进行求解系统，直到满足用户要求为止。这个研发的过程需要反复进行，所以工作量非常大，手工求解是不太容易的，直接在时域里面编程求解也很困难，所以需要用传递函数或传递函数模型4,5这种方法来求解系统。Mathematica软件，它的功能包括三个方面：数值计算、图形和符号演算；它是一个交互式的计算系统，在于用户相互传递、交换信息数据的过程中完成各种计算；它也是个很容易扩充的系统，除了为用户提供大量的系统函数和数学函数之外，还有一套功能特别强大的程序设计语言，使用这种语言，用户可以方便的定义自己的函数、编制自己的程序，以实现各种功能，完成各种复杂的工作6。它还可以从实际系统制作出来之前，预先对系统进行仿真分析，对系统进行适当的实时修正或者按照仿真的最佳效果来调试及制定控制系统的参数，来提高系统的性能，减少设计系统过程中反复修改的时间，实现效率地开发系统的目标。第1页共12页1用Simulink求解连续系统的思路1.1用Matlab的Simulink求解连续系统方法Matlab的Simulink求解连续系统的方法有三类：（1）离散模块建模的方法；（2）传递函数的方法；（3）状态空间的方法。1.2用Matlab的Simulink求解连续系统的步骤（1）离散模块建模的方法：先写出系统的微分方程，然后把系统的微分方程转化为框图，再把框图用Simulink里面的分立元件模块组装起来；（2）传递函数的方法：把微分方程写出来，直接映射到复频域求出系统的S域系统函数7，把整个系统用一个传递函数模块来表示。这样的建模有三部分，输入信号部分、信号源部分和传递函数模块部分还有数据提取来显示的部分；（3）状态空间的方法：先建立微分方程，再把微分方程转化成状态方程，转化成状态方程以后，它只有四个矩阵，分别是A、B、C、D，然后整个系统就用一个状态空间模块来表示。这个模型中只有三块，输入信号部分、信号源部分、状态空间模块部分和计算的最终结果提取或显示部分。1.3Simulink方法分析求解连续系统的特点用Simulink求解连续系统，解的特点有：（1）所有的解都是数值解，也就是在给定时间点上的采样值，它得不出解析解的形式来；（2）数值解都存在误差，一定范围内，采样周期越大，误差就越大；采样周期越小，误差就越小8，但随着采样周期增大，采样总步数将急剧增加，这时误差将变大，所以说，采样间隔要适当；（3）可以画出时域波形，但随时间变化的演化规律并不能从数据里面直接获取；（4）使用图标编程，对用户来说，求解过程简单。1.4适用的范围因为解是数值解，所以适用的范围是线性系统、非线性系统和混合系统都可以进行求解9，但是这里有一个问题是系统的所有结构必须要给定，系统里面各个部件的参数值也要给定才能进行求解。系统里面不能含有不定参数，这样微分方程、传递函数、状态空间矩阵才能确定，而且里面的参数，例如，电阻、电容、电感这些元件的值都要给定，才能求数值解，因为解的特点限制了它。第2页共12页2Mathematica的TransferFunctionModel分析方法求解连续LTI系统的思路2.1用Mathematica求解LTI连续系统的方法有两种（1）传递函数模型方法先创建传输函数模型函数，myth=TransferFunctionModelHs,s给出系统的输入信号求解零状态响应，yzst_=(OutputResponsemytf,Cos100t,t/Simplify/Expand/TrigReduce/Expand)1创建状态空间函数模型，ss=StateSpaceModelmytf给出系统的初始条件求解零输入响应。yzit_=(OutputResponsemytf,-2,7,t/Expand)1（2）状态空间模型方法先创建传输函数模型函数，mytf=TransferFunctionModel5(s+2)/(s2+7s+12),s;再创建状态空间函数模型函数，ss=StateSpaceModelmytf给出系统的输入信号求解零状态响应，yzst_=(OutputResponsess,Cos100t,t/Expand/N)1给出系统的初始条件求解零输入响应。yzit_=(OutputResponsess,-2,7,t/Expand)12.2用Mathematica的TransferFunctionModel分析方法求解连续LTI系统的思路用传递函数模型的方法求解连续LTI系统，不需要写出微分方程。先用拉普拉斯变换把时域映射到S域，再把系统S域里面的也就是相域里面的方程写出来。如果是电路问题，把电感、电容都用电阻来表示，把原来的动态电路变成相域里面的电阻电路，再把信号源取拉普拉斯变换。在相域里面列的是代数方程，关注哪个元件的电流、电压，就把相电流、相电压与信号源的关系解出来。例如，在代数方程中，就可以用节点电压法、回路电流法求解，求出来后再除以输入信号的相函数，这样系统函数就求出来了。有了系统函数就有了传递函数，解微分方程式可以调用传递函数解方程的方法，这样就可得到解析解（函数解）10。但是这里需要注意的是，这种方法只能解线性系统，非线性系统是不能解的。2.3TransferFunctionModel分析方法的特点用Mathematica的TransferFunctionModel分析方法求解连续LTI系统的出发点是先用拉普拉斯变换把时域映射到S域里面，写出S域里面的方程；优点是可以得到解析函数的表达式，直接得出规律来，并且使用Mathematica的TransferFunctionModel分析方法求解二阶连续LTI系统大大提高了效率，节约了研究人员的时间，减少了人力、物力的使用，省略了繁重的数学计算过程，还可通过图示的方式简单明确的反应问题；难点是对于一阶和二阶系统是很容易求解的，但是对于是高阶系统，它的求解还是非常复杂的，并且它的输出是有错误的，这是因为高阶系统的微分方程中，系数分布范围过广，在后台所对应的一元高次代数方程无解导致的，所以说求解高阶系统还是很难的。第3页共12页3系统0输入响应计算和0状态响应的计算思路3.1低阶系统计算思路研究表明，只能分别求零输入响应和零状态响应，一次把全响应求出来是错误的3。对于简单的系统，我们可以先写出时域里面的微分方程，再取拉普拉斯变换，导出S域里面的系统函数。以一个二阶系统为例可得：（1）先构造一个用户需要的系统，并且把系统映射到S域，映射方法有三条：先把输入信号取拉普拉斯变换；再把电容映射为相域里面容抗，容抗值为1/CS，把电感映射为相域里面感抗，感抗值为LS；例如由微分方程y(t)+a*y(t)+b*y(t)=c*f(t)+d*f(t)可得出S域系统函数Hs_=(c*s+d)/(s2+a*s+b)（2）构造传递函数；mytf=TransferFunctionModelHs,s（3）准备好初始条件；y(0)=p,y(0)=q（4）先调用OutputResponse（）函数，先求零状态响应；yzst_=(OutputResponsemytf,Cos100t,t/Simplify/Expand/TrigReduce/Expand)1（5）再调用OutputResponse（）函数，求零输入响应；ss=StateSpaceModelmytfyzit_=(OutputResponsemytf,p,q,t/Expand)1（6）由零输入和零状态求得全响应；yt_=yzst+yzit（7）最后画系统响应的时域波形图。Plotyzit,t,0,10,AxesLabelt,yzi(t),PlotRangeAllPlotyzst,t,0,2,AxesLabelt,yzs(t)Plotyt,t,0,3,AxesLabelt,y(t),PlotRangeAll3.2高阶系统计算思路例如下面所示的电路系统是四阶的，要写出微分方程就比较困难，所以可以先映射到S域写出在S域里面的系统函数。图(a)时域电路图图(b)S域电路图图3.1时域、S域电路图由图为了求出系统函数，可得到:R4与C2并联的相阻抗为Z1=R4/(1+R4*C2*s)R4和C2并联与L、R3串联的相阻抗为Z2=Z1+R3+Ls在Z2的基础上，再与C1/R2并联的相阻抗为Z3=(Z2+Z2*R2*C1*s)/(1+(Z2*C1+R2*C1)*s)在Z3的基础上，再与R1并联的相阻抗为Z4=R1+Z3由此可得到U(s）=(Z1/(Z3+R1)*Us(s)保存系统结构参数及求系统函数的程序语句如下：R1=200，R2=200，R3=10，R4=800，L=0.001，C1=10-4,C2=2*10-4第4页共12页a3=R4*L*C1*C2;a2=R42*R3*C1*C2+R42*R2*C1*C2+R4*L*C1;a1=R4*R2*C1+R42*C2+R42*C1+R4*R3*C1;a0=R4;b4=R42*R1*C22*C1*L+R42*R2*C22*C1*L;b3=R42*R2*R1*C22*C1+R4*2*R1*C2*C1*L+R42*R3*R1*C22*C1+2*R4*R2*C2*C1*L+R42*R3*R2*C22*C1+R42*C22*L;b2=2*R4*R2*R1*C2*C1+R42*R1*C22+2*R4*R3*R1*C2*C1+R42*R1*C2*C1+2*R4*R3*R2*C2*C1+R42*R2*C2*C1+2*R4*C2*L+R42*R3*C22+R1*C1*L+R2*C1*L;b1=2*R4*R1*C2+2*R4*R3*C2+R42*C2+R1*C1*R2+R3*R1*C1+R4*R1*C1+R3*R2*C1+R4*R2*C1+L;b0=R1+R3+R4;Hs_=(a3*s3+a2*s2+a1*s+a0)/(b4*s4+b3*s3+b2*s2+b1*s+b0)这样就解出了系统函数，接下来就如3.1中（2）至（7）中的步骤进行求解。第5页共12页4用Mathematica的TransferFunctionModel分析方法求解连续LTI系统程序设计思路4.1求解连续LTI系统的设计思路（1）描述系统将描述系统所用到的参数值分别存入多个变量a、b、c、d、p、q中，例如，a=5,b=6,c=2，d=4，p=1，q=0,从而由微分方程y(t)+5y(t)+6y(t)=2f(t)+4f(t),得到系统函数Hs_=2(s+2)/(s2+5s+6);（2）描述系统的输入信号，例如输入信号为：f(t)=cos(100t)（3）利用S域系统函数来创建传输函数模型mytf=TransferFunctionModelHs,s;（4）调用OutputResponse()函数和StateSpace()函数求解当用StateSpace()函数求解有多个解，但需要根据输出方程并且利用状态变量和输出矩阵来构造响应；用OutputResponse()函数就可以直接求出响应，不需要构造。yzst_=(OutputResponsemytf,f(t),t/Simplify/Expand/TrigReduce/Expand)1;ss=StateSpaceModelmytfyzit_=(OutputResponsemytf,p,q,t/Expand)1yt_=yzst+yzit（5）系统响应的表示Printyzi(t)=,yzitPrintyzs(t)=,yzstPrinty(t)=,ytPlotyzit,t,0,10,AxesLabelt,yzi(t),PlotRangeAllPlotyzst,t,0,2,AxesLabelt,yzs(t)Plotyt,t,0,3,AxesLabelt,y(t),PlotRangeAll4.2二阶系统的求解程序例如，要求解微分方程为y(t)+5y(t)+6y(t)=2f(t)+4f(t)的程序为：Clearmytf,sol,y;Hs_=2(s+2)/(s2+5s+6);mytf=TransferFunctionModelHs,s;BodePlotmytf,0.3,500;yzst_=(OutputResponsemytf,Cos100t,t/Simplify/Expand/TrigReduce/Expand)1;ss=StateSpaceModelmytfyzit_=(OutputResponsemytf,1,0,t/Expand)1yt_=yzst+yzitPlotyzit,t,0,10,AxesLabelt,yzi(t),PlotRangeAll;Plotyzst,t,0,2,AxesLabelt,yzs(t);Plotyt,t,0,3,AxesLabelt,y(t),PlotRangeAll;第6页共12页5程序应用实例5.1程序应用实例一例如，要求解微分方程为y(t)+a*y(t)+b*y(t)=c*f(t)+d*f(t)（1）需要先将描述系统的参数值存入变量a、b、c、d、p、q中，a=5,b=6,c=2，d=4，p=1，q=0，即可确定微分方程，从而得到S域的系统函数Hs_=2(s+2)/(s2+5s+6)；（2）给定输入信号f(t)=cos(100t)（3）运行方法使用笔记本电脑，可以把光标放到句子最后面，不要加分号，然后按下Shift和Enter键，系统就会运行。（4）运行结果波特图：图5.1幅频波特图中横坐标为角频率，纵坐标为20lg(|H(i)|)，图5.2相频波特图中横坐标为角频率，纵坐标为ArgH(i)。图5.1幅频频波特图图5.2相频波特图零输入响应：第7页共12页图5.3零输入响应零状态响应：图5.4零状态响应全响应：图5.5全响应第8页共12页零状态响应的解析表达式如下：状态空间矩阵为：零输入响应的解析表达式如下：全响应的解析表达式如下：5.2程序应用实例二例如电路图如下的高阶系统：图5.6电路图(1)电路图中R1=1hm,R2=2hm,R3=2hm,R4=1hm,L=1H,C1=0.1F,C2=0.2F,Us(t)=10Cos(8*t)+5Cos(12*t)+2Cos(40*t)+0.5Cos(200*t)由图为了求出系统函数，可得到:R4与C2并联的相阻抗为Z1=R4/(1+R4*C2*s)R4和C2并联与L、R3串联的相阻抗为Z2=Z1+R3+Ls在Z2的基础上，再与C1/R2并联的相阻抗为Z3=(Z2+Z2*R2*C1*s)/(1+(Z2*C1+R2*C1)*s)在Z3的基础上，再与R1并联的相阻抗为Z4=R1+Z3由此可得到U(s）=(Z1/(Z3+R1)*Us(s)Hs_=(a3*s3+a2*s2+a1*s+a0)/(b4*s4+b3*s3+b2*s2+b1*s+b0)R1,R2,R3,R4,L,C1,C2=200,200,10,800,0.001,10-4,2*10-4;a3=R4*L*C1*C2;a2=R42*R3*C1*C2+R42*R2*C1*C2+R4*L*C1;a1=R4*R2*C1+R42*C2+R42*C1+R4*R3*C1;a0=R4;b4=R42*R1*C22*C1*L+R42*R2*C22*C1*L;b3=R42*R2*R1*C22*C1+R4*2*R1*C2*C1*L+R42*R3*R1*C22*C1+2*R4*R2*C2*C1*L+R42*R3*R2*C22*C1+R42*C22*L;第9页共12页b2=2*R4*R2*R1*C2*C1+R42*R1*C22+2*R4*R3*R1*C2*C1+R42*R1*C2*C1+2*R4*R3*R2*C2*C1+R42*R2*C2*C1+2*R4*C2*L+R42*R3*C22+R1*C1*L+R2*C1*L;b1=2*R4*R1*C2+2*R4*R3*C2+R42*C2+R1*C1*R2+R3*R1*C1+R4*R1*C1+R3*R2*C1+R4*R2*C1+L;b0=R1+R3+R4;（2）给定输入信号ft_=10Cos(8*t)+5Cos(12*t)+2Cos(40*t)+0.5Cos(200*t)（3）运行方法使用笔记本电脑，可以把光标放到句子后面，然后按下Shift和Enter键，系统就会运行。（4）运行结果系统传输函数：图5.7系统传输函数输入信号：图5.8输入信号零状态响应：第10页共12页图5.9零状态响应系统函数的极点和系统函数如下：零状态响应的解析表达式为：输出矩阵值为：零输入响应的解析表达式如下：分析输出结果，里面存在正指数幂，可知这种方法得到的解存在问题。这是因为高阶系统的微分方程中，系数分布范围过广，在后台所对应的一元高次代数方程无解导致的。第11页共12页6结语通过连续系统TransferFuctionModel分析方法程序设计的研究，发现了TransferFuctionModel方法的重要性，运用它解二阶连续系统，是非常方便的，减轻了人工运算的难度，避免了复杂的数学计算，很大程度上，提高了效率。但是对于高阶系统，它还是非常复杂的，并且得不到正确的解。这是因为高阶系统的微分方程中，系数分布范围过广，在后台所对应的一元高次代数方程无解导致的，对于这个问题怎么解决还需要我们进一步研究。由此次研究，让我对传输函数模型有了更深的认识，知道了它的用法及用它给人们带来的益处。使用传递函数模型对时间序列进行分析,最大的意义在于提高预测精度。在利用一个时间序列在某时刻的有效观察值去预报在某个未来时刻该序列的值时,若能够确定出输入时间序列及输入时间序列和输出时间序列之间的动态关系,则可以输入时间序列和输出时间序列的过去值对输出时间序列作出预报11。第12页共12页致谢在导师的指导下，我圆满的完成了这篇毕业论文。在论文的写作过程中遇到了无数的困难，都是在老师和同学的帮助下度过的。特别要感谢我的论文指导老师龙姝明教授，他对我进行了无私的帮助和指导，不厌其烦的帮助我进行论文的修改。除此之外，在校图书馆查阅资料的时候，图书馆里面的老师也给我提供了很大的帮助。在这里我向帮助和指导过我的各位老师表示最衷心的感谢！第13页共12页参考文献1崔本亮.线性时不变(LTI)系统分析方法讨论J.南阳师范学院,2009,2(3):127-127.2刘玉莹,张隽.浅谈连续LTI系统数学模型的几种求解方法J.江西科技师范学院学报,2006,12(6):116-116.3吴大正,杨林耀,张永瑞等.信号与线性系统分析(第四版)M.高等教育出版社,2005:40-50.4李松,张杰,贺国光.城市交叉路口交通流传递函数模型研究J.河北大学,天津大学,2010:236-239.5陈垒.连续系统的传输函数模型辨别J.江南大学硕士学位论文,2012,3:5-5.6孔祥强.Mathematica软件在数学教学中的应用探索J.菏泽学院数学系,2012,6(6):141-141.7张登奇,张璇.连续时间系统的S域分析及实现MATLABJ.湖南理工学院学报,2012,6(2):26-27.8郭美娜.传递函数模型在股市分析中的应用J.华南理工大学硕士学位论文,2012,5(1):13-15.9HONGHua-jie,LIZhi-qiang,FANShi-xun.Multi-ratetrackingcontrolleranalysisanddesignfortargettrackingsystemsJ.JournalofCentralSouthUniversity,2013:3049-3056.10Zhen-Yu-Huo,ZhuYang,Yan-JunPang.Set-point-relatedIndirectIterativeLearningControlforMulti-inputMulti-outputSystemsJ.InternationalJournalofAutomationandcomputing,2012:266-273.11景立伟.传递函数模型及其在医院业务收入中的应用J.山西医科大学,2006,5:1-1.第14页共12页附录1用TransferFuctionModel方法求解二阶连续系统的程序：Clearmytf,sol,y;Hs_=2(s+2)/(s2+5s+6);mytf=TransferFunctionModelHs,s;BodePlotmytf,0.3,500yzst_=(OutputResponsemytf,Cos100t,t/Simplify/Expand/TrigReduce/Expand)1mytf=TransferFunctionModel2(s+2)/(s2+5s+6),s;yzit_=(OutputResponsemytf,1,0,t/Expand)1yt_=yzst+yzitPlotyzit,t,0,10,AxesLabelt,yzi(t),PlotRangeAllPlotyzst,t,0,10,AxesLabelt,yzs(t)Plotyt,t,0,10,AxesLabelt,y(t),PlotPoints100,P