各种模型识别方法的比较
此示例显示了系统识别工具箱™中可用的几种识别方法。我们从模拟实验数据开始,并使用几种估计技术从数据中估计模型。本例中说明了以下估计例程:水疗中心
,党卫军
,特遣部队
,arx
,oe
,armax
而且bj
.
系统描述
一个噪声损坏的线性系统可以描述为:
y = u + he
在哪里y
是输出和u
输入是,而e
表示未测量(白)噪声源。G
系统的传递函数是和吗H
给出了噪声干扰的描述。
估计的方法有很多种G
而且H
.本质上它们对应于不同的参数化这些函数的方法。
定义模型
“系统识别工具箱”为用户提供了模拟从物理过程中获得的数据的选项。让我们用一组给定的系数来模拟来自IDPOLY模型的数据。
B = [0 1 0.5];A = [1 -1.5 0.7];m0 = idpoly(A,B,[1 -1 0.2],“t”, 0.25,“变量”,“问^ 1”);采样时间为0.25 s。
第三个论点[1 -1 0.2]
给出作用于系统上的扰动的特征。执行这些命令会产生以下离散idpoly模型:
m0
m0 =离散时间ARMAX模型:A(q)y(t) = B(q)u(t) + C(q)e(t) A(q) =1 - 1.5 q^-1 + 0.7 q^-2 B(q) = q^-1 + 0.5 q^-2 C(q) =1 - q^-1 + 0.2 q^-2采样时间:0.25秒参数化:多项式阶数:na=2 nb=2 nc=2 nk=1自由系数数:6使用"polydata"、"getpvec"、"getcov"表示参数及其不确定性。现状:由直接构建或改造而产生。不估计。
在这里问
表示移位算子,因此A(q)y(t)是y(t) - 1.5 y(t-1) + 0.7 y(t-2)的简称。模型展示m0
表明它是一个ARMAX模型。
这种特殊的模型结构被称为ARMAX模型,其中AR(自回归)指的是a多项式,MA(移动平均)指的是噪声c多项式,X是“额外”输入B(q)u(t)。
用一般的传递函数表示G
而且H
,该模型对应一个参数化:
G(q) = B(q)/A(q)
而且H(q) = C(q)/A(q)
,有共同点
模拟模型
我们生成一个输入信号u
并模拟模型对这些输入的响应。的idinput
命令可用于创建模型的输入信号iddata
命令将把信号打包成合适的格式。的sim卡
命令可以用来模拟输出,如下所示:
prerng = rng(12,“v5normal”);U = idinput(350,苏格兰皇家银行的);生成一个长度为350的随机二进制信号U = iddata([], U,0.25);创建一个IDDATA对象。采样时间为0.25秒。Y = sim(m0,u,“噪音”)模拟模型对此的响应
y = 350个样本的时域数据集。采样时间:0.25秒名称:m0输出单位(如指定)y1
%输入加高斯白噪声e%,根据模型描述m0rng (prevRng);
需要注意的一个方面是信号u
而且y
都是IDDATA对象。接下来,我们收集这个输入-输出数据以形成一个iddata对象。
Z = [y,u];
要获得数据对象的信息(现在包含输入和输出数据样本),只需在MATLAB®命令窗口中输入它的名称:
z
z = 350个样本的时域数据集。采样时间:0.25秒名称:m0输出单位(如设置)y1输入单位(如设置)u1
绘制输入的前100个值u
和输出y
,使用情节
在iddata对象上执行命令:
H = gcf;集(h,“DefaultLegendLocation”,“最佳”) h.Position = [100 100 780 520];情节(z (1:10 0));
为了估计的目的,只使用一部分数据是一个很好的做法,泽
并保存另一部分,以便以后验证估计的模型:
Ze = z(1:200);Zv = z(201:350);
进行光谱分析
现在我们已经获得了模拟数据,我们可以估计模型并进行比较。让我们从光谱分析开始。我们调用水疗中心
命令,该命令返回频率函数和噪声频谱的频谱分析估计。
GS = spa(ze);
频谱分析的结果是频率的复值函数:频率响应。它被打包到一个名为IDFRD
对象(已识别的频率响应数据):
GS
GS = IDFRD模型。包含1个输出(s)和1个输入(s)的频率响应数据,以及输出扰动的频谱。响应数据和扰动谱可在128个频率点上获得,范围从0.09817 rad/s到12.57 rad/s。采样时间:0.25秒输出通道:“y1”输入通道:“u1”状态:使用SPA对时域数据“m0”进行估计。
为了绘制频率响应,习惯上使用波德图波德
或bodeplot
命令:
clf h = bodeploy (GS);% bodeploy返回一个plot句柄,而bode没有Ax =轴;轴([0.1 10 ax(3:4)])
估计GS
是不确定的,因为它是由噪声损坏的数据形成的。谱分析方法也提供了它自己对这种不确定性的评估。要显示不确定度(例如3个标准差),我们可以使用showConfidence
命令在情节句柄上h
由前一个bodeplot
命令。
showConfidence (h, 3)
该图表示,尽管频率响应(蓝线)的名义估计不一定准确,但有99.7%的概率(正态分布的3个标准差)真实响应位于阴影(浅蓝色)区域内。
估计参数状态空间模型
接下来让我们估计一个默认的线性模型(不指定特定的模型结构)。它将被计算为使用预测误差方法的状态空间模型。我们使用党卫军
本例中的函数为:
M = sest(ze)%模型的顺序将自动选择
m =连续时间识别状态空间模型:dx/dt = A x(t) + B u(t) + K e(t) y(t) = C x(t) + D u(t) + e(t) A = x1 x2 x1 -0.007167 1.743 x2 -2.181 -1.337 B = u1 x1 0.09388 x2 0.2408 C = x1 x2 y1 47.34 -14.4 D = u1 y1 0 K = y1 x1 0.04108 x2 -0.03751参数化:FREE形式(所有系数在A, B, C中自由)。参数及其不确定性使用“idssdata”、“getpvec”、“getcov”表示。状态:使用SSEST对时域数据“m0”进行估计。与估计数据拟合:75.08%(预测焦点)FPE: 0.9835, MSE: 0.9262
在这一点上,状态空间模型矩阵没有告诉我们很多。我们可以比较模型的频响米
与光谱分析的结果进行了比较GS
.请注意,GS
离散时间模型是while吗米
是一个连续时间模型。它是可以使用的党卫军
也可以计算离散时间模型。
bodeploy (m,GS) ax =轴;轴([0.1 10 ax(3:4)])
我们注意到这两个频响是接近的,尽管它们是用非常不同的方式估计的。
简单传递函数的估计
有各种各样的线性模型结构,都对应于描述u和y之间关系的线性差分或微分方程。不同的结构都对应于对噪声影响的各种建模方法。这些模型中最简单的是传递函数和ARX模型。传递函数的形式为:
Y(s) = NUM(s)/DEN(s) U(s) + E(s)
其中NUM(s)和DEN(s)是分子和分母多项式,Y(s)、U(s)和E(s)分别是输出、输入和误差信号(Y(t)、U(t)和E(t))的拉普拉斯变换。NUM和DEN可以通过表示零和极点数的顺序来参数化。对于给定数量的极点和零,我们可以估计传递函数使用特遣部队
命令:
MTF = tfest(ze, 2,2)有两个零点和两个极点的%传递函数
mtf =从输入“u1”输出“日元”:-0.05428 s ^ 2 - 0.02386 + 29.6 ------------------------------- s ^ 2 + 1.361 + 4.102连续时间确定传递函数。参数化:极点数:2零数:2自由系数数:5参数及其不确定度使用“tfdata”、“getpvec”、“getcov”。状态:使用TFEST对时域数据“m0”进行估计。拟合估计数据:71.69% FPE: 1.282, MSE: 1.195
ARX模型表示为:A(q)y(t) = B(q)u(t) + e(t),
或者从长远来看,
这个模型的系数是线性的。系数 , 、…… , , . .可以有效地计算使用最小二乘估计技术。具有规定结构的ARX模型的估计值-两个极点,一个零和输入上的单个滞后可得到([na nb nk] = [2 2 1]):
Mx = arx(ze,[2 2 1])
mx =离散时间ARX模型:A(z)y(t) = B(z)u(t) + e(t) A(z) =1 - 1.32 z^-1 + 0.5393 z^-2 B(z) = 0.9817 z^-1 + 0.4049 z^-2采样时间:0.25秒参数化:多项式阶数:na=2 nb=2 nk=1自由系数数:4使用"polydata"、"getpvec"、"getcov"表示参数及其不确定性。状态:在时域数据“m0”上使用ARX估计。拟合估计数据:68.18%(预测焦点)FPE: 1.603, MSE: 1.509
比较估计模型的性能
这些模型(mtf
,mx
)比默认的状态空间模型好或差米
?一种方法是用验证数据的输入来模拟每个模型(无噪声)zv
并将相应的模拟输出与集合中的实测输出进行比较zv
:
比较(zv, m, mtf, mx)
拟合是由模型解释的输出变化的百分比。很明显米
这是更好的模式吗mtf
接近。离散时间传递函数可以用任意一种方法估计特遣部队
或者是oe
命令。前者创建的是IDTF模型,后者创建的是输出误差结构(y(t) = B(q)/F(q)u(t) + e(t))的IDPOLY模型,但两者在数学上是等价的。
Md1 = tfest(ze,2,2,“t”, 0.25)%两个极点和两个零(算作z^-1多项式的根)
md1 =从输入“u1”到输出“y1”:0.8383 z^-1 + 0.7199 z^-2 ---------------------------- 1 - 1.497 z^-1 + 0.7099 z^-2采样时间:0.25秒离散时间识别传递函数。参数化:极点数:2零数:2自由系数数:4参数及其不确定度使用“tfdata”、“getpvec”、“getcov”。状态:使用TFEST对时域数据“m0”进行估计。拟合估计数据:71.46% FPE: 1.264, MSE: 1.214
Md2 = oe(ze,[2 2 1])
md2 =离散时间OE模型:y(t) = [B(z)/F(z)]u(t) + e(t) B(z) = 0.8383 z^-1 + 0.7199 z^-2 F(z) =1 - 1.497 z^-1 + 0.7099 z^-2采样时间:0.25秒参数化:多项式阶数:nb=2 nf=2 nk=1自由系数数:4参数及其不确定性使用"polydata"、"getpvec"、"getcov"。状态:在时域数据“m0”上使用OE估计。拟合估计数据:71.46% FPE: 1.264, MSE: 1.214
比较(zv, md1, md2)
这些模型md1
而且md2
为数据提供相同的匹配。
残留分析
深入了解模型质量的另一种方法是计算“残差”:即那部分e
在y = Gu + He
这不能用模型来解释。理想情况下,这些应该与输入不相关,也不相互相关。残差的计算和他们的相关性质显示渣油
命令。该函数可用于时域和频域的残差评估。首先让我们在时域内获得Output-Error模型的残差:
渣油(zv md2)%表示残差分析结果
我们看到残差和投入之间的交叉相关在置信区域,说明不存在显著相关性。动力学的估计G
应该被认为是足够的。的(自动)相关性e
很重要,所以e
不能被视为白噪音。这意味着噪声模型H
是不够的。
估计ARMAX和Box-Jenkins模型
现在让我们计算一个二阶ARMAX模型和一个二阶Box-Jenkins模型。ARMAX模型具有与仿真模型相同的噪声特性m0
而Box-Jenkins模型允许更一般的噪声描述。
Am2 = armax(ze,[2 2 2 1])%二阶ARMAX模型
am2 =离散时间ARMAX模型:A(z)y(t) = B(z)u(t) + C(z)e(t) A(z) =1 - 1.516 z^-1 + 0.7145 z^-2 B(z) = 0.982 z^-1 + 0.5091 z^-2 C(z) =1 - 0.9762 z^-1 + 0.218 z^-2采样时间:0.25秒参数化:多项式阶数:na=2 nb=2 nc=2 nk=1自由系数数:6使用“polydata”、“getpvec”、“getcov”表示参数及其不确定性。状态:在时域数据“m0”上使用ARMAX估计。与估计数据拟合:75.08%(预测焦点)FPE: 0.9837, MSE: 0.9264
Bj2 = bj(ze,[2 2 2 2 1])%二阶BOX-JENKINS模型
bj2 =离散时间BJ模型:y(t) = [B(z)/F(z)]u(t) + [C(z)/D(z)]e(t) B(z) = 0.9922 z^-1 + 0.4701 z^-2 C(z) =1 - 0.6283 z^-1 - 0.1221 z^-2 D(z) =1 - 1.221 z^-1 + 0.3798 z^-2 F(z) =1 - 1.522 z^-1 + 0.7243 z^-2采样时间:0.25秒参数化:多项式阶数:nb=2 nc=2 nd=2 nf=2 nk=1自由系数数:8参数及其不确定性使用"polydata"、"getpvec"、"getcov"。状态:在时域数据“m0”上使用BJ估计。拟合估计数据:75.47%(预测焦点)FPE: 0.9722, MSE: 0.8974
比较估计模型-模拟和预测行为
现在我们已经估计了这么多不同的模型,让我们进行一些模型比较。这可以通过比较之前的模拟输出来实现:
clf比较(zv am2 md2, bj2, mx, mtf, m)
我们还可以比较各种估计模型预测输出的能力,比如提前一步:
比较(zv、am2 md2、bj2 mx, mtf, m, 1)
模型残差分析bj2
表明预测误差没有任何信息-它不是自相关的,也与输入不相关。这表明Box-Jenkins模型的额外动态元素(C和D多项式)能够捕捉噪声动态。
渣油(zv bj2)
频率函数比较
为了比较生成模型的频率函数,我们再次使用bodeplot
命令:
CLF opt = bodeoptions;opt.PhaseMatching =“上”;W = linspace(0.1,4*pi,100);bodeplot (md2, GS, m, mx, mtf am2, bj2, w,选择);传奇({“温泉”,“党卫军”,“ARX”,“助教”,“OE”,“ARMAX”,“北京”},“位置”,“西方”);
还可以对估计模型的噪声谱进行分析。例如,我们在这里比较了噪声谱,ARMAX和Box-Jenkins模型与状态空间和频谱分析模型。为此,我们使用光谱
命令:
谱(GS, m, mx, am2 bj2, w)传说(“温泉”,“党卫军”,“ARX”,“ARMAX”,“北京”);
估计模型与真实系统的比较
在这里,我们用真实系统验证估计模型m0
我们用来模拟输入和输出数据。让我们将ARMAX模型的频率函数与真实系统进行比较。
波德(m0 am2 w)
这些反应几乎是一致的。噪声谱也可以进行比较。
谱(m0 am2 w)
让我们再看看零点点图:
H = iopzplot(m0,am2);
可以看出,真实系统(蓝色)和ARMAX模型(绿色)的极点和零点非常接近。
我们还可以计算零点和极点的不确定性。要绘制估计极点和零点周围对应3个标准差的置信区域,请使用showConfidence
或者从情节的上下文(右键单击)菜单中打开“置信区域”特性。
showConfidence (h, 3)
%
我们看到,真正的蓝色零点和极点都在绿色的不确定区域内。