系统描述

噪声破坏线性系统可描述为:

y = ug + H

在哪里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指的是“eXtra”输入B(q)u(t)。

用一般传递函数来表示G而且H，该模型对应一个参数化:

G (q) = B (q) / (q)而且H (q) = C (q) / (q)，有共同点

模拟模型

我们产生一个输入信号u并模拟模型对这些输入的响应。的idinput命令可用于创建模型和iddata命令将信号打包成合适的格式。的sim卡命令可以用来模拟输出，如下所示:

prevRng = 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));

只使用一部分数据进行评估是一种很好的做法，泽并保存另一部分以稍后验证估计的模型:

泽= z (1:200);zv = z (201:350);

进行光谱分析

现在我们已经得到了模拟数据，我们可以估计模型并进行比较。让我们从光谱分析开始。我们调用水疗中心命令，该命令返回频率函数和噪声频谱的频谱分析估计。

GS = spa(泽);

频谱分析的结果是频率的复值函数:频率响应。它被打包到一个叫做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 =轴;轴(10 ax (3:4) [0.1])

估计GS是不确定的，因为它是由噪声损坏的数据形成的。光谱分析方法也提供了它自己的这种不确定性的评估。要显示不确定度(例如3个标准差)，我们可以使用showConfidence命令在情节手柄上h返回的bodeplot命令。

showConfidence (h, 3)

图中显示，尽管频率响应(蓝线)的名义估计不一定准确，但有99.7%的概率(正态分布的3个标准差)，真实响应在阴影(浅蓝色)区域内。

估计参数状态空间模型

接下来让我们估计一个默认的(不指定特定的模型结构)线性模型。它将被计算为使用预测误差方法的状态空间模型。我们使用党卫军函数在本例中:

m = ss(泽)模型的顺序将自动选择。

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"。状态:使用时域数据m0上的sest进行估计。与估计数据拟合:75.08%(预测焦点)FPE: 0.9835, MSE: 0.9262

在这一点上，状态空间模型矩阵没有告诉我们太多。我们可以比较模型的频率响应米与光谱分析的估计结果相吻合GS．请注意,GS离散时间模型是while吗米是一个连续时间模型。它是可以使用的党卫军也可以计算离散时间模型。

bodeploy (m,GS) ax =轴;轴(10 ax (3:4) [0.1])

我们注意到这两个频率响应非常接近，尽管它们的估计方法非常不同。

简单传递函数的估计

线性模型结构种类繁多，都对应于描述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)2个零和2个极点的%传递函数

mtf =从输入“u1”输出“日元”:-0.05428 s ^ 2 - 0.02386 + 29.6  ------------------------------- s ^ 2 + 1.361 + 4.102连续时间确定传递函数。参数化:极点数:2零数:2自由系数数:5参数及其不确定性使用"tfdata"， "getpvec"， "getcov"。状态:使用时域数据“m0”上的TFEST进行估计。与估计数据拟合:71.69% FPE: 1.282, MSE: 1.195

ARX模型表示为:A(q)y(t) = B(q)u(t) + e(t)，

或者在火车,

这个模型的系数是线性的。系数 ${\mathit{一个}}_1^{}$， ${\mathit{一个}}_2^{}$、…… ${\mathit{b}}_1^{}$， ${\mathit{b}}_2^{}$, . .可以使用最小二乘估计技术有效地计算。对具有规定结构的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"表示参数及其不确定性。状态:使用ARX对时域数据“m0”进行估计。拟合估计数据:68.18%(预测焦点)FPE: 1.603, MSE: 1.509

估计模型的性能比较

这些模特是(mtf，mx)优于或劣于默认状态空间模型米？一种方法是用验证数据的输入来模拟每个模型(无噪声)zv并将相应的模拟输出与集合中的实测输出进行比较zv：

比较(zv, m, mtf, mx)

拟合度是由模型解释的输出变化的百分比。很明显米这是更好的模式吗mtf接近。离散时间传递函数可以用特遣部队或者是oe命令。前者创建IDTF模型，后者创建输出-误差结构的IDPOLY模型(y(t) = B(q)/F(q)u(t) + e(t))，但它们在数学上是等价的。

md1 =特遣部队(泽,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"。状态:使用时域数据“m0”上的TFEST进行估计。拟合估计数据: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模型具有与仿真模型相同的噪声特性m0Box-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"表示参数及其不确定性。状态:使用ARMAX对时域数据“m0”进行估计。与估计数据拟合: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 *π,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)

％

我们看到真正的蓝色零点和极点都在绿色的不确定区域内。