Matlab做状态空间辨识2_长笛人倚楼Gloria

http://blog.sina.com.cn/u/2214742602

首页博文目录关于我

个人资料

微博

加好友发纸条

写留言加关注

博客等级：
博客积分：

博客访问：
关注人气：
获赠金笔：0支
赠出金笔：0支
荣誉徽章：

正文字体大小：大中小

Matlab做状态空间辨识2

(2013-11-17 16:10:30)

标签：

ident

n4sid

ssest

系统辨识

状态空间

分类：信号处理

转载自了凡春秋USTC https://chunqiu.blog.ustc.edu.cn/?p=334

命令行辨识状态空间模型

准备工作：

构造iddata对象（输入输出数据对象）或frd、idfrd对象（频率响应数据对象）；

数据预处理，如对时域数据做去趋势；

选择阶次。

状态空间辨识有两种方式，决定于你对系统了解的先验知识。

黑箱辨识：

指定模型阶次，选择性的对一些模型矩阵结构属性进行配置。这种辨识方式，你可以使用ssest或n4sid函数以数据和模型阶次作为主要函数参数来辨识。额外的属性要以名字-值对的形式来指定，如指定模型模型采样时间、直通分量存在性、噪声成分是否存在等。你不能直接配置各系统矩阵的元素。

结构化辨识：

你需要首先创建一个包含系统矩阵初始值的idss模型结构，在此结构中，使用Structure属性来指定参数约束。如指定系统矩阵某个元素固定或设定其范围。详细参数化配置可以参考ssform函数说明。

在配置好约束后，就可以此idss模型作为ssest的输入进行辨识了。注意n4sid不能结构化辨识。

另外注意，这里的结构化辨识其实就是灰箱辨识，属于一种简单的一种，复杂的灰箱辨识可以使用idgrey和idnlgrey模型实现；结构化辨识中不能指定系统矩阵中各元素之间的关系，即各元素之间是独立的，对于元素之间有依赖性的结构化辨识属于复杂灰箱辨识，可尝试使用grayest估计器。

命令行辨识的实现：

使用ssest辨识

m = ssest(data,n,opt,Name,Value)

其中，data为估计的数据，n为系统阶次，opt包含了状态空间估计的配置参数，包括初始条件、输入偏移、估计中心、搜索方法等，其详细配置见ssest函数。

使用n4sid辨识

m = n4sid(data,n,opt,Name,Value)

除非指定采样时间，ssest默认辨识连续时间模型；而n4sid辨识离散模型。

连续模型与离散模型的选择：

默认情况下，ssest估计连续时间模型。如果使用非零采样周期的数据，可以使用如下命令估计离散模型

model = ssest(data,nx,'Ts',data.Ts);

如果使用的是连续时间频域数据，则不能辨识离散模型。

默认地，n4sid以数据的采样周期辨识离散时间模型，如果想辨识连续时间模型，可以用如下命令

model = n4sid(data,nx,'Ts',0);

对于状态空间辨识，你可以指定是否辨识 \(D\) 、 \(K\) 、 \(X0\) 矩阵，分别表示输入-输出直通分量、噪声模型和初始状态。

\(D\) ：

默认情况下， \(D\) 矩阵是不辨识的，除静态模型外其值固定为0。

黑箱辨识中，使用Feedthrough名和值对输入参数来设置直通分量的存在与否。如两输入模型中，第二输入直通，则使用如下命令设置

model = n4sid(data,n,'Feedthrough',[false true]);.

结构辨识中，通过配置init_sys.Structure.d来设置，其中init_sys是一个代表理想模型结构的idss结构。如强制第 \(i\) 个输入为直通，则

init_sys.Structure.d.Value(:,i) = 0;

init_sys.Structure.d.Free = true;

init_sys.Structure.d.Free(:,i) = false;

第一行指定矩阵 \(D\) 的第 \(i\) 列的值为0；第二行指定 \(D\) 矩阵的所有元素都是自由的，需辨识的；第三行指定 \(D\) 的第 \(i\) 列的值在辨识中是固定不变的。

