系统辨识

数字控制理论与设计上机报告

一、背景知识

辨识是通过测量得到的动态系统输入/输出数据确定系统的模型。设计和实现高性能的控制系统需要系统模型。多数实际情况中，我们需要能够直接辨识动态（控制）模型的方法。主要采用的是辨识离散时间参数化动态模型的方法。

二、实验目的。

1、学会使用递归最小二乘法（RLS）与递归增广最小二乘法（RELS）进行系统辨识。

2、理解上述两种方法的异同并验证辨识得到的模型，加深对数字控制理论的理解。

三、实验原理

1、系统辨识

本次实验采用的是Ⅰ型辨识方法——基于白化预测误差的辨识方法。

辨识采样离散时间参数化动态模型的方法是最适合设计和调试数字控制系统的方法。系统辨识是一种确定系统动态模型的实验方法，它包含有四个步骤：

⑴在设定的实验条件下采集输入/输出数据；

⑵选择或估计模型结构（复杂度）；

⑶估计模型参数；

⑷验证辨识模型（结构和参数值）。

一个完整的系统辨识算法必须包含以上四个步骤，每个步骤使用的方法取决于期望模型的类型（参数化的还是非参数化的、连续时间的还是离散时间的）以及实验条件（如对噪声的假设、开环或闭环辨识）。模型验证是判断辨识模型是否有效的关键一步。由于获取高质量辨识模型的参数估计算法和实验方法并非唯一，并且获得的模型也不是总能够通过验证测试。在这种情况下需要重新考虑辨识算法、模型复杂度和实验条件，系统辨识应视为一种反复过程。

2、验证模型

主要验证基于白化预测误差的方法辨识得到的模型，验证方法的准则如下：

⑴ 选择的对象+扰动结构是否正确，即能否代表实际系统； ⑵ 选择的多项式()()()1

1

1

A q

B q

C q ---、、的阶数和d （时延）的值是否正确。

验证方法基于以上准则，由以下几步构成：

⑴ 为辨识模型创建一个I/O 文件（使用和原系统一样的输入序列）； ⑵ 为辨识模型创建一个预测误差文件（最少100个采样值）； ⑶ 测试预测误差序列（即预测残差序列）的白度（不相关性）。

四、实验内容

1、Session LS1-1

打开T0.mat 文件，将已经测量得到的数据T0.mat 装载入工作空间。 输入plot(u)及plot(y)观察系统输入u 及输出y 的波形，如下图所示：

input u

然后输入以下指令对输入输出数据进行归一化处理： >> y=y-mean(y) >> u=u-mean(u)

因为用于辨识的模型结构应能表达动态关系，亦即在工作的附近，输入变化如何影响输出的函数关系。因此，为了辨识准确，应从I/O 数据中消除（对应工作点的）直流分量。分为两个阶段完成：①计算I/O 数据的平均值；②从I/O 数据中减去平均值，得到新的I/O 数据。因此需要对I/O 数据进行上述语句的处理。

打开函数文件夹下的rls.m 文件，依据Session LS1-1中所给的条件，输入 >> na=2;nb=2;d=1; >> [B,A]=rls(y,u,na,nb,d)

而rls.m 文件中函数的格式为[B,A]=rls(y,u,na,nb,d,Fin,lam1,lam0)，其中输入参数后三个这里都没有指定，因此舍去，使用程序中默认的缺省值。得到辨识出的结果如下：

B =

0 0 1.0000 0.5000 A =

1.0000 -1.5000 0.7000 也就是()1

2

310.5B q

q

q ---=+，()1121 1.50.7A q q q ---=-+，说明我们能够将此

output y

系统模型辨识出来。

通过同样的步骤再打开T1.mat 数据文件，与T0.mat 文件不同之处在于其是基于以下ARMAX 模型获取：

()()()()()()111d A q y t q B q u t C q e t ----=+

其中，多项式()1

A q

-和()1B q -与前面相同，()e t 为新增的随机扰动，是由计算机

产生的近似白噪声，()1

C q -的阶数为2C

n

=。

最后Matlab 运算得到的结果为： B =

0 0 0.9831 0.6512 A =

1.0000 -1.4035 0.6066 也就是()1

2

30.98310.6512B q

q

q ---=+，()1121 1.40350.6066A q q q ---=-+。与

上述辨识存在差异。这是由于文件T1受噪声的高度干扰，RLS 方法将产生有偏差的参数估计，因此需要对辨识结果的质量进行验证，同时验证与辨识都是基于同样的输入/输出数据集进行。

验证所用的函数[wlossf,ulossf,wrni,urni,wyhat,uyhat]=olvalid(B,A,C,y,u)封装在olvaid.m 文件中，其输入参数中带有C ，但是之前的辨识模型中并没有考虑C ，也就是并没有得到C ，因此我们默认C=1，输入以下指令：

>> [wlossf,wrni]=olvalid(B,A,1,y,u) 得到白噪声图像如下：

