Matlab系统辨识尝试之详细过程1
前面介绍了Matlab系统辨识工具箱的一些用法，这里拿一个直观的例子来尝试工具箱的具体用法。
比较长，给个简单目录吧：
1.辨识的准备
2.辨识数据结构的构造
3.GUI辨识
4.辨识效果
5.对固有频率的辨识
6.结构化辨识
7.灰箱辨识
8.加入kalman滤波的灰箱辨识
1.辨识的准备
在辨识前，首先要根据自己辨识的情况，确定要辨识的状态空间模型的一些特点，如连续还是离散的；有无直通
分量（即从输入直通到输出的分量）；输入延迟；初始状态等。
了解了这些情况就可以更快速的配置辨识时的一些设
置选项。
2.辨识数据结构的构造
使用原始数据构造iddata结构：
data=iddata(y,u,Ts);
这里以一个弹簧质量系统的仿真为例
代码如下，其中用到了函数MDOFSolve，这在之前的博文介绍过（/?p=183），拿来用即可。
如果发现运行有错误，可以将MDOFSolve函数开头的一句
omega2=real(eval(omega2));
注释掉。
%弹簧质量系统建模
clc
clear
close all
m=200;
k=980*1000;
c=1.5*1000;
m1=1*m;
m2=1.5*m;
k1=1*k;
k2=2*k;
k3=k1;
%%由振动力学知识求固有频率
M=[m10;0m2];
K=[k1+k2-k2;-k2k3+k2];
[omega,phi,phin]=MDOFSolve(M,K);
fprintf('固有频率:%fHz\n',subs(omega/2/pi));
%%转化到状态空间
innum=2;
outnum=2;
statenum=4;
A=[0100;
-(k1+k2)/m10k2/m10;
0001;
k2/m20-(k3+k2)/m20];
B=[00;
1/m10;
00;
01/m2];
C=[1000;
0010];
D=zeros(outnum,innum);
K=zeros(statenum,innum);
mcon=idss(A,B,C,D,K,'Ts',0);%连续时间模型
figure
impulse(mcon)
%%信号仿真，构造数据供辨识
n=511;%输入信号长度
Ts=0.001;
t=0:Ts:(n-1)*Ts;
u1=idinput(n,'prbs');%输入1为伪随机信号
u2=zeros(n,1);%输入2为空
u=[u1u2];
simdat=iddata([],u,Ts);%形成输入数据对象
e=randn(n,2)*1e-7;
simopt=simOptions('AddNoise',true,'NoiseData',e);%添加噪声yn=sim(mcon,simdat,simopt);%加噪声仿真
y=sim(mcon,simdat);%无噪声仿真
figure
for i=1:outnum
subplot(outnum,1,i)
plot(t,y.OutputData(:,i))
hold on
plot(t,yn.OutputData(:,i),'r')
axis tight
title(sprintf('输出%d',i))
legend({'无噪声仿真','含噪声仿真'})
end
%保存输入输出数据，供后续辨识
data=iddata(y.OutputData,simdat.InputData,Ts);
datan=iddata(yn.OutputData,simdat.InputData,Ts);
运行后，变量data中保存了无噪声的系统仿真输入输出数据，datan中为含噪声的仿真数据。
产生的图形如下固有频率:9.915822Hz
固有频率:22.853200Hz
3.GUI辨识
辨识可以从GUI开始，GUI的辨识对应的代码可以导出，方便以后更便捷的调整参数。
3.1输入ident打开系统辨识工具箱
3.2数据导入到GUI，如下
3.3原始数据预处理，主要是去趋势detrend：
3.4开始状态空间辨识
3.5先做阶次选择，红框标记的为建议设置，其他设置根据自己的模型特性做配置
点击Estimate做阶次估计
可以看到4阶后落差明显，所以4阶是最好的阶次，插入此辨识模型，查看辨识结果如下，可以看到辨识结果很好（因为原始数据中没加入噪声）。
在辨识主窗口双击辨识结果，可以看到辨识结果，点击Present可以将结果显示出来，Diary and Notes里有辨识的相关代码，可以拿来进一步做修改。
上面有个N4Horizon参数设置，它表示N4SID算法中使用的前向和后向预测界，意义如下：
它是一个包含3个元素的向量[r sy su]，其中r是最大前向预测界，表示使用多达r步的前向预测器；sy、su表示预测器使用过去输出和过去输入的个数。
这些值的设置对辨识结果有实质性的影响，但没有简单的准则指导如何选择它们。
如果指定N2Horizon为一个k*3的矩阵，则可以测试不同参数的效果，结果会给出最佳预测模型。
如果设置成一个标量，则表示3个值相同。
如果自己不知如何设置这些参数，可以让N2Horizon='auto'，这样软件会根据AIC准则选择sy和su的值。
'auto'也是默认值。
3.6辨识的改善
后面可以尝试调整参数改善辨识效果，如使用PEM迭代辨识方法、改变辨识的状态方程形式等。
下面总结一些常用的改善辨识效果的方法：
3.6.1PEM或N4SID方法的选择
N4SID为子空间方法，可以先用来估计得到一个初始模型；PEM为最小化预测误差迭代方法（得到最大似然估计），对应的函数为ssest，可以对初始估计进行改善。
PEM会默认先运行N4SID得到初始估计矩阵再对初始估计矩阵进行迭代优化，所以它运行时间较长。
N4SID运行比SSEST快，但精度和鲁棒性均不如SSEST，需要指定额外的难于调整的参数。
3.6.2Focus的选择
这在之前博文介绍过，这里简单提一下：
Prediction：使用噪声模型H的逆来衡量不同频率范围内拟合的相对重要性。
这与最小化一步前向预测相对应，适合于在一段短时间间隔的拟合。
适合于输出预测应用。
Simulation：使用输入的频谱来衡量不同频率范围内拟合的相对重要性，而不使用噪声模型来衡量。
适合于输出仿真应用。
Stability：适合于估计最稳定的模型。
Filter：通过弹出Estimation Focus对话框来指定一个用户滤波器。
干扰模型由估计数据决定。
3.6.3Initial state
指定对待初始状态的方法。
这个在使用PEM方法辨识时作为迭代的初始状态使用。
如果辨识得到不精确的解，可以尝试指定特定的方法来处理初始状态而不是选择自动。
Auto：基于估计的数据自动选择Zero、Estimate或Backcast。
Zero：如果初始状态对估计误差的效果微不足道，可以设置为0来最优化方法的效果。
Estimate：把初始状态作为一个未知的参数向量从数据中辨识。
Backcast：用后向滤波方法（最小二乘拟合）估计初始状态。
3.6.4K的辨识
当不确定如何配置K时，可以先不辨识它，只辨识系统矩阵。
在得到系统模型后，再使用ssest配置K为自由参数辨识改善模型。
如
init_sys=ssest(data,n,'DisturbanceModel','none');
init_sys.Structure.k.Free=true;
sys=ssest(data,init_sys);
init_sys为无噪声的动态模型。
4.辨识效果
上面讲的过程中是对无噪声数据处理的，效果很好，这里对各辨识算法的效果做对比。
对无噪声数据，n4sid和ssest辨识的效果如下
可以看到对无噪声数据两种算法的效果都非常好，都能达到准确的辨识。
对含噪声数据，n4sid和ssest辨识的效果如下
误差曲线图如下
可以看到n4sid算法对含噪声的数据辨识相对较差，而ssest辨识效果很好，基本把原系统辨识出来了。