系统辨识的Matlab实现方法(手把手)

最近在做一个项目的方案设计，应各位老总的要求，只有系统框图和器件选型可不行，为了凸显方案设计的高大上，必须上理论分析，炫一下“技术富”，至于具体有多大实际指导意义，那就不得而知了！本人也是网上一顿百度，再加几日探索，现在对用matlab实现系统辨识有了一些初步的浅薄的经验，在此略做一小节。必须要指出的是，本文研究对象是经典控制论理最简单最常用的线性时不变的siso系统，而且是2阶的哦，也就是具有如下形式的传递函数：本文要做的就是，对于有这样传递函数的一个系统，要辨识得到其中的未知数T，£!!这可是控制系统设计分析的基础哦，没有系统模型，啥理论、算法都是白扯，在实际工程中非常重要哦！经过总结研究，在得到系统阶跃响应实验数据之后（当然如果是其他响应，也有办法可以辨识，在此还是只讨论最简单的阶跃响应实验曲线，谁让你我是菜鸟呢），利用matlab至少可以有两种方法实现实现（目前我只会两种，呵呵）！一、 函数法二、 GUI系统辨识工具箱下面分别作详细介绍！函数法看官别着急，先来做一段分析（请看下面两排红*之间部分），这段分析是网上找来的，看看活跃一下脑细胞吧，如果不研读一下，对于下面matlab程序，恐怕真的就是一头雾水咯！G（s）可以分解为:1G(s)= ”T2(s+o)(s+o)121[—]3广T昇占冃其中，3=1 ]2T3、13、1则有都是实数且均大于零。3+31223312传递函数进一步化为：因此,辨识传递函数就转化为求解312当输入为单位阶跃函数时,对上式进行拉普拉斯反变换,得系统时域下的单位阶跃响应为：33①一①21即1—y(t)= 2 e-31t— 1 e-叱t①一①213-321令3=k3 (k〉1),得12对上式两边取以e为底的对数得当t时,当t时,ln1一e一(k一1)32tk则上式化简为该式的形式满足直线方程其中，y*y*(t)」n(1-y(t)]a=-3k2,b=lnE(k>1)通过最小二乘算法实现直线的拟合，得到a,b的值，即可得到3、3的值，进而12可得系统的传递函数。Matlab程序代码:clccloseallt=[135791113151719];y=[0.1490860.58900670.8306170.9339900.9739800.9910950.9958680.9986800.9994900.999850];y2=log(1-y);plot(t,y2,'*');gridonpm=polyfit(t,y2,1)value=polyval(pm,t);holdonplot(t,value,'r')title('\fontname{黑体}\fontsize{20}y(t)二at+b')w2=-pm(1)w1=w2/(1-exp(-pm(2)))T=1/sqrt(w1*w2)theta=(w1+w2)/(2*sqrt(w1*w2))z=[];p=[-w1-w2];k=w1*w2;sys=zpk(z,p,k)figure(2)step(sys,[0:0.5:20]);axis([02001.2])holdonplot(t,y,'r*')打开matlab，新建一个Function，把上述程序段拷进去，保存，运行~~~~~~~~~，运行结果：系统的传递函数为很顺利吧？先高兴一个！问题接着马上就来了，上面这个例子，这个传递函数的极点刚好都是负实数，因此辨识得很顺利，但是如果系统是欠阻尼系统，也就是如果传递函数的根是复数，那么上述函数段，就无能为力咯，会出现说“matlab无法处理增益为复数情况之类 ”例如对于这个系统，若果用simulink做一下阶跃响应，再把实验数据代入上述函数段，那就不行咯！怎么办呢，只能另辟蹊径了！二、（SystemIdentificationTool）系统辨识工具箱早听说matlab博大精深，神通广大了，于是乎我确定肯定有更简单、直观、强大的工具来完成这小儿科把戏。查资料琢磨之后，我做了个小实验，在simulink里验证了该种方法。该方法的大原则是：在确定了系统的输入输出数据两个列向量NX1形式，如果是1XN,会提示出错）之后，设计好一定的辨识原则（比如说是2阶？3阶？，传递函数是零极点形式，还是带阻尼形式，等等），然后就交给强大的matlab，得到辨识结果。Stepbystep,plz!Step1、建立模型获取系统输入输出数据图1图1系统的输入是阶跃信号，用Scope1监视，并输出到workspace（这步不会的自己百度哦），采样周期是0.1s，得到输入变量u（101X1的矩阵）；本人在系统的阶跃响应上叠加了一白噪声，当然也可以不加噪声，加了噪声就是期望更真实的模拟实际情况，白噪声参数设置见图2图2同样在Scope2监视，也将结果输出到workspace,得到响应数据y（同样也是101X1的矩阵）Step2、进入辨识工具箱&设置辨识规则直接在commandwindow输入ident,回车，进入辨识工具箱图3图3点击import下拉菜单，选时域数据timedomaindata，见图4图4在下图5红色圈区域输入之前得到的系统输入和输出数据，u和y图5在下图6绿色圈内输入数据的一些信息，因为之前模型中，阶跃起点我是放在0s处的，这里也设置0,如果前面模型仿真是1s,这里应该也是1s;采样时间是0.1s,根据实际情况设置统一哦图6设置完之后，点击import此时界面变成图7图7如果在下图8勾选红框这个选项,就会出现我们刚才设定输入输出数据的曲线,如图9所示,其他勾选项是频域的分析和显示,暂不用它。图8图9看看与我们实际设置的输入输出是否符合,如果符合,那么我们离成功就不远咯,如果发现异常,那再好好检查一遍,直到确保数据导入没有问题！下面两段红色斜杠之间的内容,对于本实验,可以直接跳过,看一下对后续复杂模型的处理有好处哦,也算全面熟悉一下工具。///////////////////////////////////////////////////////////到这接着选preprocess也就是对数据进行预处理了，下拉菜单中有很多种处理方法和手段，有这个心思的人可以挨个试一下功能。图10图10预处理的对象是workingdata中的数据，每进行一种预处理在左边就会有新的数据生成，这时只要将新的数据移动到workingdata的那个方框，就可以将workingdata换成你所想处理的数据了，可以这样多次进行处理，得到你最终想用来辨识的数据和用于验证的数据（不需要的数据可以拖到那个trash里面删除，就是回收站了，也可以从回收站中找回的）接下来就是辨识了，首先把辨识用的数据拖到workingdata那个方框，再把验证的数据拖到validationdata那个方框，这个validationdata就是最原始数据稍作处理得到的一个更接近理论模型的对象数据，在这实验里，本人用的就是默认数据，也就是不做任何preprocess处理。///////////////////////////////////////////////////////////点下拉菜单estimate，选你想要的模型，在弹出的对话框中设定参数，图11图11这里边的模型种类比较多，有线性的、非线性的、状态空间的、经典传递函数形式的等等，我们所选的就是图中红色方框processmodels，单击。弹出如下界面图n传o递□函00数□00□00[-InfInf]Fit:(14-Tp1◎ProcessModelsParameterKnownValueinitialGuessBoundsMode=lTransferFunctionAuto[-InfInfJKexp(-Tdl-s)Auto[0001InfJ[0.001Inf][0.001Ini]RdlesAutoInitialGuessqAuto-selectedFromexistingmo<lel:IntegratorUa-erddisfinaidVslua—>lnrtislGusbbDisturbanceMedel:Focus:Options...Iteration]DisplayStopIterations[Zero极点个数及其他n传o递□函00数□00□00[-InfInf]Fit:(14-Tp1◎ProcessModelsParameterKnownValueinitialGuessBoundsMode=lTransferFunctionAuto[-InfInfJKexp(-Tdl-s)Auto[0001InfJ[0.001Inf][0.001Ini]RdlesAutoInitialGuessqAuto-selectedFromexistingmo<lel:IntegratorUa-erddisfinaidVslua—>lnrtislGusbbDisturbanceMedel:Focus:Options...Iteration]DisplayStopIterations[Zero极点个数及其他是否有零点、有延迟、有积分？根据勾选，传函自动调整None▼Initialstate:Auto▼Simulation▼CDvariance:Estimate▼ImprovementP1DEstimateCloseHelpName:图12在本实验中，我们做如下设置，见图13红框标出部分。图13然后点击最下方Estimate，就有模型生成了。图14图14先勾选上图绿色框选项，看到了什么？给出了拟合率，bestfits98.72有木有？？传递函数具体的数值，双击上图红色框，见图15？图15在右边的数据栏中；也就是modelviews中了，下面有很