使用系统识别工具箱构建结构化和用户定义的模型
这个例子展示了如何估计用户定义的模型结构中的参数。这种结构由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 = ax + bu + K e y = cx + du + e
我们将做一个直流电机的模型。电机的动力学是众所周知的。如果我们选择x1作为角位置，x2作为角速度，很容易建立一个忽略干扰的状态空间模型:(见Ljung(1999)中的例4.1:
| 0 1 | | 0 | d/dt x = | | x + | | u | 0 -th1 | | th2 |
| 10 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 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”上使用SSEST估计。拟合估计数据:[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;dcmodel2.Structure.a.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”参数及其不确定性。状态:在时域数据“z”上使用PEM估计。拟合估计数据:[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 |
| 10 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) %返回dc电机的状态空间矩阵，其中%时间常数Tau (Tau = par)和已知静态增益G。样本%时间为Ts。% %如果输入参数Ts %为零，该文件返回连续时间表示。如果Ts>0，则返回离散时间表示。要使使用此文件的IDGREY模型意识到这种灵活性，请设置Structure的%值。FcnType属性到'cd'。这种灵活性对于估计和模拟所需的连续域和离散域之间的转换非常有用。参见IDGREY, IDDEMO7。版权所有:The MathWorks, Inc. t = par(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.25赋给辅助变量G(增益)和采样时间。
Aux = 0.25;DCMM = idgrey(“motorDynamics”par_guess,“cd”辅助0);
时间常数现在由
dcmm =灰色(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”表示。状态:在时域数据“z”上使用GREYEST估计。拟合估计数据:[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元素表示为:
地点:
因此它是一个多项式1 /问的程度nakl．
类似地，B(q)是一个ny | nu矩阵，其kj-元为:
因此有一个延迟nkkj从输入数字j输出数字k．最常见的创建方法是使用arx命令。订单被指定为:Nn = [na nb nk]与na作为一个ny-by-ny矩阵的kj进口是nakj；注而且nk定义相似。
让我们用数据测试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"。状态:在时域数据“z”上使用ARX估计。拟合估计数据:[97.87;83.44]%(预测焦点)FPE: 0.002157, MSE: 0.1398
结果,dcarx1，被存储为IDPOLY模型，并且应用所有前面的命令。例如，我们可以显式地列出arx多项式:
dcarx1.a
Ans = 2x2单元阵列{[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)u(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"。状态:在时域数据“z”上使用ARX估计。拟合估计数据:[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模型为快速估计具有用户指定结构的多输出模型提供了另一种选择。
您也可以从以下列表中选择一个网站:
