系统辨识课程报告

z(k ) (k ) v(k )
（2.2）
上式中， (k ) 是观测数据向量。辨识的问题则是如何求 a i 和 bi 。当模型结构已经选 定，即式（2.2）中的 na 和 nb 假定是已知的（通常 na〉nb） 。但是有时，为了系统分析方 便起见，也可以设 na=nb=n，这样并不会失去研究问题的普片性。 在（2.2）式中，有
（3.5）
将上式代人（3.4）式中得
2 L log L( ML ) const log v 2
（3.6）

故
AIC (n) 2const L log v 4 n
2

2
（参数个数=2n）
（3.7）
其等价式为
AIC (n) L log v 4 n







AIC (n) 中起主导作用。因此 AIC (n) 在 n0 处形成了一个最小值，如图 3.1 所示。
3.2 AIC 法确定线性定常系统的阶 考虑如下模型


z (k ) ai z (k i) bi u (k i) v(k )
i 1 i 1
n
n
(3.1)
伪随机序列可由线性移位寄存器网络产生。该网络由 r 级串联的双态器件，移位脉冲 产生器和模 2 加法器组成，下面以 4 级移位寄存器为例，说明伪随机序列的产生。例如， 初始状态是 0001，那么 X1=0，X2=0，X3=0，X4=1。如果反馈逻辑为 X= X3⊕ X4，对于 初始状态为 0001， 经过一个时钟节拍后， 各级状态自左向右移到下一级， 未级输出一位数， 与此同时模 2 加法器输出值加到移位寄存器第一级，从而形成移位寄存器的新状态，下一 个时钟节拍到来又继续上述过程。未级输出序列就是伪随机序列。其产生的伪随机序列为 X=100110101111000100110101111000… ，这是一个周期为 15 的周期序列。改变反馈 逻辑的位置及数量还可以得到更多不同的序列输出。如图 1 所示。


参数估计量 满足
lT ( Z l l ) 0

（2.9）
这就是以向量矩阵形式表示的正规方程， 它由具有 2l+1 个线性方程式所组成。 可解出
(lT l ) 1lT Z l
2.3 最小二乘估计的递推算法归结
T l 1 l K l 1 [ z l 1 l 1 l ] Pl l 1 K l 1 1 lT1 Pl l 1 Pl 1 [ I K L 1 lT1 ]Pl
2 其中，u(k) ,z(k)为系统的输入输出变量，v(k)是服从 N (0, v ) 分布的不相关随机噪
声，上述模型可写成最小二乘矩阵格式
z n H n n Vn
则有
（3.2）
L 2
L( n ) (2 )
2 v
1 T exp ( z H ) ( z H ) n n n n n n 2 2 v
参考文献
杨承志，孙棣华，张长胜.系统辨识与自适应控制.2003 年 7 月第一版.重庆大学出版社.
《系统辨识》课程报告
学号：2007073124
附：MATLAB 程序
1 最小二乘一次完成算法 M序列子函数：
function M=U X1=1;X2=0;X3=1;X4=0;X5=1;X6=0; %移位寄 存器输入Xi初态（101010），Yi为移位寄存器的 各级输出 m=300; %置M序列的总长度 for i=1:m Y6=X6; Y5=X5; Y4=X4; Y3=X3; Y2=X2; Y1=X1; X6=Y5; X5=Y4;X4=Y3; X3=Y2; X2=Y1; X1=xor(Y5,Y6); %异或运算 if Y6==0 U(i)=-1; else U(i)=Y6; end end M=U; J1=0; for x=2:300 J1=J1+[z(x)+z(x-1)*c1(1)-u(x-1)*c1(2 )].^2; end JJ1=J1+z(1).^2; AIC(1)=30*log(JJ1/300)+4*1; %若为二阶 z=zeros(1,300); for k=4:300 z(k)=1.5*z(k-1)-0.7*z(k-2)-0.1*z(k-3 )+1.1*u(k-1)+1.5*u(k-2)+1.7*u(k-3); %用理想输出值作为观测值 end for i=3:300 HL2(i,:)=[-z(i-1) -z(i-2) u(i-1) u(i-2)]; ZL2(i,:)=[z(i)]; end c2=inv(HL2'*HL2)*(HL2'*ZL2); J2=0; for x=3:300 J2=J2+[z(x)+z(x-1)*c2(1)+z(x-2)*c2(2 )-u(x-1)*c2(3)-u(x-2)*c2(4)].^2; end JJ2=J2+z(2).^2+z(1).^2; AIC(2)=30*log(JJ2/300)+4*2; %若为三阶 z=zeros(1,300); for k=4:300 z(k)=1.5*z(k-1)-0.7*z(k-2)-0.1*z(k-3 )+1.1*u(k-1)+1.5*u(k-2)+1.7*u(k-3); %用理想输出值作为观测值 end for i=4:300 HL3(i,:)=[-z(i-1) -z(i-2) -z(i-3) u(i-1) u(i-2) u(i-3)]; ZL3(i,:)=[z(i)]; end c3=inv(HL3'*HL3)*(HL3'*ZL3); J3=0; for x=4:300 J3=J3+[z(x)+z(x-1)*c3(1)+z(x-2)*c3(2 )+z(x-3)*c3(3)-u(x-1)*c3(4)-u(x-2)*c 3(5)-u(x-3)*c3(6)].^2; end JJ3=J3+z(3).^2+z(2).^2+z(1).^2; AIC(3)=30*log(JJ3/300)+4*3; %若为四阶
（3.3）
故
log L( n )
L L 1 log 2 log v2 ( z n H n n ) T ( z n H n n ) 2 2 2 v2
（3.4）
式中，L 为数据长度。根据 ML 原理，有
T 1 T ML ( H n H n ) H n z n 2 1 T v ( z n H n ML ) ( z n H n ML ) L