另外，也可以再ssform中以参数对的形式配置直通。

\(K\) ：

\(K\) 代表系统的噪声矩阵，如模型的噪声成分

\[\dot x = Ax + Ke \\ y_n = Cx + e\]

对于频域数据，无噪声模型辨识， \(K\) 被设置为0。对于时域数据， \(K\) 默认以黑箱的设置辨识。

黑箱辨识中，使用DisturbanceModel参数对来表示噪声模型是设置为0（如果值设为'none'）还是以自由参数辨识（值设为'estimate'）。如

model = n4sid(data,n,'DisturbanceModel','none')

在结构辨识中，通过配置init_sys.Structure.k来设置。你可以固定一些 \(K\) 矩阵系数的值或设置一些值的上下界范围。如，只辨识一个两输出模型的 \(K\) 矩阵的第一列

kpar = init_sys.Structure.k;

kpar.Free(:,1) = true;

kpar.Free(:,2) = false;

kpar.Value(:,2) = 0; % 第二列的值固定为0

init_sys.Structure.k = kpar;

还可以使用ssform来配置。

当不确定如何配置 \(K\) 时，可以先不辨识它，只辨识系统矩阵。在得到系统模型后，再使用ssest配置 \(K\) 为自由参数辨识改善模型。如

init_sys = ssest(data, n,'DisturbanceModel','none');

init_sys.Structure.k.Free = true;

sys = ssest(data, init_sys);

init_sys为无噪声的动态模型。

要在一个已存在的模型中设置 \(K\) 为0，可以如下配置

m.Structure.k.Value = 0;

m.Structure.k.Free = false;

\(X0\) ：

初始状态向量 \(X0\) 一般作为模型辨识的副产品，即在n4sid和ssest函数中返回的第二个输出参数为 \(X0\) 。可以使用InitialState选项来选择如何控制初始状态。使用n4sid的n4sidOptions和ssest的ssestOptions来创建辨识选项。如在用n4sid辨识中保持初始状态为零

opt = n4sidOptions;

opt.InitialState = 'zero';

[m,X0] = n4sid(data,n,opt);

返回的 \(X0\) 为长度为 \(n\) 的零向量。

在使用ssest迭代辨识中，初始状态有Auto、Zero、Estimate和Backcast四个选项，具体说明见上面GUI中的介绍，在命令行中使用方法如下

opt = ssestOptions('InitialState','estimate');

m = ssest(data, 4, opt)

当使用多组实验数据辨识时， \(X0\) 矩阵的列数要和实验数据的一样。

辨识自由参数状态空间模型

自由辨识中，状态空间各矩阵 \(A\) 、 \(B\) 、 \(C\) 、 \(D\) 、 \(K\) 的参数是自由，即它们的所有元素在辨识过程中都是可调整的。由于 \(A\) 、 \(B\) 、 \(C\) 是自由的，所以状态空间实现的基是在满足最优条件下自动选择的。

如果对离散状态空间模型的内部结构知识了解不多，可以使用下面语句做快速尝试辨识：

m = ssest(data)

此语句辨识1到10阶的连续状态空间模型。

要辨识指定阶数的模型，使用如下语法

m = ssest(data,n)

迭代方法ssest是以n4sid的结果为初始模型的，所以你可以直接使用n4sid代替ssest来辨识：

m = n4sid(data,n)

辨识标准型状态空间模型

标准型辨识表示矩阵 \(A\) 、 \(B\) 、 \(C\) 是稀疏的，即除一些行或列的值可辨识外，其他元素固定为0或1。支持的标准型辨识如下：

伴随型：特征多项式的系数出现在 \(A\) 矩阵的最右列；

模型分解型： \(A\) 矩阵是块对角线的形式；

可观标准型：自由参数只出现在 \(A\) 的一些行。

对于模型分解型，其块对角化元素是对称的，如果使用结构化辨识，这些对称性将不会保留，虽然得到的模型仍然是块对角的。

命令行辨识标准型的语句如下

m = ssest(data, n, 'Form', 'canonical')

上面语句默认辨识连续时间模型，辨识离散模型使用如下语句

