系统辨识实验报告

一、相关分析法
(1)实验原理
图1 实验原理图
本实验的原理图如图1。

过程传递函数()G s 中12120,8.3, 6.2K T Sec T Sec ===；输入变量()u k ，输出变量()z k ，噪声服从2(0,)v N σ，0()g k 为过程的脉冲响应理论
值，ˆ()g k 为过程脉冲响应估计值，()g k %为过程脉冲响应估计误差。

过程输入()u k 采用M 序列，其输出数据加白噪声()v k 得到输出数据()z k 。

利
用相关分析法估计出过程的脉冲响应值ˆ()g k ，并与过程脉冲响应理论值0()g k 比较，得到过程脉冲响应估计误差值()g k %。

M 序列阶次选择说明：首先粗略估计系统的过渡过程时间T S (通过简单阶跃响应)、截止频率f M (给系统施加不同周期的正弦信号或方波信号，观察输出)。

本次为验证试验，已知系统模型，经计算Hz T T f M 14.01
2
1≈=，s T S 30≈。

根据式M
f t 3
.0≤
∆及式S T t N ≥∆-)1(，则t ∆取值为1，此时31≥N ，由于t ∆与N 选择时要求完全覆盖，则选择六阶M 移位寄存器，即N =63。

(2)编程说明
图2 程序流程图
(3)分步说明 ① 生成M 序列：
M 序列的循环周期63126=-=N ，时钟节拍1t Sec ∆=，幅度1a =，移位寄存器中第5、6位的内容按“模二相加”，反馈到第一位作为输入。

其中初始数据设为{1,0,1,0,0,0}。

程序如下：
过程仿真
得到理论输出数据()z k
计算脉冲响应估计值
计算互相关函数，得到脉冲响应估计值
计算脉冲响应估计误差
计算脉冲响应理论值，得到脉冲响应估
计误差
人机对话
噪声标准差：sigma ；生成数据周期数：r
生成数据
生成M 序列()u k ；生成白噪声序列()v k
② 生成白噪声序列： 程序如下：
③ 过程仿真得到输出数据：
如图2所示的过程传递函数串联，可以写成形如1212
11
()
1/1/K G s TT s T s T =
++，
其中112
K
K TT =。

图2 过程仿真方框图
程序如下：
④ 计算脉冲响应估计值：
互相关函数采用公式)()(1
)(10
k i y i x N
r k R N r i xy +⋅⋅=
∑-⋅=，互相关函数所用的数据
是从第二个周期开始的，其中r 为周期数，取1-3之间。

则脉冲响应估计值为：
])([1
)(ˆc k R k k g xy -=，t
a N N k ∆⋅⋅+=2
)1(。

补偿量)1(-=N R c xy 。

程序如下：
⑤ 计算脉冲响应估计值： 脉冲响应的理论值由式12//012
()[
]k t T k t T K
g k e e T T -∆-∆=
--可计算得到。

这时可得到过程脉冲相应估计误差0ˆ()()()g
k g k g k =-%。

脉冲响应估计误差为： ∑∑===
N
k N
k g k g k g
1
2
01
2
))(())(~(δ
程序如下：
(4)数据记录
①当噪声标准差sigma=0.1，生成数据周期r 为2时：脉冲响应估计误差为0.0121。

脉冲响应估计曲线为图3所示。

②当噪声标准差sigma=0.5，生成数据周期r 为1时：脉冲响应估计误差为 0.0347。

脉冲响应估计曲线为图4所示。

③当噪声标准差sigma=0.5，生成数据周期r 为3时：脉冲响应估计误差为 0.0258。

脉冲响应估计曲线为图5所示。

④当噪声标准差sigma=1，生成数据周期r 为3时：脉冲响应估计误差为 0.0279。

脉冲响应估计曲线为图6所示。

图3 sigma=0.1，r=2时脉冲响应估计曲线图4 sigma=0.5，r=1时脉冲响应估计曲线
图5 sigma=0.5，r=3时脉冲响应估计曲线图6 sigma=1，r=3时脉冲响应估计曲线
(5)结果分析
实验中可以看到脉冲响应估计的曲线与理论曲线的重合度还是比较高的，脉冲响应估计误差也比较小，实验证明相关分析法的估计效果还是不错的。

同时，经过实验可以得出结论：固定数据周期r，给定不同的噪声标准差sigma可以发现，噪声的方差越大，也就是信噪比越大，估计的效果越不好；固定噪声标准差sigma，选择不同的数据生成周期r可以发现，数据周期越大，估计的周期越多，估计的效果越好。

二、最小二乘法
1. 基本最小二乘（离线辨识）
ζφθζ+=+⎥⎥⎥⎦
⎤
⎢⎢⎢⎣⎡=b a U Y Y ΛM ][
残差为： )(ˆ)()(k y
k y k e -= 最小二乘目标：残差平方和最小（一阶导为0，二阶导>0）。

[][]
正定，且φφφφφθT T T LS
Y 1ˆ-=⇒ 从上式看出，逆存在才有解，满足条件的u(k)：(1) 伪随机；(2) 白噪声；(3)有色随机信号。

程序如下：
结果如下：
result1 =[ -0.8287; 0.1275; -0.0024; 1.9884; -1.2723]
2. 递推最小二乘（在线辨识）
RLS )
1()]()([)()
()1()(1)
()1()()]1(ˆ)()()[()1(ˆ)(ˆ--=-+-=
--+-=k P k k K I k P k k P k k k P k K k k k y k K k k T T T φφφφθφθθ
修正项老估值新估值+=+N
N θθˆ)(ˆ)(1 为了启动RLS ，需给初值：I c P 20,0ˆ==θ。

计算框图见书P66。

程序如下：
结果如下：
result2 =[ -0.8284; 0.1274; -0.0024; 1.9883; -1.2717]
图7 递推最小二乘法参数过渡过程
数据饱和：（1）原因：01→+N K ，不再起修正作用，引起误差变大。

（2）为了克服数据饱和现象，可以用降低老的数据影响的方法：
①渐消记忆法（遗忘因子法）
②限定记忆法（固定窗法）
当)(k ζ为不相关序列，最小二乘有一致性与无偏性，但)(k ζ往往为相关序列，为克服最小二乘有偏估计的缺点，引入辅助变量法和广义最小二乘法，增广最小二乘法等。

3. 渐消记忆法（遗忘因子法）
α
φ
φ
φ
α
φ
θ
φ
θ
θ
)1
(
)]
(
)
(
[
)
(
)
(
)1
(
)
(
)
(
)1
(
)
(
)]
1
(ˆ)
(
)
(
)[
(
)1
(ˆ
)
(ˆ
-
-
=
-
+
-
=
-
-
+
-
=
k
P
k
k
K
I
k
P
k
k
P
k
k
k
P
k
K
k
k
k
y
k
K
k
k
T
T
T
一般99
.0
~
95
.0
=
α
程序如下：
结果如下：
result3 =[ -0.6862; 0.1114; 0.0640; 2.0356; -0.9429]
图8 渐消记忆法参数过渡过程
4. 限定记忆法（固定窗法）
程序如下：
结果如下：（较前三种方法偏差较大）
result4 =[ -1.0190; 0.1725; -0.4531; 2.1092; -1.3724]
5. 辅助变量法（IV ） (1).辅助变量Z
Y Z Z Z Z Y Z Z Z Z Y Z Z Q Q Z E Z E T T T T T T T T T T T 111)(ˆ)()()(噪声不相关
Z ,0)(---=-=⇒+===φθ
ζ
φφθζφθφφζ其中，辅助变量估计值由强烈相关与为非奇异阵，，与
(2).计算步骤：
①先根据实测数据最小二乘求粗略θ
ˆ（为有偏估计） ②)(ˆ)()(ˆ)(ˆ)(ˆ11k y k u z b k y z a
，得到辅助变量代入--= ③n )()(ˆ)(),(Z k u k y
k u k y 构造，用构造用φ ④n T n T Y Z Z =θφ (3).递推：
RIV )
1()]()([)()
()1()()
()1()()]1(ˆ)()()[()1(ˆ)(ˆ--=-+-=
--+-=k P k z k K I k P k z k P k k z k P k K k k k y k K k k T T T φαθφθθ
初始条件：I c P 20,0ˆ==θ
缺点: P0的选择非常敏感，一个改进方法是，用递推最小二乘辨识算法作为
启动方法，然后转换到辅助变量法。

程序如下：
结果如下：
result5 =[ -0.9369; 0.1207; -0.0254; 1.9781; -1.5331]
图9 辅助变量法参数过渡过程
6. 广义最小二乘法（GLS ）
(1)广义最小二乘法的基本思想：
由于)(k ξ在n+k 个采样周期的时差范围内具有自相关性，从而使θ的最小二乘估计为有偏的，所以引入一个所谓成形滤波器（白化滤波器），把相关噪声)(k ξ转化成白噪声)(k ε。

如果知道有色噪声序列的相关性：)()()()(11k z d k z c εξ--= 令)()()(111
---=z d z c z f ，有)()(1)(1k z f k εξ-= 广义最小二乘法(GLS)是建立在最小二乘法(LS)的基础上的。

基本最小二乘法只是广义最小二乘法在1)(1=-z f 时的特例。

(2)广义最小二乘法计算步骤：
广义最小二乘法的关键问题是如何用比较简单的方法找到成形滤波器的系数。

其计算是逐次逼近法。

①应用输入输出数据按最初模型求出θ的最小二乘估计。

这个估计值是不精确的，它只是被估参数的一次近似。

②计算残差e (k )，并拟合成形滤波器的模型：
)(ˆ)1(ˆ)(ˆ)(ˆ)1(ˆ)(y )(101n k y b k u b k u b n k y a k y a k k e n
n -------++-+=ΛΛ 得到 T m
T T f f f e f ]ˆˆˆ[][ˆ211Λ=ΩΩΩ=- 其中 ⎥⎥⎦
⎤⎢⎢⎣⎡+--+-+--=Ω)()1()1()(N m n e N n e m n e n e ΛM M Λ ③应用所得的成形滤波器，对输入输出数据滤波：
⎩⎨⎧-++-+=-++-+=)
(ˆ)1(ˆ)()()(ˆ)1(ˆ)()(11m k u f k u f k u k u m k y f k y f k y k y m m ΛΛ 其中，m 为噪声模型的阶，一般事先不知道，实际经验表明指定m 为2或3可以得到比较满意的输出。

④按新的输入、输出模型求出参数θ的第二次估计值θ
ˆ。

结果如下：
Result6 =[ -0.6538; 0.2926; 0.0454; 1.8776; -1.6042]
7. 广义递推最小二乘法（GLS ）
广义最小二乘法的递推计算过程可分成两个部分：
(1)按递推最小二乘法(RLS)，随着N 的增大，不断计算N
θˆ(逐步接近于无偏)
和N
f ˆ(逐步使噪声白化)； (2)在递推过程中，N θˆ和N
f ˆ是时变的，则过滤信号)(),(k y k u 及残差)(k e
是由时变系统产生，要不断计算。

因而，递推广义最小二乘法由两组普通的递推最小二乘法组成，它们是通过滤波算法联系起来的：
⎪⎪⎪⎩⎪⎪⎪⎨⎧ΦΦ=+-=+=-+=-+++++++++++++++1)1(1)1(1)1(11)1(1)1()1(11)1(11)1(1)1(111)1(11][11)ˆ(ˆˆT N N N T N N T N N N N N N N T N N N N N T N N N N N P P P P P P P P K y K ϕϕϕϕϕϕϕθϕθθ ⎪⎪⎪⎩⎪⎪⎪⎨⎧ΩΩ=+-=+=-+=-+++++++++++++++1)2(1)2(1)2(11)2(1)2()2(11)2(11)2(1)2(111)2(11]
[11)ˆ(ˆˆT N N N T N N T N N N N N N N T N N N N N T N N N N N P P P P P P P P K f e K f f ωωωωωωωω 结果如下：
result7 =[ -1.0164; 0.1754; -0.0159; 2.0056; -1.6438] 噪声传递系数的估计结果：[-0.0307; 0.0900]
图10 广义递推最小二乘参数过渡过程
8.增广矩阵法(ELS/RELS)(增广最小二乘法)
增广矩阵法是把观测矩阵适当增大，使得有偏估计的程度得到一定改善。

]
)()1()()()()1([N N n N u N n u N y N n y T N εεφΛΛΛ-++--+-= 这一方程结构适用于递推最小二乘法，但向量T N φ中)(k ε是未知的。

解决这
个矛盾的一个方法是用估计值)(ˆk ε
代替)(k ε，即： ])(ˆ)1(ˆ)()()()1([ˆN N n N u N n u N y N n y T N
ε
εφΛΛΛ-++--+-=其中，1ˆˆ)()(ˆ--+=+N T N N n y N n θφε，)(ˆN n +ε初始值可以取为零。

按照与前面类似的递推方法，得到增广矩阵法参数估计。

增广矩阵法在实际中有广泛的应用，收敛情况也比较好。

结果如下：
result8 =[ -0.8636; 0.0747; -0.0162; 1.9729; -1.1420]
图10 增广最小二乘参数过渡过程。