（2.10）
（2.11）
3 根据信息的准则估计模型的阶次
1974 年， Akaike 提出了一个基于 Kullback-Leible 信息量判据来定阶的赤池信息准则 AIC（Akaike Information Criterion） 。 3.1 AIC 准则意义
AIC (n)

n0
图 3.1
n
AIC (n) -- n 图
利用数据序列{z(k)}和{u(k)},极小化下列准则函数
l
T
J z (k ) (k )
T k 1


2
（2.6）

使 J(θ )=min 的θ 估计值记作 ，称之为参数的最小二乘估计值。 2.2 最小二乘问题的求解 用表示根据 l 次数采样数据所求得的参数的估计值 。

《系统辨识》课程报告
及 AIC (n) 值，找到使 AIC (n) min 的 n 作为 n 0 。


（3.8）
具体的定阶用法是：对不同阶次首先使用极大似然法估计参数，然后计算似然函数值
4 例子及其 MATLAB 程序实现
待辨识系统的阶次大于等于三阶；输入采用六位 M 序列，且用子函数实现；数据长度 大于等于三百；包含系统阶次的辨识。
被辨识的系统
图 2.1
SISO 辨识系统示意图
《系统辨识》课程报告
学号：2007073124
假设被辨识系统为一个单输入单输出（SISO）离散时间动态系统，如图 2.1 所示。其 系统数学关系式可用如下随机差分方程描述
z (k ) ai z (k i) bi u (k i) v(k )
其中
（2.5）
Z l [ z (1), z (2),, z (l )]T Vl [v(1), v(2),, v(l )]T
u (0) u (1 nb ) (1) z (0) z (1 na ) (2) z (1) z (2 n ) u (1) u (1 nb ) a l (l ) z (l 1) z (l na ) u (l 1) u (l nb )
《系统辨识》课程报告
学号：2007073124
待辨识系统为：
z(k ) 1.5 * z(k 1) 0.7 * z(k 2) 0.1* z(k 3) 1.1* u(k 1) 1.5 * u(k 2) 1.7 * u(k 3)
实验结果： c = -1.5000 0.7000 0.1000 1.1000 1.5000 1.7000 jc = 3 结果分析： 实验结果和理论值是完全相等的。由于实验中的观测值是使用的理想观测值，并没有 引入误差因素。系统的阶次辨识中，运用 AIC 信息准则，在求得的几个系数中比较出最小 值，从而得出系统阶次。
i 1 i 1
na
nb
k 1,2,3,, l
(2.1)
式中
u(k)—输入变量； y(k)—系统输出变量； z(k)—系统量测输出变量； v(k)—表示均值为零的随机噪声项。