连续系统simulink传递函数建模分析方法程序设计初稿

连续系统simulink传递函数建模分析方法与程序设计【摘要】：建立LTI连续系统的时域微分方程和复频域(S域)系统函数的数学模型、Simulink模型，编写求解连续系统数值解即可视化MATLAB程序。选择典型的RLC三阶电路系统，运用所建立的仿真模型和程序求解电路的零输入、零状态和全响应，并对其进行理论分析研究。【关键字】：连续系统；传递函数；simulink模型；MATLAB程序设计Themodelinganalysismathsandprogramdesignofthecontinuoussystemssimulinktransferfunction【Abstract】：EstablishingthemathematicmodelingofLTIcontinuoussystemstimedomaindifferentialequationsandcomplexfrequencydomain(S-domain)systemandthemodelingofSimulink,PreparedtosolvethenumericalsolutionofacontinuoussystemcanvisualizetheMATLABprogram.ChooseatypicalthirdorderRLCcircuitsystem,Theuseofsimulationmodulesandcircuitsarecreatedtosolvezeroinput,zerostateandtheentireresponse,analysizingandresearchingitstheory.【Keywords】：continuoussystem;transferfunction;simulinkmodel;MATLABprogamming目录引言.21探索快速描述LTI连续系统的传递函数方法及频谱特性.21.1理论分析方法.21.2LTI连续系统时域映射到复频域(S域)的必要性及S域的方法和思路.21.3LTI连续系统时域映射复频域的方法及所需响应.31.4定义系统函数.42.探索利用simulink传递函数仿真模型文件求解连续线性系统数值解的思路、技巧和方法.52.1连续系统电路的时域映射到复频域（S域）.52.2创建simulink传递函数和状态空间仿真模型文件.63.连续系统数值解及其可视化Matlab程序.73.1程序设计的思路与技巧.73.2程序语句.74方案例证.84.1S域电路系统函数和连续系统数值解的思路.84.2simulink仿真模型及程序求解电路响应.114.2.1传递函数模型：.114.2.2状态空间模型：.134.3程序运用实例展示.165结束语.20致谢.20参考文献.21第0页共22页引言现代社会的发展，诸多领域应用的系统都是连续系统，如科研、生产实践、产品和仪器检测等。其连续系统中的控制电路都应该满足在一定的频带范围内，具有一定的放大和延迟功能。连续系统的解法有解析解和数值解两种，相比而言，连续系统的时域解析解法虽然便于理论分析系统响应的变化趋势和系统特性，但实际系统总是多输入多输出的高阶系统，它们的解析微分方程书写困难，时域响应求解极为困难，出错率也较高，即便较低阶系统的解析方程能够得到，其求解也较复杂，耗时耗力。而连续系统的区间数值解法就比较容易，本质上用的是迭代解法，总是能够方便、快速地的得到，之后如果企图观察其响应随时间演化的趋势，可用数值解画出其波形来观察，甚至必要时做数据拟合寻找区间解的拟合函数也是人可能的，而且数值解法还可以求一定区间上的非线性问题。将数值解和传递函数、状态空间simulink仿真建模及MATLAB程序结合，解出LTI连续系统高阶电路的零输入、零状态和全响应，并画出各个状态的时域响应及频谱波形，以使我们更好的研究信号随时间或频率的演化规律。第1页共22页1探索快速描述LTI连续系统的传递函数方法及频谱特性1.1理论分析方法LTI连续系统分析方法有时域分析和变换域分析：时域分析：对于给定的激励，根据描述系统响应与激励之间关系的微分方程求得其响应的方法，一般情况下用卷积积分。我们借助于连续时间LTI系统的单位冲激响应表示如下：*ytxhtdxth例：一个二阶电路微分方程为：，信号源为56ftfe设初值，求解系统的传递函数。01,x(1)先解特征方程：2560得：23(2)求h(t)再利用拉普拉斯逆变换写出写出系统函数：Hs2156若求解系统零输入响应就利用经典解求；零状态响应为所求激励与激励响应的卷积积分；全响应为零输入响应和零状态响应之和。*zsxhtf变换域分析：利用拉普拉斯变换将时域映射到S域，然后根据代数方程输出比输入得系统函数，将系统函数与仿真建模、程序结合就可以求出系统的零输入响应、零状态响应和全响应的波形。再令，程序中调用fft函数对系统进行频谱分析画出相频图和幅频图。sj1.2LTI连续系统时域映射到复频域(S域)的必要性及S域的方法和思路线性时不变（LTI）连续系统的时域分析方法，即对于给定的激励，根据描述系统的响应与系统之间关系的微分方程求得其响应的方法，其主要方法为经典解。但是利用经典解在求解微分方程的基础上讨论其零输入、零状态和全响应比较复杂，高阶解更是困难。在应用傅里叶积分变换求解LTI连续系统微分方程中，有许多重要函数不满足傅里叶积分变换的绝对可积条件，例如常数、单位阶跃函数以及正余弦函数等。这使得傅里叶变换求解微分方程有很大的限制，而且傅里叶积分变换比较烦琐。在利用傅里叶分析研究LTI系统时，将只局限于系统的冲激响应有傅里叶变换的情况。在LTI连续系统的分析和研究中，拉普拉斯变换是一种特别有用的分析工具，它将描述系统的时域微积分方程变换为S域的代数方程，便于运算和求解；同时它将系统的初始状态自然地包含于象函数方程中，即可分别求得零输入响应、零状态响应和系统的全响应。综上，从时域、傅里叶变换和拉普拉斯变换三种方法解微分方程中得出，拉普拉斯变换是最方便易行的。拉普拉斯变换的重要应用之一是对于LTI系统的分析与表征。对于LTI系统，利用积分变换给时域变量求拉普拉斯变换。拉普拉斯变换的作用直接来源0stXexd于卷积性质x1(t)*x2(t)X1(s)X2(s),根据这一性质就可以得到，一个LTI系统输LS入和输出的拉普拉斯变换是通过乘以单位冲激响应的拉普拉斯变换联系起来的，即Y（s）=H（s）X（s），X（s）、Y（s）和H（s）分别是系统输入、输出和单位冲击响应的230,tthtthtte第2页共22页拉普拉斯变换。传递函数为：，所以S域求解系统函数很容易。YHSF初值设置：若微分方程为n阶，则n-1，n-2,n-32,1阶的初始状态可设初值。1.3LTI连续系统时域映射复频域的方法及所需响应LTI连续系统的数学建模是常系数微分方程，用拉普拉斯变换求解微分方程：设LTI系统的f(t),响应y（t），描述n阶系统的微分方程的一般形式可写为（1.3-1）()()00nmijijttfab式中，系数均为实数，设系统的初始状态为,(,1)ija、（=,1）0y,。(1)0y)n令，。根据时域微分定理，y（t）及其各阶导数的拉普拉(tYsL()ftFs斯变换为（）（1.3-1()()0()iiipytsYsyi=0,1n2）如果f(t)是t=0时接入的，则在时f(t)及其各阶导数均为零，即t。因而f(t)及其各阶导数的拉普拉斯变换为()0(,1)jfjm(1.3-3)()jjftsFL取式（1.3-1）的拉普拉斯变换并将式（1.3-2）、式（1.3-3）代入得1()000()()nimipjijasYsybsF即（1.3-1()000()()nnimipjiijssys4）由上式可解得（1.3-()()MsBYFsA5）式中，是方程（1.3-1）的特征多项式；，多项式和0()niiAsa0()mjjsb()As的系数仅与微分方程的系数、有关；，它也是()Biajb1()0()nipiMay第3页共22页s的多项式，其系数与和响应的各初始状态有关而与激励无关。ia()0py由式（1.3-5）可以看出，其第一项仅与初始状态有关而与输入无关，因而是零输入响应的象函数，记为；其第二项仅与激励有关而与初始状态无关，因而是零状态响()ziyt()ziYs应的象函数，记为。于是式（1.3-5）可写为zszs（1.3-()()()zizsMBsFA6）式中，。取上式逆变换，得系统的全响应()ziMsYA()zsBYFs（1.3-()()zizsyttyt7）1.4定义系统函数描述n阶LTI系统的微分方程一般可写为：（1.4-()()00nmijijaytbft8）设是时接入的，则其零状态响应的象函数为()ft（1.4-()zsBYFsA9）式中为激励的象函数，、分别为()Fs()ft()s（1.4-0()niimjjAaBsb10）它们很容易根据微分方程写出。系统零状态响应的象函数与激励的象函数之比称为系统函数，用表示，()zsY()Fs()Hs即（1.4-()()zsBHFA11）第4页共22页由描述系统的微分方程容易写出该系统的系统函数，反之亦然。由式（1.4-()Hs11）以及式（1.4-10）可见，系统函数只与描述系统的微分方程系数、有关，即只()Hsiajb与系统的结构、元件参数等有关，而与外界因素（激励、初始状态等）无关。第5页共22页2.探索利用simulink传递函数仿真模型文件求解连续线性系统数值解的思路、技巧和方法2.1连续系统电路的时域映射到复频域（S域）首先列出三种元件（R、L、C）的时域和S域的关系如表2-1：表2-1电路元件的S域模型电阻电感电容基本关系utRit1tiditutL001Ltxii010tcutixduCti串联形式UsRI0LsIi01CuUsIsS域模型并联形式1sURI01LisUsI0CsUsuI搭建时域电路系统：图2.2时域二阶RLC串联电路根据上述表中的S域模型,将时域电路变成S域电路并进行分析求解：sL第6页共22页图2.3S域电路模型根据KVL定律的复频域形式及元件电流关系的复频域形式，可求得代数方程为：010cLuRIsIiIsUsC代入电感与电容的初值代数方程为：IsIIs系统函数为:1()HRsLC2.2创建simulink传递函数和状态空间仿真模型文件运行Matlab2012b,点击Home-SimulinkLibrary或在提示符后运行simulink，进入simulink环境，点击simulinklibraryBrowser上的File，再选择New-Model进入模型搭建编辑环境。分别从simulinklibrary中的Sources,Continuous,Sinks各目录中找出脉冲方波信源，传递函数，状态空间模型，输出等模块，并拖入模型文件编辑窗口。用信号流线联接各模型建立仿真模型。点击编辑窗File-Save,给文件命名xxx.mdl并保存模型文件。需要输出到内存的数据都需要将其送入一个专门的Out模块，所有输出数据按Out口号被组织到输出数组变量中。设置系统结构各模块参数:a.信号源的设置信号源的脉冲类型、幅度、周期、高低电平和相位根据具体要求具体设置；因为传递函数不能求解出连续系统的零输入响应，而状态空间模型可以设置系统初始条件，所以把传递函数转换成状态空间模型求解LTI电路系统的零输入响应，所以状态空间根据所求响应设置幅度。b.传递函数和状态空间模型的设置打开传递函数模块，根据系统函数分子分母多项式中变量S各幂次的系数，幂次由高到低排列来设置模块中的分子系数和分母系数，一般把分子分母系数定义成符号变量，传递函数不能设置初始值所以只能求零状态响应；打开状态空间模型，运用程序中系统函数转换成状态空间的矩阵系数，设置所得矩阵系数为符号变量到模块中，然后设置初值也为符号变量。运行仿真前，要在Matlab命令窗口提示符后或程序中为符号变量赋数值。点击模型文件窗口的Simulation-Modelconfigurationparameters，选择solver项，设置solver中的：仿真时间段tf，采样周期Ts，步长类型(定步长Fix-step)，解算器类型(ode型)。点击模型文件窗口的Simulation-Modelconfigurationparameters，选择Datainput/output项，设置数据选项：在Savetoworkspace下面的Time,States,Output框上点击打钩，并将后面的变量名改为简单变量名t,x,y,系统必定输出状态变量数据到程序所定义的内存数组中。Format选为Array。取消Limitdatapointstolast1000前面的钩，取消对话框下面各处的钩。取消显示模块中Limitdatatolast1000前的钩。每当模型文件发生变化时，都是应该及时保存其变化。模型文件扩展名保存为mdl。第7页共22页3.连续系统数值解及其可视化Matlab程序3.1程序设计的思路与技巧设计程序前，以防产生不必要的干扰，首先清理磁盘和界面，关闭所有窗口，然后设计程序。根据所知的系统函数，其函数分子分母系数用符号变量代替，把符号变量按设置要求设置到模块中，因为有时连续系统系数比较复杂，不能总是方便的设置到仿真模块中，所以要在程序中设置个符号变量代替所有系数。这样，当函数系数变化时，我们就不用重新打开仿真模块另行设置，只要在程序中变换就可以了。首先对系统函数作频谱分析，得知输入信号的通频带，根据通频带设置设置采样时间，这样能达到预期的效果。电子元件的量值在程序中直接定义；设置仿真时间段和采样周期，采样的周期要在频带范围内，且要满足周期的个数为整数；设置时间变量t，状态变量x和输出变量y，运行程序要调用模型文件的语句t,x,y=sim(模型文件主名,0,tf)。因为要求零输入响应、零状态响应，所以模型文件名不同，变量的定义要标明，例如变量可以定义分t1,x1,y1依次类推；输出的变量要确保是个值，所以计算时要把有的变量转置，然后存在二维数组中。画连续系统波形图要调用plot，画出来为连续的折线图；调用fft函数画频谱图，频谱分析画图要为杆状图是离散的要用stem函数画图，这样更能分析出系统随频率变化的特性。这种程序的设计更体现了其可读性，调试方便性和良好的通用性；运行程序时，点击simulink模型文件编辑窗口上的“三角”，运行模型文件对应程序代码。3.2程序语句a.t,x,y=sim(模型文件名,0,tf)在模型文件函数中定义时间变量t、状态变量x、输出变量y；tf为截止时间，0,tf为仿真时区。sim函数括号里必须为时间段。b.系统函数分子系数num=，分母系数den=，中系数用空格隔开;A,B,C,D=tf2ss（num,den）调用函数tf2ss，把传递函数转换成状态空间，A、B、C、D为状态空间矩阵系数；c.W=(0:采样时间:截止时间)*2*pi;w为角频率，在截止时间范围内对信号频谱采样；d.H=reshape(freqresp(mytf,W),1,length(W);在频带范围内，对矩阵w进行重塑得H；e.T=4;Ts=T/100;fs=1/Ts;tf=5*T-Ts;设置输入信号周期，采样周期、采样频率、采样时区范围；f.信号1=(2/时间长度)*fft(信号);对信号进行fft变换；g信号=abs(傅里叶变换的信号(1:Ln2);频谱分析对信号取模值；hph1=angle(Y(cr1);对经fft变换的信号Y求其相位存入变量ph1;i.subplot(所画图总行数、每行列数、本图所处行数);plot(折线图)/stem（杆状图）(时间,所求变量,-（线图）/.(点图);xlabel(所求变量横轴单位);ylabel(所求信号:定义名);第8页共22页4方案例证电阻、电感及电容是电路的基本元件，在交流电或电子技术中，常需要利用电阻、电感及电容元件组成不同的电路，用来改变输入正弦信号和输出正弦信号之间的相位差，可以构成各种振荡、选频电路、滤波器等。具有电阻电感电容的无源二端网络RLC三阶时域电路，其中：R=0.4，L=1H,C=2F。现已知电压为输入，电压为输出。1()ut2()ut图4.1RLC时域电路时域电路映射到S域：图4.2RLC复频域电路4.1S域电路系统函数和连续系统数值解的思路a.若电容电感初始值为零，则S域电路为图4.3所示求系统函数图4.3RLC复频域电路如图4.3所示，由复频域电路图建立复频域代数方程:11OCsUUsLRC其等效阻抗：21第9页共22页01()2sLRCZs输出象函数：20()2OCRUsUsZ12OCssLR1122sUssLRLRCC最后的输出象函数整理得：213211212LUsUsRLRCLRssC由输出象函数得系统函数：32112121LHsRLRCLRssC代入电阻、电感、电容的值得：32406958Hssb.若电容电压初始值不为零，则根据系统函数列时域微分方程，利用拉普拉斯变换，代入所给的初始值求解其连续系统的零输入响应、零状态响应和全响应。设输入为，初始值,，1()ftut(0)cyup(0)Lyiq(0)Lyir根据传递函数列出时域电路微分方程：121212121RRCLRRtytytytftCL第10页共22页令，12169RLa12258RCLb，p=2,q=0,r=0；40dC40kL把系数代入微分方程得：abd、ytatbytdtkft给方程两边同取拉普拉斯变换得：3220000sYyssYybsYydYskFabdyabak即得全响应的象函数为：23200ssykFsYabd23232abyssssabsd由上式得复频域的零输入响应的象函数为：23200ziyayYssbsd零状态响应的象函数为：32zskFYas代入系数值和初始值得零输入象函数为：2318056694ziss零状态象函数为：32580zsFYs因为给定的输入为周期脉冲，脉冲宽度,所以其象函数为T21Tsse则零状态象函数为：第11页共22页232140695840TszseYs对上式取逆变换得零状态响应：-=zizszizszizsytYytt111LL+4.2simulink仿真模型及程序求解电路响应4.2.1传递函数模型：传递函数模型中的传递模块不能设初始值所以只能求解零状态响应首先创建仿真文件建立仿真模型,文件命名为yf301.mdl并保存模型文件信号源的脉冲类型选择为Timebased,幅度为1，周期为T，程序中给出周期值，高低电平各占50%，相位0。第12页共22页打开传递函数模块，设置参数：分子系数k和分母系数1abd分别是分子、分母多项式中变量s各幂次的系数，幂次由高到低的顺序排列。运行仿真前，要在Matlab命令窗口提示符后或程序中为符号变量a,b,d,k赋数值。点击模型文件窗口的Simulation-Modelconfigurationparameters，选择solver项，设置solver中的：仿真时间段tf，采样周期Ts步长类型(定步长Fix-step)，解算器类型(ode型)。第13页共22页点击模型文件窗口的Simulation-Modelconfigurationparameters，选择Datainput/output项，设置数据选项：在Savetoworkspace下面的Time,States,Output框上点击打钩，并将后面的变量名改为简单变量名t,x,y,系统必定输出状态变量数据到内存数组yzs中。Format选为Array。取消Limitdatapointstolast1000前面的钩，取消对话框下面各处的钩。取消显示模块中Limitdatatolast1000前的钩。每当模型文件发生变化时，都是应该及时保存其变化，模型文件扩展名有mdl。4.2.2状态空间模型：用状态空间模型求解系统零输入第14页共22页创建仿真文件建立仿真模型，文件命名为yf302.mdl并保存模型文件求连续系统的零输入响应时，信号源的脉冲类型选择为Timebased,幅度为0，周期为T，高低电平各占50%，相位为0。打开状态空间模块，设置参数：矩阵系数为A、B、C、D，初始状态设置为符号变量p、q、r，运行仿真前，要在Matlab命令窗口提示符后或程序中为符号变量第15页共22页a、b、d、k、p、q、r赋数值。点击模型文件窗口的Simulation-Modelconfigurationparameters，选择solver项，设置solver中的：仿真时间段tf，采样周期Ts步长类型(定步长Fix-step)，解算器类型(ode型)。点击模型文件窗口的Simulation-Modelconfigurationparameters，选择Datainput/output项，设置数据选项：在Savetoworkspace下面的Time,States,Output框上点击打钩，并将后面的变量名改为t2,x2,y2,系统的输出变量数据存到yzi数组中,系统第16页共22页的输出状态变量数据存到内存数组中。Format选为Array。取消Limitdatapointstolast1000前面的钩，取消对话框下面各处的钩。取消显示模块中Limitdatatolast1000前的钩。每当模型文件发生变化时，都是应该及时保存其变化。执行程序时要调用模型文件，模型文件的语句t2,x2,y2=sim(模型文件主名,0,tf),点击simulink模型文件编辑窗口上的“三角”，运行模型文件对应程序代码，求解系统响应-零输入响应。4.3程序运用实例展示程序如下所示：clc;clearall;closeall;R1=4;L1=0.5;c1=0.05;L2=1;R2=60;a=(L1*R2+L2*R1+L1)/(L1*L2);b=(R1*c1+L2+c1*R1*R2)/(c1*L1*L2);d=(1+R2)/(L1*L2*c1);k=R2/(L1*L2*c1);num=k;den=1abd;mytf=tf(num,den);W=(0:0.01:5)*2*pi;H=reshape(freqresp(mytf,W),1,length(W);plot(W,abs(H);xlabel(w/(Rad/sec);ylabel(:H(w);第17页共22页T=4;Ts=T/100;fs=1/Ts;tf=5*T-Ts;t,x,y=sim(yf301,0,tf);figure;in=y(:,1);yzs=y(:,2);Ln=length(y);Ln2=floor(Ln/2);IN=(2/Ln)*fft(in);IN(1)=IN(1)/2;Y=(2/Ln)*fft(yzs);Y(1)=Y(1)/2;AY=abs(Y(1:Ln2);F=(0:(Ln2-1)*fs/Ln;AIN=abs(IN(1:Ln2);subplot(511);plot(t,in,.);xlabel(t/sec);ylabel(:in(t);subplot(512);stem(F,AIN,.);xlabel(f/Hz);ylabel(:A(f);subplot(513);plot(t,yzs,-);xlabel(t/sec);ylabel(:yzs(t);subplot(514);stem(F,AY,.);xlabel(f/Hz);ylabel(:A(f);A,B,C,D=tf2ss(num,den);p=2;q=0;r=0;x0=p;q;r;t2,x2,y2=sim(yf302,0,tf);subplot(515);plot(t2,y2,g-);xlabel(t);ylabel(yzi(t)Am=max(AY);th=0.01;cr1=find(AY=th*Am);AY1=AY(cr1);F1=F(cr1);sp=F1,AY1;ph1=angle(Y(cr1);yzs2=yzs;forj=1:Lnyzs2(j)=cos(2*pi*t(j)*F1+ph1)*AY1;endfiguresubplot(311);plot(t,in,.);xlabel(t/sec);ylabel(:in(t);subplot(312);plot(t,yzs,r-);xlabel(t/sec);ylabel(:yzs(t);subplot(313