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

最近在做一个项目的方案设计，应各位老总的要求，只有系统框图和器件选型可不行，为了凸显方案设计的高大上，必须上理论分析，炫一下“技术富”，至于具体有多大实际指导意义，那就不得而知了！本人也是网上一顿百度，再加几日探索，现在对用matlab 实现系统辨识有了一些初步的浅薄的经验，在此略做一小节。

必须要指出的是，本文研究对象是经典控制论理最简单最常用的线性时不变的siso 系统，而且是2阶的哦，也就是具有如下形式的传递函数：
1
21)(2
2++=Ts s T s G ξ 本文要做的就是，对于有这样传递函数的一个系统，要辨识得到其中的未知数T ， ξ!!这可是控制系统设计分析的基础哦，没有系统模型，啥理论、算法都是白扯，在实际工程中非常重要哦！ 经过总结研究，在得到系统阶跃响应实验数据之后（当然如果是其他响应，也有办法可以辨识，在此还是只讨论最简单的阶跃响应实验曲线，谁让你我是菜鸟呢），利用matlab 至少可以有两种方法实现实现（目前我只会两种，呵呵）！
一、函数法
二、GUI 系统辨识工具箱
下面分别作详细介绍！
一、 函数法
看官别着急，先来做一段分析（请看下面两排红*之间部分），这段分析是网上找来的，看看活跃一下脑细胞吧，如果不研读一下，对于下面matlab 程序，恐怕真的就是一头雾水咯！
*******************************************************************************
G(s)可以分解为：)
)((1)(212ωω++=s s T s G
其中， []
[]
1
1
1
1
2221--=-+=ξξωξξωT
T
1ω、2ω都是实数且均大于零。

则有：
211
ωω=T ，2
12
12ωωωωξ+=
传递函数进一步化为：
)
)(()(212
1ωωωω++=
s s s G 因此，辨识传递函数就转化为求解1ω、2ω。

当输入为单位阶跃函数时，对上式进行拉普拉斯反变换，得系统时域下的单位阶跃响应为：
t
t
e
e
t y 21
2111
221)(ωωωωωωωω---+
--
=
即 t
t
e
e
t y 21
2111
22)(1ωωωωωωωω----
-=
-
令1ω=2ωk )1(>k
，得
t
k t e
k e k k t y 221
11)(1ωω-----=-
⎥⎦
⎤⎢⎣⎡--=---t k t e k e k k 2)1(2111ωω 对上式两边取以e 为底的对数得
[]⎥⎦⎤
⎢⎣⎡-+--=---t k e k t k k t y 2)1(211ln 1ln )(1ln ωω 当∞→t 时，⎥⎦
⎤⎢⎣⎡---t k e k 2)1(11ln ω0
→，则上式化简为
[]t k k t y 21
ln )(1ln ω--=-
该式的形式满足直线方程
b at t y +=)(*
其中，
)(*
t y =[])(1ln t y -，1
ln ,2-=-=k k
b a ω)1(>k
通过最小二乘算法实现直线的拟合，得到a ，b 的值，即可得到1ω、
2ω的值，进而可得系统的传递函数。

***************************************************
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*')
打开matlab，新建一个Function,把上述程序段拷进去，保存，运行~~~~~~~~~，运行结果：
系统的传递函数为
)
4797.0)(126.1(54034
.0)(++=S S S G
很顺利吧？先高兴一个！问题接着马上就来了，上面这个例子，这个传递函数的极点刚好都是负实数，因此辨识得很顺利，但是如果系统是欠阻尼系统，也就是如果传递函数的根是复数，那么上述函数段，就无能为力咯，会出现说“matlab 无法处理增益为复数情况之类······” 例如
)
1)(1(2
)(j S j S S G -+++=
对于这个系统，若果用simulink 做一下阶跃响应，再把实验数据代入上述函数段，那就不行咯！怎么办呢，只能另辟蹊径了！
二、（System Identification Tool）系统辨识工具箱
早听说matlab博大精深，神通广大了，于是乎我确定肯定有更简单、直观、强大的工具来完成这小儿科把戏。

查资料琢磨之后，我做了个小实验，在simulink里验证了该种方法。

该方法的大原则是：在确定了系统的输入输出数据（两个列向量N×1形式，如果是1×N，会提示出错！）之后，设计好一定的辨识原则（比如说是2阶？3阶？，传递函数是零极点形式，还是带阻尼形式，等等），然后就交给强大的matlab，得到辨识结果。

Step by step，plz！
Step1、建立模型获取系统输入输出数据
图1
图1系统的输入是阶跃信号，用Scope1监视，并输出到workspace （这步不会的自己百度哦），采样周期是0.1s，得到输入变量u（101×1的矩阵）；本人在系统的阶跃响应上叠加了一白噪声，当然也可以不加噪声，加了噪声就是期望更真实的模拟实际情况，白噪声参数设置见图2
图2
同样在Scope2监视，也将结果输出到workspace，得到响应数据y（同样也是101×1的矩阵）
Step 2、进入辨识工具箱&设置辨识规则
直接在command window 输入 ident，回车，进入辨识工具箱图3
图3
点击import下拉菜单，选时域数据time domain data，见图4
图4
在下图5红色圈区域输入之前得到的系统输入和输出数据，u和y
图5
在下图6绿色圈内输入数据的一些信息，因为之前模型中，阶跃起点我是放在0s处的，这里也设置0，如果前面模型仿真是1s，这里应该也是1s；采样时间是0.1s，根据实际情况设置统一哦
图6
设置完之后，点击import此时界面变成图7
图7
如果在下图8勾选红框这个选项，就会出现我们刚才设定输入输出数据的曲线，如图9所示，其他勾选项是频域的分析和显示，暂不用它。

图8。