使用系统识别工具箱构建结构化和用户定义的模型

打开脚本

这个例子展示了如何在用户定义的模型结构中估计参数。这种结构由IDGREY(线性状态空间)或IDNLGREY(非线性状态空间)模型指定。我们将研究如何分配结构，固定参数和创建它们之间的依赖关系。

实验数据

我们将研究由(模拟的)直流电机产生的数据。我们首先加载数据:

负载dcmdatayu

矩阵y包含两个输出:日元电机轴的角度位置和y2是角速度。有400个数据样本，采样时间为0.1秒。输入包含在向量中u．它是电机的输入电压。

Z = iddata(y,u,0.1);IDDATA对象z.InputName =“电压”；z.OutputName = {“角”；“AngVel”};情节(z (: 1:))

图:测量数据:电压对角度

情节(z (: 2:))

图:测量数据:电压对角速度

模型结构选择

d/dt x = A x + B u + K e y = C x + du + e

我们将制作直流电机的模型。电机的动力学是众所周知的。如果我们选择x1作为角位置，x2作为角速度，很容易建立一个忽略干扰的状态空间模型，其性质如下:(见Ljung(1999)中的例4.1:

| 0 1 | | | d/dt x = | | x + | | u | 0 -th1 | | th2 |

| 1 0 | y = | | x | 0 1 |

的参数th1这是电机的逆时间常数吗th2是这样的th1 / th2从输入到角速度的静态增益。(详见Ljung(1987))th1而且th2与电机的物理参数有关)。我们将根据观测到的数据估计这两个参数。上面描述的模型结构(参数化状态空间)可以在MATLAB®中使用IDSS和IDGREY模型表示。这些模型允许您使用实验数据进行参数估计。

标称(初始)模型的规范

如果我们猜测th1=1和th2 = 0.28，我们得到名义或初始模型

A = [0 1;0 1];A(2,2)的初始猜测为-1B = [0;0.28);%对B(2)的初始猜测为0.28C =眼睛(2);D = 0 (2,1);

我们将其打包到一个IDSS模型对象中:

ms = idss(A,B,C,D);

该模型的特征是矩阵、它们的值、哪些元素是自由的(待估计的)以及它们的上限和下限:

ms.Structure.a

ans =名称:'A'值:[2x2 double]最小值:[2x2 double]最大值:[2x2 double]空闲值:[2x2 logical]比例:[2x2 double]信息:[2x2 struct] 1x1 param。连续

ms.Structure.a.Value ms.Structure.a.Free

Ans = 0 1 0 -1 Ans = 2x2逻辑阵列1 1 1 1

使用IDSS模型的自由(独立)参数规范

所以我们现在应该标记，只有A(2,2)和B(2,1)是需要估计的自由参数。

ms.Structure.a.Free = [0 0;0 1];ms.Structure.b.Free = [0;1);ms.Structure.c.Free = 0;%使用标量展开ms.Structure.d.Free = 0;ms.Ts = 0;这将模型定义为连续的

初始模型

女士初始型号

=连续时间状态空间模型确定女士:dx / dt = x (t) + B u e (t) + K (t) y (t) = C x (t) + D u (t) + e (t) = (x1, x2) x1 0 1 x2 0 1 B = u1 x1 0 x2 0.28 C y2 = (x1, x2)日元1 0 0 1 D = u1 y1 y2 0 K = y₁y2 x1 0 0 x2 0 0参数化:结构化形式(一些固定系数A、B、C)。引线:没有干扰组件:没有很多免费的系数:2使用“idssdata”、“getpvec”、“getcov”参数及其不确定性。现状:直接建造或改造而成。不估计。

IDSS模型自由参数的估计

现在计算参数的预测误差(最大似然)估计为:

dcmodel = sest(z,ms,ssestOptions)“显示”，“上”));dcmodel

dcmodel =连续时间状态空间模型发现:dx / dt = x (t) + B u e (t) + K (t) y (t) = C x (t) + D u (t) + e (t) = (x1, x2) x1 0 1 x2 0 x1 0 x2 1.002 C = -4.013 B =电压(x1, x2)角1 0 AngVel 0 1 D =电压角0 AngVel 0 K =角AngVel x1 0 0 x2 0 0参数化:结构化形式(一些固定系数A、B、C)。引线:没有干扰组件:没有很多免费的系数:2使用“idssdata”、“getpvec”、“getcov”参数及其不确定性。状态:在时域数据“z”上使用sest进行估计。拟合估计数据:[98.35;84.42]% FPE: 0.001071, MSE: 0.1192

参数的估价值与模拟数据时的值(-4和1)非常接近。为了评估模型的质量，我们可以通过与实际输入的模型进行模拟，并将其与实际输出进行比较。

比较(z, dcmodel);

例如，我们现在可以画出零点和极点以及它们的不确定区域。我们将画3个标准差对应的区域，因为模型是相当准确的。注意，原点的极点是绝对确定的，因为它是模型结构的一部分;从角速度到位置的积分器。

clf showConfidence (iopzplot (dcmodel), 3)

现在，我们可以做一些修改。a矩阵的1,2元素(固定为1)告诉我们这一点x2是的导数x1．假设传感器没有校准，所以可能有一个未知的比例常数。为了包括这样一个常数的估计，我们只是“放松”(1、2)和重新评估:

Dcmodel2 = dcmodel;free (1,2) = 1;Dcmodel2 = pem(z, Dcmodel2);

得到的模型是

dcmodel2

dcmodel2 =连续时间状态空间模型发现:dx / dt = x (t) + B u e (t) + K (t) y (t) = C x (t) + D u (t) + e (t) = (x1, x2) x1 0.9975 0 x2 0 x1 0 x2 1.004 C = -4.011 B =电压(x1, x2)角1 0 AngVel 0 1 D =电压角0 AngVel 0 K =角AngVel x1 0 0 x2 0 0参数化:结构化形式(一些固定系数A、B、C)。引线:没有干扰组件:没有很多免费的系数:3使用“idssdata”、“getpvec”、“getcov”参数及其不确定性。状态:用PEM对时域数据“z”进行估计。拟合估计数据:[98.35;84.42]% FPE: 0.001077, MSE: 0.1192

我们发现估计(1、2)接近于1。为了比较这两个模型，我们使用比较命令:

比较(z, dcmodel dcmodel2)

使用IDGREY对象的耦合参数模型规范

假设我们准确地知道直流电机的静态增益(从输入电压到角速度，例如从先前的阶跃响应实验。如果静态增益为G，电机的时间常数为t，则状态空间模型为

1 | 0 | | 0 | d / dt x = | | x + | | u | 0 1 / t | | G / t |

|1 0| y = | | x |0 1|

与G已知，不同矩阵的项之间是有相关性的。为了描述这一点，先前使用的带有“Free”参数的方法是不够的。因此，我们必须编写一个MATLAB文件，它生成a、B、C和D，以及可选的K和X0矩阵作为输出，对于每个给定的参数向量作为输入。它还将辅助参数作为输入，因此用户可以更改模型结构中的某些内容，而不必编辑文件。在这种情况下，我们让已知的静态增益G作为这样一个论点进入。已写入的文件名为“motorDynamics.m”。

类型motorDynamics

function [A,B,C,D,K,X0] = motorDynamics(par,ts,aux) % motorDynamics表示电机动力学的ODE文件。% % [A,B,C,D,K,X0] = motorDynamics(Tau,Ts,G) %返回具有%时间常数Tau (Tau = par)和已知静态增益G的直流电机的状态空间矩阵。样本%时间为Ts。% %如果输入参数Ts %为零，该文件返回连续时间表示。如果Ts>0，则返回离散时间表示。要使使用此文件的IDGREY模型意识到这种灵活性，请设置Structure的%值。FcnType属性设置为'cd'。这种灵活性对于估计和模拟所需的连续和离散域之间的转换非常有用。参见IDGREY, IDDEMO7。版权所有The MathWorks, Inc. (1);G = aux(1);A = [0 1;0 -1/t];B = [0;G/t]; C = eye(2); D = [0;0]; K = zeros(2); X0 = [0;0]; if ts>0 % Sample the model with sample time Ts s = expm([[A B]*ts; zeros(1,3)]); A = s(1:2,1:2); B = s(1:2,3); end

现在，我们创建了一个对应于此模型结构的IDGREY模型对象:假设的时间常数将为

Par_guess = 1;

我们还给辅助变量赋值0.25G(增益)和采样时间。

Aux = 0.25;DCMM = idgrey(“motorDynamics”par_guess,“cd”辅助0);

时间常数现在由

dcmm = greyest(z,dcmm,greyestOptions)“显示”，“上”));

这样我们就直接估计出了电机的时间常数。它的值与先前的估计很吻合。

dcmm

dcmm =连续时间线性灰箱模型定义为“motorDynamics”功能:dx / dt = x (t) + B u e (t) + K (t) y (t) = C x (t) + D u (t) + e (t) = (x1, x2) x1 0 1 x2 0 x1 0 x2 1.001 C = -4.006 B =电压(x1, x2)角1 0 AngVel 0 1 D =电压角0 AngVel 0 K =角AngVel x1 0 0 x2 0 0模型参数:Par1 = 0.2496参数化:颂歌功能:motorDynamics(用参数表示连续和离散方程)干扰组件:初始状态:ODE函数参数化自由系数数:1参数及其不确定性用"getpvec"， "getcov"表示。状态:使用GREYEST对时域数据“z”进行估计。拟合估计数据:[98.35;84.42]% FPE: 0.00107, MSE: 0.1193

有了这个模型，我们现在可以像以前一样继续测试各个方面。所有命令的语法都与前一种情况相同。例如，我们可以将idgrey模型与其他状态空间模型进行比较:

比较(z, dcmm dcmodel)

他们显然非常亲密。

多元ARX模型的估计

工具箱的状态空间部分还处理多变量(多个输出)ARX模型。通过多元arx模型，我们的意思是:

A(q) y(t) = B(q) u(t) + e(t)

这里A(q)是一个ny |ny矩阵，它的项是延迟算符1/q中的多项式。k-l元素表示为:

$ $现代{kl} (q) $ $

地点:

$a_{kl}(q) = 1 + a_1 q^{-1} + ....+ a_{nakl} q^{-nakl} q$

因此它是一个多项式1 /问的程度nakl．

同样，B(q)是一个ny | nu矩阵，其kj元为:

$$b_{kj}(q) = b_0 q^{-nkk}+b_1 q^{- nkj -1}+…+ b_ {nbkj}问^ {-nkkj-nbkj} $$

因此有一个延迟nkkj从输入数字j输出数字k．创建这些文件的最常用方法是使用arx命令。订单指定为:Nn = [na nb nk]与na作为一个ny-by-ny矩阵的kj进口是nakj；注而且nk定义相似。

让我们用dc数据测试一些arx模型。首先，我们可以简单地建立一个二阶模型:

Dcarx1 = arx(z，“na”, 2, 2, 2, 2,“注”(2, 2),“朝鲜”(1, 1))

dcarx1 =离散ARX模型:模型输出“角”:一个(z) y_1 (t) = - ai (z) y_i (t) + B (z) u (t) + e_1 (t) (z) = 1 - 0.5545 z ^ 1 - 0.4454 z ^ 2 A₂(z) = -0.03548 z ^ 1 - 0.06405 z ^ 2 B (z) = 0.004243 z ^ 1 + 0.006589 z ^ 2模型输出“AngVel”:一个(z) y_2 (t) = - ai (z) y_i (t) + B (z) u (t) + e_2 (t) (z) = 1 - 0.2005 z ^ 1 - 0.2924 z ^ 2 A_1 (z) = 0.01849 z ^ 1 - 0.01937 z ^ 2 B (z) = 0.08642 z ^ 1 + 0.03877 z ^ 2样品时间:0.1秒参数化:多项式订单:na=[2 2;2 2] nb=[2;2] nk=[1;1]自由系数数:12参数及其不确定性使用"polydata"， "getpvec"， "getcov"。状态:使用ARX对时域数据“z”进行估计。拟合估计数据:[97.87;83.44]%(预测焦点)FPE: 0.002157, MSE: 0.1398

结果,dcarx1，存储为IDPOLY模型，并且应用前面的所有命令。例如，我们可以通过以下方法显式列出arx多项式:

dcarx1.a

Ans = 2x2 cell array {[1 -0.5545 -0.4454]} {[0 -0.0355 -0.0640]} {[0 0.0185 -0.0194]} {[1 -0.2005 -0.2924]}

作为单元格数组，例如dcarx1的{1,2}元素。a是之前描述过的多项式a(1,2)关于y2和y1。

我们也可以测试一个我们知道的结构日元是通过过滤得到的y2通过一阶滤波器。(角度是角速度的积分)。我们还可以假设从输入到输出2的一阶动力学:

Na = [1 1;0 1];Nb = [0;1);Nk = [1;1);Dcarx2 = arx(z，[na nb nk])

dcarx2 =离散时间ARX模型:输出“Angle”的模型:A(z)y_1(t) = - A_i(z)y_i(t) + B(z)u(t) + e_1(t) A(z) = 1 - 0.9992 z^-1 A_2(z) = -0.09595 z^-1 B(z) = -1 B(t) + e_2(t) A(z) = 1 - 0.6254 z^-1 B(z) = 0.08973 z^-1采样时间:0.1秒参数化:多项式阶:na=[1 1;0 1] nb=[0;1] nk=[1;1]自由系数:4使用"polydata"， "getpvec"， "getcov"表示参数及其不确定性。状态:使用ARX对时域数据“z”进行估计。拟合估计数据:[97.52;81.46]%(预测焦点)FPE: 0.003452, MSE: 0.177

为了比较得到的不同模型，我们使用

比较(z, dcmodel dcmm dcarx2)

最后，我们可以对不同模型的输入和输出的bodeploy进行比较波德:第一次输出:

Dcmm2 = idss(dcmm);%转换为IDSS用于子引用波德(dcmodel (1, 1),“r”dcmm2 (1,1),“b”dcarx2 (1,1),‘g’）

第二个输出:

波德(dcmodel (2, 1),“r”dcmm2 (1),“b”dcarx2 (1),‘g’）

前两个模型或多或少是完全一致的。由于非白方程误差噪声引起的偏差，arx模型不太好。(我们在模拟中有白测量噪声)。

结论

预选结构的模型估计可以使用系统识别工具箱执行。在状态空间形式中，参数可以固定为已知值，也可以限制在规定的范围内。如果需要指定参数或其他约束之间的关系，则可以使用IDGREY对象。IDGREY模型评估用户指定的MATLAB文件，用于估计状态空间系统参数。多变量ARX模型为快速估计具有用户指定结构的多输出模型提供了另一种选择。