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

图7 如果在下图 8 勾选红框这个选项， 就会出现我们刚才设定输入输出数 据的曲线， 如图 9 所示， 其他勾选项是频域的分析和显示， 暂不用它。
图8
图9 看看与我们实际设置的输入输出是否符合，如果符合，那么我们 离成功就不远咯，如果发现异常，那再好好检查一遍，直到确保数据 导入没有问题！ 下面两段红色斜杠之间的内容，对于本实验，可以直接跳过，看 一下对后续复杂模型的处理有好处哦，也算全面熟悉一下工具。 /////////////////////////////////////////////////////////// 到这接着选 preprocess 也就是对数据进行预处理了， 下拉菜单中 有很多种处理方法和手段，有这个心思的人可以挨个试一下功能。图 10
最近在做一个项目的方案设计，应各位老总的要求，只有系统框 图和器件选型可不行， 为了凸显方案设计的高大上， 必须上理论分析， 炫一下“技术富” ，至于具体有多大实际指导实现系统 辨识有了一些初步的浅薄的经验，在此略做一小节。 必须要指出的是，本文研究对象是经典控制论理最简单最常用的 线性时不变的 siso 系统，而且是 2 阶的哦，也就是具有如下形式的 传递函数：
G(s)可以分解为： G ( s )

1 T 2 ( s 1 )( s 2 )
其中，
1 2 1 T 1 2 2 1 T
1



