matlab系统辨识尝试之详细过程1,Matlab系统辨识尝试之详细过程1

Matlab系统辨识尝试之详细过程1

2018-11-27

前面介绍了Matlab系统辨识工具箱的一些用法，这里拿一个直观的例子来尝试工具箱的具体用法。比较长，给个简单目录吧：1. 辨识的准备

2. 辨识数据结构的构造

3. GUI辨识

4. 辨识效果

5. 对固有频率的辨识

6. 结构化辨识

7. 灰箱辨识

8. 加入kalman滤波的灰箱辨识

1.辨识的准备

在辨识前，首先要根据自己辨识的情况，确定要辨识的状态空间模型的一些特点，如连续还是离散的；有无直通分量(即从输入直通到输出的分量)；输入延迟；初始状态等。了解了这些情况就可以更快速的配置辨识时的一些设置选项。

2.辨识数据结构的构造

使用原始数据构造iddata结构：

data = iddata(y, u, Ts);

这里以一个弹簧质量系统的仿真为例

代码如下，其中用到了函数MDOFSolve，这在之前的博文介绍过(http://chunqiu.blog.ustc.edu.cn/?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 = [m1 0; 0 m2];

K = [k1+k2 -k2; -k2 k3+k2];

[omega, phi, phin] = MDOFSolve(M, K);

fprintf("固有频率:%fHz\n", subs(omega/2/pi));

%% 转化到状态空间

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];

B = [0 0;

1/m1 0;

0 0;

0 1/m2];

C = [1 0 0 0;

0 0 1 0];

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 = [u1 u2];

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.1 PEM或N4SID方法的选择

N4SID为子空间方法，可以先用来估计得到一个初始模型；PEM为最小化预测误差迭代方法(得到最大似然估计)，对应的函数为ssest，可以对初始估计进行改善。PEM会默认先运行N4SID得到初始估计矩阵再对初始估计矩阵进行迭代优化，所以它运行时间较长。

N4SID运行比SSEST快，但精度和鲁棒性均不如SSEST，需要指定额外的难于调整的参数。

3.6.2 Focus的选择

这在之前博文介绍过，这里简单提一下：

Prediction：使用噪声模型H的逆来衡量不同频率范围内拟合的相对重要性。这与最小化一步前向预测相对应，适合于在一段短时间间隔的拟合。适合于输出预测应用。

Simulation：使用输入的频谱来衡量不同频率范围内拟合的相对重要性，而不使用噪声模型来衡量。适合于输出仿真应用。

Stability：适合于估计最稳定的模型。

Filter：通过弹出Estimation Focus对话框来指定一个用户滤波器。干扰模型由估计数据决定。

3.6.3 Initial state

指定对待初始状态的方法。这个在使用PEM方法辨识时作为迭代的初始状态使用。如果辨识得到不精确的解，可以尝试指定特定的方法来处理初始状态而不是选择自动。

Auto：基于估计的数据自动选择Zero、Estimate或Backcast。

Zero：如果初始状态对估计误差的效果微不足道，可以设置为0来最优化方法的效果。

Estimate：把初始状态作为一个未知的参数向量从数据中辨识。

Backcast：用后向滤波方法(最小二乘拟合)估计初始状态。

3.6.4 K的辨识

当不确定如何配置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辨识效果很好，基本把原系统辨识出来了。

免责声明：本文仅代表文章作者的个人观点，与本站无关。其原创性、真实性以及文中陈述文字和内容未经本站证实，对本文以及其中全部或者部分内容文字的真实性、完整性和原创性本站不作任何保证或承诺，请读者仅作参考，并自行核实相关内容。

http://www.pinlue.com/style/images/nopic.gif