由图可知其归一化自相关函数()N R i 在1i =时超出了实践值的限制。也就是说估计的参数有偏差，验证结果也不理想（()10.15N R >），预测残差不近似白噪声，因此，需要考虑其他的对象+扰动结构。

下面将采用RELS 方法再次对数据T1.mat 进行系统辨识与验证。

打开rels.m 文件，通过函数名[B,A,C]=rels(y,u,na,nb,nc,d,Fin,lam1,lam0)可以发现输入比RLS 方法中多了nc ，输出多了个矩阵C 来形容干扰误差。在Matlab 命令行输入：

>> [B,A,C]=rels(y,u,na,nb,na,d)

输入中的后一个na 表示的是多项式C 的阶次nc 与na 相等，得到的模型为： B =

0 0 0.9516 0.5479 A =

1.0000 -1.4981 0.6942 C =

1.0000 -0.8383 0.0415

也就是()1230.95160.5479B q q q ---=+，()112

1 1.49810.6942A q q q ---=-+，

0.51 1.52

2.53

3.54

4.55

i

W R N (i )

Whiteness test for the OL identification

()11210.83830.0415C q q q ---=-+。此时我们辨识得到的模型则比较准确，再进行

白噪声模型验证可得：

此时从图形中能够直观的看出其归一化自相关函数()N R i 仅在1i =时超过了理论值，但是没有超出实践值的限制，因此认为此模型是准确可靠的。

以上在进行系统辨识前我们都假定已知了多项式()()()

1

1

1

A q

B q

C q ---、、的阶次以及时延d 大小，然而在实际操作中，我们并不一定能清楚的知道这几者的大小，因此需要对模型的阶数进行估计。

[N0,nA0,nB0,d0,Tab]=estorderiv(y,u,nmax,sel)，nmax 表示的是多项式中的最大次数，一般选取20，对于一般的实际系统这么高的次数足够了。对于输出Tab 表示的是各个阶次的可能性，同时()0max ,a b N n d n d =++。

>> [N0,nA0,nB0,d0,Tab]=estorderiv(y,u,20,0) 输入上述语句后得到仿真结果为：

N0 = 3 nA0 = 2 nB0 = 2

0.51 1.52

2.53

3.54

4.55

i

W R N (i )

Whiteness test for the OL identification

d0 = 1

同时我们可以通过观测图像()IV J N criterion ，可以看到在N=3时取得一个最小值。

然后我们可以通过这个求得的多项式阶次以及时延大小进行系统辨识。 >> na=nA0; nb=nB0; d=d0;

后面可以运用RLS 以及RELS 方法进行系统辨识便进行验证。

下面再来看下辨识出的系统响应情况。 系统传递函数如下： Transfer function: 0.9581 z^-2 + 0.542 z^-3 ------------------------------ 1 - 1.5 z^-1 + 0.696 z^-2 阶跃响应曲线：

Instrumental Variables estimation of system order (N)

Order N

J I V (N ) c r i t e r i o n

频率响应曲线：

同样的我们可以通过roots(A)、pzmap(model)等语句看其特征根以及零极点分布等情况。 2、Session LS1-2

根据I/O 数据直接辨识四个实际系统的参数,包括一个空气加热器、一个直流电机、一个挠性传动以及一个挠性机械臂。

051015

20253035

Step Response

Time (sec)

A m p l i t u d

e

-20-1001020

30M a g n i t u d e (d B )

-225

-180-135-90-45

0P h a s e (d e g )

Bode Diagram

Frequency (rad/sec)

1)aeroc: 128 I/O samples corresponding to an air heater. Input: voltage applied to

a power amplifier that feeds the heating resistor, Output: temperature of the

output air. The input signal is a PRBS( shift register with N = 6 and frequency divider p=2). The sampling period is 5 s.

在Matlab中打开数据文件aeroc.mat，将空气加热器数据装载入工作空间。然后输入以下指令对输入输出数据进行归一化处理：

>> y=y-mean(y)

>> u=u-mean(u)

打开函数文件夹下的estorderiv.m文件，运行后在command窗口输入

>> [N0,nA0,nB0,d0]=estorderiv(y,u,20,0)

得到的运行结果为：

N0 =

2

nA0 =

2

nB0 =

1

d0 =

J N criterion，可以看到在N=2时取得一个最小值。

观测图像()

IV

随后运用RLS方法进行系统辨识，得到的系统模型为：

B =

0 0.1778 A =

1.0000 -0.9210 0.1865 也就是()1

1

0.1778B q q

--=，()1

1

210.92100.1865A q

q

q ---=-+。进行白噪声模

型验证可得：

可以看出其归一化自相关函数()N R i 在i 取不同值时均未超出理论限制，因此可以认为此模型是可靠的。

输入>> model=filt(B,A,5)，得到系统传递函数如下： Transfer function: 0.1778 z^-1 --------------------------------- 1 - 0.921 z^-1 + 0.1865 z^-2 阶跃响应曲线如下：

频率响应曲线如下：

2)mot3c: 256 I/O samples corresponding to a dc motor. Input: voltage applied to a

