加载中…
Matlab系统辨识工具箱之后处理
(2013-11-18 17:48:49)
标签：
compareresidpesimpredict |
分类： 信号处理 |
转载自了凡春秋USTC https://chunqiu.blog.ustc.edu.cn/?p=386
% 模型辨识后续处理
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;
C = [1 0 0 0;
K = randn(statenum, innum)*1e-2;
x0 = [1 2 5 0]; % 初始状态
R = [0.125274,0.116642;0.116642,0.216978];
%% 信号仿真
% 输入信号长度
n = 1000;
% 模型对象
mcon = idss(A,zeros(size(K)),C,zeros(2),K,x0,'Ts',0,'NoiseVariance', R); % 连续时间模型
% 零输入
simdat0 = iddata([], zeros(n, 2), ts); % 时间序列模型做仿真，貌似只能这样指定空输入用sim仿真
% 噪声
eo = randn(n, 2);
simopt = simOptions('AddNoise', true, 'NoiseData', eo);
% 仿真
yn = sim(mcon, simdat0, simopt);
figure
for i = 1:outnum
end
%% n4sid时间序列辨识
opt = n4sidOptions('N4Horizon',[15 39 39]); % N4Horizon可以由辨识工具箱给出一个良好的估计
[syssid, x00] = n4sid(yn, 4, 'Ts', ts, opt);
%% 辨识结果查看
% 1.sim函数做时间序列初始条件仿真（无输入控制信号）
% 对于ARMA时间序列模型，可以加噪声（可参考sim帮助下Tips的代码）；对状态空间
% 模型，不能加噪声，想加噪声可以按上面信号仿真部分的代码来完成。
% 下面是辨识模型的初始条件仿真（初始条件在syssid中），只要给一个零输入即可
simdat = iddata([], zeros(n, 2), ts);
figure, sim(syssid, simdat);
% 2.predict函数做辨识模型预测
figure, predict(syssid, yn, 1);
yp1 = predict(syssid, yn, 1);
yp2 = predict(syssid, yn, 2);
figure
subplot(211)
plot([yp1.OutputData(:,1) yp2.OutputData(:,1)])
title('预测-输出1')
legend('1步预测', '2步预测')
subplot(212)
plot([yp1.OutputData(:,2) yp2.OutputData(:,2)])
title('预测-输出2')
legend('1步预测', '2步预测')
% 3.compare函数做辨识模型的效果对比
figure
compare(yn, syssid); % 省略第三个参数表示仿真对比，有第三个参数为K步预测对比
figure
compare(yn, syssid, 1); % 1步预测对比，syssid中的Fit to estimation data是1步预测fit
% 4.pe、resid函数做辨识模型预测误差分析
figure
pe(syssid, yn, 2); % 计算K步预测误差
e = resid(syssid, yn); % 计算1步预测误差，等价于e = pe(syssid, datan, 1);
figure
subplot(211)
plot(e.OutputData(:,1))
title('1步预测误差-输出1')
subplot(212)
plot(e.OutputData(:,2))
title('1步预测误差-输出2')
% resid的函数帮助中提到，resid返回的预测误差e可以用来做误差模型建模，从而了
% 解辨识残余信息的模型结构。为什么不对仿真误差建模，这点我也不清楚。
% 5.forecast函数做辨识模型预报
yf = forecast(syssid, yn, 1000);
yf1 = forecast(syssid, yn(900:end), 1000); % 只用最后的100个数据做预报
t1 = 1:size(yn, 1);
t2 = t1(end)+(1:size(yf, 1));
figure
subplot(211)
plot(t1,yn.OutputData(:,1), t2, yf.OutputData(:,1), 'r', t2, yf1.OutputData(:,1), 'g')
title('模型预报-输出1')
legend('Measured', 'Forecasted-All', 'Forecasted-100')
subplot(212)
plot(t1,yn.OutputData(:,2), t2, yf.OutputData(:,2), 'r', t2, yf1.OutputData(:,2), 'g')
title('模型预报-输出2')
legend('Measured', 'Forecasted-All', 'Forecasted-100')
% 使用全部数据预报和部分数据预报，结果相差不大，与预报算法有关。
模型仿真
模型预测
模型、数据对比
预测误差
模型预报
喜欢
0
赠金笔