md = ssest(data, n, 'Form', 'canonical', 'Ts', data.Ts)

辨识结构参数状态空间模型

结构参数化辨识允许用户通过设置一些参数为固定值来剔除一些辨识参数。这种方法是很有用的，如果你的状态空间矩阵是从物理原理得来的，那么一些参数将有其物理意义，是固定的，这时可以使用此种方法做辨识。

结构化辨识需要两个步骤：

1.使用idss命令指定状态空间矩阵的结构和辨识参数的初始值；

2.使用ssest辨识未知参数。

要在GUI中实现结构化辨识，需要先在命令行中创建模型结构，然后倒入到GUI中辨识。

定义离散时间状态空间结构，语句如下

m = idss(A,B,C,D,K,'Ts',T)

A,B,C,D,K指定矩阵初始值，T指定采样间隔。

相似的，定义连续时间模型，语句如下

m = idss(A,B,C,D,K,'Ts',0)

在创建了模型结构后，你需要指定哪些参数是要辨识的，哪些参数使用设定值。这可以使用模型的Structure属性来设置，需要设置的详细属性如下

Value：参数值；

Minimum：辨识过程中参数的最小值；

Maximum：辨识过程中参数的最大值；

Free：布尔值，指定此参数是否要辨识。true为需要辨识，false表示固定为Value中值；

Scale：参数值的比例，它不在辨识中使用；

Info：保存参数的单位和标签信息，它是一个包含Lable和Unit的结构。

举个例子，假如你创建了一个如下A矩阵的状态空间模型m

A = [2 0; 0 3]

如果你想固定A(1,2) = A(2,1)=0，则使用如下语句

m.Structure.a.Value(1,2) = 0;

m.Structure.a.Value(2,1) = 0;

m.Structure.a.Free(1,2) = false;

m.Structure.a.Free(2,1) = false;

另外，快速配置参数、指定直通分量、干扰动态特性，可以使用ssform。

在结构化辨识中，对需要辨识的参数可以使用符合物理意义的值作为辨识初始值。但由于这可能导致迭代过程收敛到局部极小，所以建议尝试几组不同的初始化值辨识。如果想使用随机初始化值，可以使用init命令。另外，当模型结构包含不同量级的参数时，建议缩放变量使它们具有相近的量级。

迭代搜索算法中使用数值微分来计算预测误差对参数的梯度，步长的大小由nuderst命令，默认步长是参数绝对值*10的-4次方与10的-7次方的大者。如果想指定不同的步长，可通过nuderst来编辑。

辨识结构化离散时间状态空间模型的例子

这个例子中，我们辨识如下离散方程中的未知参数 \((\theta_1,\theta_2,\theta_3,\theta_4,\theta_5)\) ：

假设未知参数 \((\theta_1,\theta_2,\theta_3,\theta_4,\theta_5)\) 的名义值为-1,2,3,4,5。

辨识此离散模型的步骤如下：

1.构造参数矩阵，以参数名义值来初始化参数

A = [1,-1;0,1];

B = [2;3];

C = [1,0];

D = 0;

K = [4;5];

2.构造状态空间模型对象

m = idss(A,B,C,D,K);

3.指定不想辨识的参数

S = m.Structure;

S.a.Free(1,1) = false;

S.a.Free(2,:) = false;

S.c.Free = false;

m.Structure = S;

如果想固定初始状态为已知的0值，可以在InitialState选项中设置

opt = ssestOptions;

opt.InitialState = 'zero';

4.估计模型

m = ssest(data, m, opt)

data为包含时域数据或频域数据的iddata对象。迭代搜索以各系统矩阵的名义值开始搜索。

辨识结构化连续时间状态空间模型的例子

这个例子中，我们辨识如下连续时间模型中的未知参数 \((\theta_1,\theta_2,\theta_3)\)：

这个方程代表一个电机模型， \(y_1(t) = x_1(t)\) 是电机轴的角位置， \(y_2(t) = x_2(t)\)是角速度。参数 \(-\theta_1\) 是电机时间常数的逆， \(-\theta_2/\theta_1\) 是输入到角速度的静态增益。

