使用系统识别工具箱™构建结构化和用户定义模型
这个例子展示了如何在用户定义的模型结构中估计参数。这种结构由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 + d u + e
我们将做一个直流电机的模型。这马达的动力学是众所周知的。如果我们选择x1作为角位置，x2作为角速度，很容易建立一个忽略干扰的状态空间模型:(见Ljung(1999)的例子4.1:
| 0 1 | | 0 | 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模型对象中:
= idss女士(A, B, C, D);
该模型的特征是其矩阵、其值、哪些元素是自由的(待估计的)以及这些元素的上下限:
ms.Structure.a
ans = Name: 'A' Value: [2x2 double] Minimum: [2x2 double] Maximum: [2x2 double] Free: [2x2 logical] Scale: [2x2 double] Info: [2x2 struct] 1x1 param. Name: [2x2 struct] 1x1 param. Name: [2x2 double]连续
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. struct .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)。引线:没有干扰组件:none自由系数个数:2使用“idssdata”、“getpvec”、“getcov”表示参数及其不确定性。状态:由直接构造或转换产生。不估计。
IDSS模型自由参数的估计
参数的预测误差(最大似然)估计现在由:
dcmodel = ss (z,女士,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)。引线:参数及其不确定性使用“idssdata”、“getpvec”、“getcov”。状态:使用SSEST对时域数据“z”进行估计。拟合估计数据:[98.35;84.42]% FPE: 0.001071, MSE: 0.1192
这些参数的估价值非常接近模拟数据时使用的值(-4和1)。为了评估模型的质量，我们可以用实际输入的by模拟模型，并与实际输出进行比较。
比较(z, dcmodel);
例如，我们现在可以画出零点和极点以及它们的不确定区域。我们将画出3个标准差对应的区域，因为模型比较准确。注意，原点的极点是绝对确定的，因为它是模型结构的一部分;从角速度到位置的积分。
clf showConfidence (iopzplot (dcmodel), 3)
现在，我们可以做各种修改。a矩阵(固定为1)的1,2元告诉我们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”来表示参数及其不确定性。状态:使用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(τ,Ts, G) %返回的直流电机的状态空间矩阵%时间常数τ(τ= par)和已知的静态增益G .示例Ts %时间。% %这个文件返回连续时间表示如果输入参数Ts %是零。如果Ts>0，则返回一个离散时间表示。要使%使用此文件的IDGREY模型意识到这种灵活性，请设置Structure的%值。FcnType属性为'cd'。这种灵活性对于估算和模拟所需的连续域和离散域之间的转换非常有用。% %参见IDGREY, IDDEMO7。版权所有1986-2015 The MathWorks, Inc. t = par(1);G =辅助(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(增益)和采样时间。
辅助= 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参数化:ODE函数:motorDynamics(参数化连续和离散时间方程)扰动分量:由ODE函数参数化初始状态:由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元素表示为:
地点:
因此它是一个多项式1 /问的程度nakl．
同理，B(q)是一个ny |矩阵，其kj元素为:
因此有一个延迟nkkj从输入数j输出数量k．最常用的创建方法是使用arx命令。订单具体说明如下:n = [na nb nk]与na作为一个ny-by-ny矩阵的kj进口是nakj；注和nk类似的定义。
我们用dc数据测试一下arx模型。首先，我们可以简单地建立一个一般的二阶模型:
dcarx1 = arx (z,“na”, 2, 2, 2, 2,“注”(2, 2),“朝鲜”(1, 1))
A(z)y_1(t) = - A_i(z)y_i(t) + B(z)u(t) + e_1(t) A(z) = 1 - 0.5545 z^-1 - 0.4454 z^-2 A_2(z) = -0.03548 z^-1 - 0.06405 z^-2 B(z) = 0.004243 z^-1 + 0.006589 z^-2A(z)y_2(t) = - A_i(z)y_i(t) + B(z)u(t) + e_2(t) A(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通过一阶滤波器。(角是角速度的积分)。然后我们还可以假设从输入到输出的一阶动态:
Na = [1 1;0 1];Nb = [0;1);Nk = [1;1);Dcarx2 = arx(z，[na nb nk])
dcarx2 =离散ARX模型:模型输出“角”:一个(z) y_1 (t) = - ai (z) y_i (t) + B (z) u (t) + e_1 (t) (z) = 1 - 0.9992 z ^ 1 A₂(z) = -0.09595 z ^ 1 B (z) = 0模型输出“AngVel”:一个(z) y_2 B (t) = (z) u (t) + e_2 (t) (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)
最后，我们可以比较不同模型的输入和输出的bodeplots波德:首先输出:
dcmm2 = ids (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模型为使用用户指定的结构快速估计多输出模型提供了另一种选择。
你也可以从以下列表中选择一个网站:
