系统辨识实验报告
一、相关分析法
(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 ][
残差为： )(?)()(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. 限定记忆法（固定窗法）
程序如下：