电机在 \(t = 0\) 时是静止的，但其角度位置 \(\theta_3\) 是未知的。假设位置参数的近似名义值为 \(\theta_1 = -1\) 、 \(\theta_2 = 0.25\) 。位置测量的误差方差为0.01，角速度测量的方差为0.1。

辨识此连续模型的步骤如下：

1.构造参数矩阵、初始化参数名义值

A = [0 1;0 -1];

B = [0;0.25];

C = eye(2);

D = [0;0];

K = zeros(2,2);

x0 = [0;0];

2.构造连续时间状态空间模型对象

m = idss(A,B,C,D,K,'Ts',0);

3.指定在辨识中固定的值

S = m.Structure;

S.a.Free(1,:) = false;

S.a.Free(2,1) = false;

S.b.Free(1) = false;

S.c.Free = false;

S.d.Free = false;

S.k.Free = false;

m.Structure = S;

m.NoiseVariance = [0.01 0; 0 0.1];

初始状态是部分未知的，可以用InitialState选项来设置 \(X0\)

opt = ssestOptions;

opt.InitialState = idpar(x0);

opt.InitialState.Free(2) = false;

4.辨识模型

m = ssest(data, m, opt)

其中data是包含时域或频域数据的iddata对象。这个例子中可以使用dcmotordata.mat中的数据，它在matlabroot/toolbox/ident/iddemos/data/下。

5.如果想仿真的，可以构造采样周期为T=0.1的仿真数据，如下

e = randn(300,2);

u = idinput(300);

simdat = iddata([],u,'Ts',0.1);

simopt = simOptions('AddNoise', true, 'NoiseData', e)

y = sim(m,simdat,simopt)

通过辨识状态空间模型来得到等价的ARMAX和OE模型

你可以使用状态空间模型来辨识等价的ARMAX和OE多输出模型。对于ARMAX，需要指定状态空间辨识估计 \(K\) 矩阵，而对于OE模型，设置 \(K=0\) 。最后将得到的状态空间模型通过idpoly转化为相应的多项式格式。

下面的例子展示了如何利用状态空间方法辨识ARMAX和OE格式的模型。

首先加载测量数据，辨识状态空间模型

load iddata1 z1

mss_noK = n4sid(z1, 2,'DisturbanceModel','none');

mss = n4sid(z1,2);

mss_noK是二阶无干扰模型的状态空间模型，mss是有噪声成分的二阶状态空间模型。两者都使用z1数据来辨识。

将状态空间模型转换为多项式模型：

mOE = idpoly(mss_noK);

mARMAX = idpoly(mss);

上面的转化将是参数方差信息丢失。你可以重新计算方差：使用相同的估计数据做零迭代更新；使用translatecov做高斯估计公式。

使用零迭代更新方法重估计mOE和mARMAX的多项式模型参数如下

opt = polyestOptions;

opt.SearchOption.MaxIter = 0;

mOE = polyest(z1,mOE);

mARMAX = polyest(z1,mARMAX);

使用零迭代的方法，模型参数保持不变，只有方差信息得到更新。

除此之外，使用translatecov将估计的模型转化为多项式模型，代码如下

fcn = @(x)idpoly(x);

mOE = translatecov(fcn, mss_noK)

mARMAX = translatecov(fcn, mss)

因为polyest和translatecov使用不同的算法，所以得到的方差数据可能是不同的。使用如下命令查看模型参数的不确定性

present(mOE)

present(mARMAX)

另外，你也可以使用 \(K=0\) 的状态空间模型作为Hammerstein-Wiener的初始辨识，这可以提高模型的拟合效果，详见相关帮助。

来源： <https://chunqiu.blog.ustc.edu.cn/?p=334>

来自为知笔记(Wiz)

阅读┊ 收藏 ┊ 喜欢 ▼ ┊打印┊举报/Report

前一篇：Matlab做状态空间辨识1

后一篇：Matlab 灰箱辨识

新浪BLOG意见反馈留言板　欢迎批评指正