Matlab 系统辨识工具箱之时间序列辨识_长笛人倚楼Gloria

http://blog.sina.com.cn/u/2214742602

首页博文目录关于我

个人资料

微博

加好友发纸条

写留言加关注

博客等级：
博客积分：

博客访问：
关注人气：
获赠金笔：0支
赠出金笔：0支
荣誉徽章：

正文字体大小：大中小

Matlab 系统辨识工具箱之时间序列辨识

(2013-11-18 17:15:35)

标签：

时间序列

n4sid

ssest

compare

pe

分类：信号处理

这篇博文是对前面介绍的辨识方法的一个小总结和应用，并转化成代码形式进行示例，其中很多函数在辨识中比较常用，如n4sid，ssest，compare，pe等，此外，还对离散系统和连续系统的A矩阵进行对比，将离散的进行转换得到系统模态。

以下转载自了凡春秋USTC https://chunqiu.blog.ustc.edu.cn/?p=384

之前写过Matlab系统辨识工具箱的博文，介绍过状态空间辨识、结构辨识、灰箱辨识等，实际上工具箱还有一个常用的辨识工具---时间序列辨识，试了试，效果蛮好的，这里总结一下！

还是以之前的弹簧质量系统为例来做一个仿真测试。时间序列辨识就是只基于输出信号的辨识，如弹簧质量系统在一个非稳态初始条件下释放，系统将出现振荡。如果这时测出系统的位移输出，如何从这个输出数据得到弹簧质量系统的信息呢？如估计系统的固有频率、阻尼比、振型等。而这可以利用时间序列辨识实现。代码如下

% 弹簧质量系统建模时间序列辨识

clc

clear

close all

ts = 0.001;

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 = [m1 0; 0 m2];

%% 转化到状态空间

% x_{k+1} = A x_k + K e_k

% y_k = C x_k + e_k

% cov(e_k) = R

innum = 2;

outnum = 2;

statenum = 4;

A = [0 1 0 0;

-(k1+k2)/m1 0 k2/m1 0;

0 0 0 1;

k2/m2 0 -(k3+k2)/m2 0];

C = [1 0 0 0;

0 0 1 0];

K = randn(statenum, innum)*1e-2;

x0 = [1 2 5 0]; % 初始状态

R = [0.125274,0.116642;0.116642,0.216978];

% 计算系统固有频率、阻尼比

fprintf('系统固有频率、阻尼比：\n');

[f, xi] = CalcModeParameters(A, 0);

%% 信号仿真

% 输入信号长度

n = 1000;

% 模型对象

mcon = idss(A,zeros(size(K)),C,zeros(2),K,x0,'Ts',0,'NoiseVariance', R); % 连续时间模型

% 空输入

simdat0 = iddata([], zeros(n, 2), ts);

% 噪声

eo = randn(n, 2);

simopt = simOptions('AddNoise', true, 'NoiseData', eo);

% 仿真

yn = sim(mcon, simdat0, simopt);

% 绘制仿真结果

figure

for i = 1:outnum

subplot(outnum, 1, i)

plot(eo(:, i), 'b')

hold on

plot(yn.OutputData(:, i), 'r')

axis tight

title(sprintf('输出%d', i))

legend({'噪声序列', '噪声仿真'})

end

%% n4sid时间序列辨识

opt = n4sidOptions('N4Horizon',[15 39 39]); % N4Horizon可以由辨识工具箱给出一个良好的估计

[syssid, x0] = n4sid(yn, 4, 'Ts', ts, opt);

fprintf('\nn4sid辨识得到的固有频率、阻尼比：\n');

[f, xi] = CalcModeParameters(syssid.A, ts);

% 辨识结果查看

figure

compare(yn, syssid); % 省略第三个参数表示仿真compare，否则表示K步预测compare

figure

pe(syssid, yn, 1); % 计算K步预测误差

%% ssest时间序列辨识

sysss = ssest(yn, 4, 'Ts', ts);

fprintf('\nssest辨识得到的固有频率、阻尼比：\n');

[f, xi] = CalcModeParameters(sysss.A, ts);

% 辨识结果查看

figure

compare(yn, sysss); % 省略第三个参数表示仿真compare，否则表示K步预测compare

figure

pe(sysss, yn, 1); % 计算K步预测误差

其中根据系统矩阵A计算模态参数的函数如下

% 据状态空间矩阵A计算模态参数--固有频率、阻尼比

% 对于离散状态空间，额外求出Psi，可用于求振型

function [f, xi, Psi] = CalcModeParameters(A, ts)

if ts == 0 % 连续时间状态空间矩阵

d = eig(A);

% 连续时不求Psi，如果想求，则用下面的转化关系求

% Ad = exp(A*ts); % 连续时间A转换为离散时间A

% [Psi, ~] = eig(A);

else % 离散时间状态空间矩阵

[Psi, D] = eig(A);

d = log(diag(D))/ts; % 离散时间特征值向连续时间特征值转化

end

len = length(d);

ftmp = zeros(len, 1);

xi = zeros(len, 1);

for i = 1:len

lambda = d(i);

ftmp(i) = abs(lambda)/2/pi; % 固有频率

xi(i) = abs(real(lambda))/abs(lambda); % 阻尼比

end

[f, ia] = unique(eval_r(vpa(ftmp, 10)));

xi = xi(ia);

for i = 1:length(f)

fprintf('固有频率%d:%fHz\n', i, f(i));

fprintf('阻尼比%d:%f\n', i, xi(i));

end

运行程序，结果如下

系统固有频率、阻尼比：

固有频率1:9.915822Hz

阻尼比1:0.000000

固有频率2:22.853200Hz

阻尼比2:0.000000

n4sid辨识得到的固有频率、阻尼比：

固有频率1:9.944764Hz

阻尼比1:0.004017

固有频率2:22.850819Hz

阻尼比2:0.000716

ssest辨识得到的固有频率、阻尼比：

固有频率1:9.911302Hz

阻尼比1:0.000138

固有频率2:22.850265Hz

阻尼比2:0.000099

可以看到辨识结果还是很准确的，而且ssest迭代辨识比n4sid辨识要准一些。相应的结果图如下，从中可以更直观地看出辨识的效果：

第一幅图为仿真添加的噪声和系统仿真信号

第二幅图为n4sid辨识结果与仿真信号的吻合程度，第三幅为n4sid的辨识误差

第四、五副图为ssest辨识的结果和误差

其实这里使用的时间序列辨识，其理论基础就是著名的随机子空间辨识方法，可以看到此算法的有效性。

你可能会注意到在n4sid辨识中，有个设置选项

opt = n4sidOptions('N4Horizon',[15 39 39]);

这个[15 39 39]参数向量是怎么来的呢？实际上，这个向量对n4sid的辨识结果影响是挺大的，如果你发现n4sid辨识结果不佳，可以调整这个参数向量来查看不同参数的效果，之前博文中有提到过这个问题。如果自己拿不定主意，可以让matlab帮你选一个不错的参数，这可以通过系统辨识工具箱GUI来完成，其实，对于常见的系统辨识问题，建议都先用GUI快速辨识，对于GUI不便完成的，如结构化辨识、灰箱辨识等，再用命令行辨识，而且matlab提供了从GUI向命令行的代码导出功能，很实用，这个[15 39 39]就是从GUI辨识结果中获取的，过程如下

在状态空间辨识中设置Orders如下，然后进行辨识

然后使用推荐的（红色）阶次，点击Insert，插入模型

这时会弹出一个Order编辑框，其下有使用的辨识算法的推荐参数设置，这个参数一般是比较不错的。

来源： <https://chunqiu.blog.ustc.edu.cn/?p=384>

来自为知笔记(Wiz)

阅读┊ 收藏 ┊ 喜欢 ▼ ┊打印┊举报/Report

前一篇：状态空间系统加噪声

后一篇：Matlab系统辨识工具箱之后处理

新浪BLOG意见反馈留言板　欢迎批评指正