1 、2 都是实数且均大于零。
则有：
1 2 T ， 1 2 2 12
1
传递函数进一步化为：
即
1 y (t )
2 1 e 1t e 2 t 2 1 2 1
1) ，得
令 1 = k 2 ( k
1 y (t )
k 2 t 1 k 2 t e e k 1 k 1
k 2t 1 ( k 1) 2t e 1 e k 1 k
G (s)
1 2 ( s 1 )( s 2 )
因此，辨识传递函数就转化为求解 1 、 2 。 当输入为单位阶跃函数时，对上式进行拉普拉斯反变换，得系统时域 下的单位阶跃响应为：
y (t ) 1
2 1 e 1t e 2 t 2 1 2 1
图4 在下图 5 红色圈区域输入之前得到的系统输入和输出数据，u 和 y
图5
在下图 6 绿色圈内输入数据的一些信息，因为之前模型中，阶跃起 点我是放在 0s 处的，这里也设置 0，如果前面模型仿真是 1s，这里 应该也是 1s；采样时间是 0.1s，根据实际情况设置统一哦
图6 设置完之后，点击 import 此时界面变成图 7
图1 图 1 系统的输入是阶跃信号， 用 Scope1 监视， 并输出量 u（101 ×1 的矩阵） ；本人在系统的阶跃响应声就是期望更真实的模拟实际情况，白噪声参数 设置见图 2
G(S )
2 (S 1 j )(S 1 j )
对于这个系统，若果用 simulink 做一下阶跃响应，再把实验数据 代入上述函数段，那就不行咯！怎么办呢，只能另辟蹊径了！
二、 （System Identification Tool）系统辨识工具箱 早听说 matlab 博大精深，神通广大了，于是乎我确定肯定有更简 单、直观、强大的工具来完成这小儿科把戏。查资料琢磨之后，我做 了个小实验，在 simulink 里验证了该种方法。该方法的大原则是： 在确定了系统的输入输出数据 （两个列向量 N×1 形式， 如果是 1×N， 会提示出错！ ） 之后， 设计好一定的辨识原则 （比如说是 2 阶？3 阶？， 传递函数是零极点形式，还是带阻尼形式，等等） ，然后就交给强大 的 matlab，得到辨识结果。Step by step，plz！ Step1、 建立模型获取系统输入输出数据
图 12 在本实验中，我们做如下设置，见图 13 红框标出部分。
图 13
然后点击最下方 Estimate，就有模型生成了。图 14
图 14 先勾选上图绿色框选项，看到了什么？给出了拟合率，best fits 98.72 有木有？？传递函数具体的数值，双击上图红色框，见图 15？
图 15
在右边的数据栏中；也就是 model views 中了，下面有很多可以 选择，每选一个就可以生成一幅对应的图，是由用于验证的数据生成 的。把模型拖到 to workspace 那个方框，再去看 workspace，多的 那个变量就是你所辨识出的模型了
对上式两边取以 e 为底的对数得
ln1 y (t ) ln
当t
k 1 2 t ln 1 e ( k 1) 2 t k 1 k
1 时， ln 1 e ( k 1)2t 0 ，则上式化简为 k
ln1 y (t ) ln k 2t k 1
图 10 预处理的对象是 working data 中的数据，每进行一种预处理在左 边就会有新的数据生成，这时只要将新的数据移动到 working data 的那个方框，就可以将 working data 换成你所想处理的数据了，可 以这样多次进行处理， 得到你最终想用来辨识的数据和用于验证的数 据（不需要的数据可以拖到那个 trash 里面删除，就是回收站了，也 可以从回收站中找回的） 接下来就是辨识了，首先把辨识用的数据拖到 working data 那个 方 框 ， 再 把 验 证 的 数 据 拖 到 validation data 那 个 方 框 ， 这 个 validation data 就是最原始数据稍作处理得到的一个更接近理论模 型的对象数据，在这实验里，本人用的就是默认数据，也就是不做任 何 preprocess 处理。 ///////////////////////////////////////////////////////////
图2 同样在 Scope2 监视，也将结果输出到 workspace，得到响应数据 y（同样也是 101×1 的矩阵） Step 2、进入辨识工具箱&设置辨识规则 直接在 command window 输入 ident，回车，进入辨识工具箱图 3
图3
点击 import 下拉菜单，选时域数据 time domain data，见图 4
注意：在某个模型或某组数据上点一下，线变细了就不会在图中显示 出来了！ ！ ！ ！ ！ ！ ！ ！
最后再一次提出，上述分析并不一定完全正确，可能有些概念并 不清晰可靠，当然还有很多功能并没有被发掘，可以确定的是解决这 个问题的两个大方向没有问题，若想做到精益求精，还需要再仔细研 究斟酌哦！
1 G (s) 2 2 T s 2Ts 1
本文要做的就是，对于有这样传递函数的一个系统，要辨识得到 其中的未知数 T， ξ!!这可是控制系统设计分析的基础哦，没有系 统模型，啥理论、算法都是白扯，在实际工程中非常重要哦！ 经过总结研究，在得到系统阶跃响应实验数据之后（当然如果是 其他响应，也有办法可以辨识，在此还是只讨论最简单的阶跃响应实 验曲线，谁让你我是菜鸟呢） ，利用 matlab 至少可以有两种方法实现 实现（目前我只会两种，呵呵） ！
一、函数法 二、GUI 系统辨识工具箱
下面分别作详细介绍！
一、 函数法 看官别着急，先来做一段分析（请看下面两排红*之间部分） ，这 段分析是网上找来的，看看活跃一下脑细胞吧，如果不研读一下，对 于下面 matlab 程序，恐怕真的就是一头雾水咯！
*******************************************************************************
Matlab 程序代码:
clc close all t=[1 3 5 7 9 11 13 15 17 19]; y=[0.149086 0.5890067 0.830617 0.933990 0.973980 0.991095 0.995868 0.998680 0.999490 0.999850]; y2=log(1-y); plot(t,y2,'*'); grid on pm=polyfit(t,y2,1) value=polyval(pm,t); hold on plot(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([0 20 0 1.2]) hold on plot(t,y,'r*')
点下拉菜单 estimate，选你想要的模型，在弹出的对话框中设定参 数,图 11
图 11 这里边的模型种类比较多，有线性的、非线性的、状态空间的、 经典传递函数形式的等等，我们所选的就是图中红色方框 process models，单击。弹出如下界面图 12！
极点个 数及其 他
传 递 函 数
是否有零点、 有 延迟、 有积分？ 根据勾选， 传函 自动调整
打开 matlab，新建一个 Function,把上述程序段拷进去，保存，运行 ~~~~~~~~~，运行结果：
系统的传递函数为
G( S )
0.54034 ( S 1.126)(S 0.4797)
很顺利吧？先高兴一个！问题接着马上就来了，上面这个例子， 这个传递函数的极点刚好都是负实数，因此辨识得很顺利，但是如果 系统是欠阻尼系统，也就是如果传递函数的根是复数，那么上述函数 段，就无能为力咯，会出现说“matlab 无法处理增益为复数情况之 类· · · · · · ” 例如