power amplifier that feeds the motor, Output: speed of the motor measured by a tachometer. The input signal is a PRBS( shift register with N = 7 and frequency

divider p=2). The sampling period is 15 ms.

操作同上，打开数据文件mot3c.mat，将直流电机数据载入工作空间，进行归一化处理后调用[N0,nA0,nB0,d0,Tab]=estorderiv(y,u,20,0)对模型阶数进行估计：

N0 =

2

nA0 =

1

nB0 =

2

d0 =

J N criterion，可以看到在N=2时取得一个最小值。

观测图像()

IV

随后运用RLS方法进行系统辨识，得到的系统模型为：

B =

0 0.2629 0.2257

A =

1.0000 -0.5403

进行白噪声验证可得：

从图形中能够直观的看出其归一化自相关函数()N R i 在1i =和i=2时不仅超过了理论值，更是超出了实践值的限制。

也就是说估计的参数有偏差，验证结果也不理想（()10.15N R >），预测残差不近似白噪声，因此，需要考虑其他的对象+扰动结构。

选择na=nc=2，nb=2，d=0，输入>> [B,A,C]=rels(y,u,na,nb,na,d)，得到模型：

B =

0 0.2612 0.2319 A =

1.0000 -0.5270 -0.0070 C =

1.0000 0.4307 0.3491 进行白噪声验证可得：

R i<0.15，模型准确可靠。其归一化自相关函数()

N

系统传递函数如下：

Transfer function:

0.2612 z^-1 + 0.2319 z^-2

------------------------------------

1 - 0.527 z^-1 - 0.006981 z^-2

阶跃响应曲线如下：

频率响应曲线：

3)poulbo1c: 256 I/O samples corresponding to a flexible transmission (see course).

Input: position of the first pulley, Output: position of the third pulley. The input signal is a PRBS ( shift register with N = 7 and frequency divider p=2). The sampling frequency is 20 Hz. There are two vibration modes

操作同上，打开数据文件poulbo1c.mat，将挠性传动数据载入工作空间，进行归一化处理后调用[N0,nA0,nB0,d0,Tab]=estorderiv(y,u,20,0)对模型阶数进行估计：

N0 =

4

nA0 =

4

nB0 =

2

d0 =

2

()

J N criterion图像如下：

IV

进行白噪声验证可得：

从图形中能够直观的看出其归一化自相关函数()N R i 在1i =和i=2时不仅超过了理论值，更是超出了实践值的限制。

也就是说估计的参数有偏差，验证结果也不理想（()10.15N R >），预测残差不近似白噪声，因此，需要考虑其他的对象+扰动结构。

选择na=nc=4，nb=2，d=2，输入>> [B,A,C]=rels(y,u,na,nb,na,d)，得到模型： B =

0 0 0 0.3027 0.4033 A =

1.0000 -1.6016 1.8654 -1.4929 0.8865 C =

1.0000 -0.7136 0.2605 0.0568 -0.0897 进行白噪声验证可得：

R i<0.15，所得模型准确可靠。其归一化自相关函数()

N

系统传递函数如下：

Transfer function:

0.3027 z^-3 + 0.4033 z^-4

-----------------------------------------------------------------

1 - 1.60

2 z^-1 + 1.865 z^-2 - 1.49

3 z^-3 + 0.8865 z^-4

频率响应曲线：

Question 1: at what frequencies are the vibration modes? What are the corresponding damping factors?

4)256 I/O samples corresponding to a very flexible robot arm. Input: axis position.

Output: robot arm end position. The input signal is a PRBS. The sampling frequency is 5 Hz. Any vibration mode over 2.5 Hz has been eliminated by an anti-aliasing filter. A priori knowledge: in the frequency region from 0 to 2.5 Hz, there are two vibration modes. The file has to be centered.

操作同上，打开数据文件Rob2.mat，将挠性机械臂数据载入工作空间，进行归一化处理后调用[N0,nA0,nB0,d0,Tab]=estorderiv(y,u,20,0)对模型阶数进行估计：

N0 =

4

nA0 =

4

nB0 =

4

d0 =

随后运用RLS 方法进行系统辨识，得到的系统模型为：

B =

0 0.0255 0.2980 0.4481 0.1325 A =

1.0000 0.0820 -0.7391 0.2491 0.7887 进行白噪声验证可得：

从图形中能够直观的看出其归一化自相关函数()N R i 在1i =、i=3及i=4时不仅超过了理论值，更是超出了实践值的限制。

也就是说估计的参数有偏差，验证结果也不理想（()10.15N R >），预测残差不近似白噪声，因此，需要考虑其他的对象+扰动结构。

选择na=nc=4，nb=4，d=0，输入>> [B,A,C]=rels(y,u,na,nb,na,d)，得到模型：

B =

0 0.0263 0.2962 0.4934 0.1704 A =

1.0000 0.2408 -0.8311 0.1691 0.9330 C =

1.0000 0.6203 -0.1426 0.2804 0.5055