zy系统辨识

《系统辨识》课程报告
题目：Matlab实现系统的阶次及模型辨识班级：工业自动化081班
姓名：吴政
学号：2008073121
日期：2011-5-28
成都信息工程学院控制工程学院
Matlab实现系统的阶次及模型辨识
摘要
基于最小二乘原理确定系统模型阶次n，常用的定阶方法有六种：残差方差定阶，AIC 准则，残差白色定阶，零极点消去定阶，行列式比定阶，Hankel矩阵定阶。

在已知系统阶次的情况下，为了辨识出系统的模型，方法有很多种，比如一般最小二乘法，遗忘因子最小二乘法，增广最小二乘法，广义最小二乘法等等。

本课程报告根据残差平方和（作为损失函数）达到最小值的直接法来辨识系统的阶次，然后，通过递推最小二乘法来辨识出系统的模型，整个过程在Matlab环境下完成。

关键词：系统辨识；最小二乘法； Matlab；高斯白噪声
Using Matlab To Identify The Order And Model Of System
Abstract
In order to determine the order of model n based on the Least Squares Method(LSM), There are six commonly used method: order identification with the residual variance, AIC criterion, residual white, Pole-Zero cancellation, the ratio of the determinant, Hankel matrix,In the context of known the system’s order, There are many ways to determine the system’s model, such as general LSM, forgetting factor LSM, extended LSM, generalized LSM,etc. The course report based on the direct method when the residual sum of squares(as a loss function) achieve the minimum to identify the order of the system,then, identify the system model with recursive LSM,.The whole process completed in the Matlab environment.
1Key words:System Identification；Least Squares Method(LSM)；Matlab；
Gaussian White Noise
引言
1.1系统辨识概述
系统辨识的实质就是从一组模型类中选择一个模型，使之能够很好的拟合所关心的实际过程的动态和静态特性。

经典的辨识方法有：阶跃响应法，脉冲响应法，频率相应法，相关分析法，谱分析法，最小二乘法和极大似然法，其中最小二乘法是一种经典的和最基本的，也是应用最广泛的方法，可用于动静态。

线性，非线性系统。

1.2 递推最小二乘概述
最小二乘一次完成算法进行系统参数辨识的时候，存在一定的限定条件，并且需要用到全部的观测数据，每采样一次就需要增添一组新的观察数据，所以引入递推最小二乘法来辨识系统参数，递推最小二乘法是用旧的估计值加上修正值得到的新的估计值，用新的测量数据对上一次的估计结果进行修正，直到估计值达到需要的精度为止。

在阶次已知的前提下，我们可以通过最小二乘法的集中算法都可以确定模型参数，但是实际上我们对阶次并不知道。

实践表明，如果假定的模型阶次不合适，则系统的数学模型就不准确，进而就会在控制系统的分析和设计中产生不良结果。

1.3 课程报告要求
对本课程报告的要求具体如下： 1、待辨识系统阶次大于等于3； 2、输入采用6位M 序列；
3、采样数据长度大于等于300；
4、包含系统的阶次辨识。

2课程原理
2.1 M 序列
M 序列除了具有伪随机二位式序列的性质外，还具有如下一些特点： ① 一个n 级移位寄存器产生的M 序列的周期长度12-=n
p N 。

② 在M 序列的一个周期中，逻辑0状态总比逻辑1状态少1个（必须避免全为逻辑0的组
合）。

所以在一个周期中，由于逻辑1和逻辑0的状态数为p N ，则逻辑1状态出现的次数为2/)1(+p N ，逻辑0状态出现的次数为2/)1(-p N 。

③ 序列中某种状态连续出现的段称为“游程”长度。

一个n 级移位寄存器产生的M 序列，
在一个周期中有1
2
-n 段“游程”。

④ 所有M 序列都具有“移位相加”的性质。

这就是说，如果将一个M 序列与将它延迟了r
为以后的序列按模2加法原则相加，则所有的新序列是延迟了q 位的原来那个M 序列。

r 和q 都是整数，而且处在下列范围内：1≤r ， q ≤1-p N 。

M 序列产生方法：工程上产生M 序列采用移位寄存器法，下图1所示为四位移位寄存器产生M 序列的方法。

图1 四位移位寄存器产生M 序列
2.2 最小二乘递推算法辨识系统参数
根据最小二乘原理，用l 次观测数据，得出参数向量θ的最小二乘估计l
θˆ ()
Z
H H H
l
T
l
l
T l
1
-∧
=
θ……..②
其中，l
θˆ表示根据l 次观测数据所得到的最小二乘值计量，下表l 表示该符号代表l 次观测数据构成的矩阵。

T l l z z z Z )](),...,2(),1([=………③
l H =⎥⎥⎥⎥
⎥⎦
⎤⎢⎢⎢⎢⎢⎣⎡--------------)(.....)1()(.....)1(.
.)2(.....)1()2(.....)1()1(.....)0()
1(.....)0(n l u l u n l z l z n u u n z z n u u n z z …….④
令1
)(-=l T l l H H P ，且l P
是一方阵，它的维数取决于未知数的个数，而与观测次数无关
…⑤，
式中1+l h 表示第1+l 次观测数据。

利用矩阵反演公式计算⑤式
[][]l
T l l l
T l l l l
l T l l l P h h P h I h P P h h P P 111
11
1
1
111+-+++-++-++-=+=…….⑥，
式中[]1
1+++l l T
l h P h I 是一个标量，其求逆只是一个简单的除法。

令[]11
11
-++++=l l T
l l h P h I γ，于是有
[]l T
l l l l l P h h P I P 1
111
++++-=γ，则最小二乘估计量l
θˆ的表达式可改写为
l
T l l l z h P 1ˆ+=θ，于是第1+l 次观测数据所得到的最小二乘值计量为 [
]
11111111ˆˆ+++++++++==l l T l l l T l l l z h P P Z h P θθ，由⑤式可知T l l T l T l h h P P
111+++-=……⑦，于是 ，[
]
l T l l l l l h z K θθθˆˆˆ1111++++-+=……⑧，其中111+++=l l l h P
K ……⑨， 方程⑧就是最小二乘估计的递推算法。

综上所述，最小二乘递推算法归结如下：
[]
[]
⎪⎪
⎩
⎪
⎪⎨⎧-=+=-+=+++++++++++l T
l l l l l T
l l l l l T l l l l l P h K I P h P h h P K h z K 111111111111ˆˆˆθθθ 2.3 损失函数法原理辨识系统阶次
2.3.1 定阶原理
计算不同阶次n 辨识结果的估计误差方差，按估计误差方差最小或最显著变化原则来确定模型阶次n 。

常用的方法有按估计误差方差最小定阶和F 检验法（工程中采用F 检验法）。

2.3.2 残差方差最小定阶步骤
系统差分方程为：
)()()()()(1
1k k u z b k y z a ε+=-- 向量形式： E H Z +=θ LS 估计：
Z
H H H T
T 1^
)(-=θ
残差： )
()()(^
^
k z k z k e -=
指标函数：
∑++==L n n k n k e L J 1
2
^)
(1 依次计算n=1,2,3,···时的指标函数Jn ，并将其绘制成曲线。

如图2所示：
图2 n
J 曲线
定阶原则：则随着n 增大，J 值是下降的。

若n0为正确的阶次，此时J 值所在的点是曲线上最大的拐点，此后J 值基本不变化或变化很小。

依上述原则，上述曲线模型阶次为3。

2.3.3 F 检验法
实际工程应用时，在定阶过程中，我们并不是取Jn 最小时n 值，作为系统模型的阶次，而是对在n 增大过程中，使Jn 显著减小的n 值感兴趣。

为了避免人为主观判断的影响，引入准则：
选取
表1 F 与阶次的变化关系
依F 检验法，系统模型阶次为3
3实验步骤及Matlab 程序流程
3.1实验步骤
第一步：工作间清零后，利用6为移位寄存器产生长度为L （L ≥300）的M 序列，作为未知系统的输入信号。

第二步：用高斯白噪声作为待辨识系统的噪声输入，待辨识系统的阶次为3阶，故初始化时，要给)3(~)1(z z 赋初值。

第三步：利用循环语句产生长度为L+2的采样输出信号。

第四步：利用采样输出数据，求出5~1阶下的残差平方和。

第五步：将不同阶次下的残差平方和画在一张图中，找出拐点最大，此后值基本不变化或变化很小的点。

该拐点所对应的阶次就是系统的阶次。

第六步：在已知系统阶次的情况下，通过递推算法即可辨识出系统参数。

第七步：通过图形观察递推算法得到的历史数据曲线。

3.2 Matlab 程序流程
)n 2(n 12n N J J J F i 1i 1i 1i 1i i 1i ---⋅-=+++++
4实验程序及结果
4.1实验程序代码
%利用残差平方和求系统的阶次
clear all;
clc; %---工作间清空%=========产生M序列===========
x=[0 1 0 1 1 0]; %----initial value n=5000; % ---脉冲数目
u=[]; %---存放M序列
for i=1:n
temp=xor(x(3),x(6));
if(x(6)==0)
u(i)=-1;
else
u(i)=1;
end
for j=6:-1:2
x(j)=x(j-1);
end
x(1)=temp;
end
%=========产生高斯白噪声======
v=sqrt(0.01)*randn(1,5000); %----方差为0.01.均值为0
%=====产生观测序列============
z=[];
z(1)=-1;
z(2)=0;
z(3)=0;
L=400;
for i=4:L+5
z(i)=1.0*z(i-1)-0.5*z(i-2)+0.2*z(i-3)+0.9*u(i-1)-0.7*u(i-2)+1.2*u(i-3)+v(i-2); end
%========模型阶次n=1=====
H1=zeros(L,2);
for i=1:L
H1(i,1)=-z(i);
H1(i,2)=u(i);
end
estimate=inv(H1'*H1)*H1'*(z(2:L+1))';
e=(z(2:L+1))'-H1*estimate;
V1=e'*e/L;
%========模型阶次n=2=====
H2=zeros(L,4);
for i=1:L
H2(i,1)=-z(i+1);
H2(i,2)=-z(i);
H2(i,3)=u(i+1);
H2(i,4)=u(i);
end
estimate=inv(H2'*H2)*H2'*(z(3:L+2))';
e=(z(3:L+2))'-H2*estimate;
V2=e'*e/L;
%========模型阶次n=3=====
H3=zeros(L,6);
for i=1:L
H3(i,1)=-z(i+2);
H3(i,2)=-z(i+1);
H3(i,3)=-z(i);
H3(i,4)=u(i+2);
H3(i,5)=u(i+1);
H3(i,6)=u(i);
end
estimate=inv(H3'*H3)*H3'*(z(4:L+3))'; e=(z(4:L+3))'-H3*estimate;
V3=e'*e/L;
%========模型阶次n=4=====
H4=zeros(L,8);
for i=1:L
H4(i,1)=-z(i+3);
H4(i,2)=-z(i+2);
H4(i,3)=-z(i+1);
H4(i,4)=-z(i);
H4(i,5)=u(i+3);
H4(i,6)=u(i+2);
H4(i,7)=u(i+1);
H4(i,8)=u(i);
end
estimate=inv(H4'*H4)*H4'*(z(5:L+4))'; e=(z(5:L+4))'-H4*estimate;
V4=e'*e/L;
%========模型阶次n=5=====
H5=zeros(L,10);
for i=1:L
H5(i,1)=-z(i+4);
H5(i,2)=-z(i+3);
H5(i,3)=-z(i+2);
H5(i,4)=-z(i+1);
H5(i,5)=-z(i);
H5(i,6)=u(i+4);
H5(i,7)=u(i+3);
H5(i,8)=u(i+2);
H5(i,9)=u(i+1);
H5(i,10)=u(i);
end
estimate=inv(H5'*H5)*H5'*(z(6:L+5))'; e=(z(6:L+5))'-H5*estimate;
V5=e'*e/L;
%====递推最小二乘法辨识系统=====
%==========初始化===========
P=10*eye(6); %----估计方差
Pstore=zeros(6,L+1);
Pstore(:,1)=[P(1,1),P(2,2),P(3,3),P(4,4),P(5,5),P(6,6)];
Theta=zeros(6,L+1); %----参数估计值，存放中间过程估值Theta(:,1)=[0.5;0.5;0.5;0.5;0.5;0.5];
K=zeros(6,L); %----增益矩阵
K=[10;10;10;10;10;10];
E=0.000005;
%==========递推求解=========
for i=4:L+3
h=[-z(i-1);-z(i-2);-z(i-3);u(i-1);u(i-2);u(i-3)];%因为hi'=(hi1,hi2,hi3,hi4,...hin) K=P*h*inv(1+h'*P*h);
Theta(:,i-2)=Theta(:,i-3)+K*(z(i)-h'*Theta(:,i-3));
P=(eye(6)-K*h')*P;
Pstore(:,i-2)=[P(1,1),P(2,2),P(3,3),P(4,4),P(5,5),P(6,6)];
e1=Theta(:,i-2)-Theta(:,i-3); %----参数当前值与前一次值的差值
e2=e1./Theta(:,i-3); %----求参数的相对变化
e3(:,i)=e2; %----把当前相对变化的列向量加入误差矩阵的最后一列
if(abs(e2)<=E)
N=i;
break; %----若参数收敛满足要求，终止计算end
end
%=======画图=======
figure(1); %----第一个图形
stem(u),grid on; %----以径的形式显示出输入信号并给图形加上网格
figure(2)
plot(1:5,[V1 V2 V3 V4 V5]);
title('利用残差的方差辨识模型的阶次');
xlabel('阶次');
ylabel('残差方差');
figure(3)
i=1:L+1;
plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:),i,Theta(5,:),i,Theta(6,:));
title('Parameter Identification with Recursive Least Squares Method');%-- 递推最小二乘法辨识参数
figure(4)
i=1:L+1;
plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:),i,Pstore(5,:),i,Pstore(6,:))
title('估计方差变化过程');
figure(5)
i=1:N;
plot(i,e3(1,:),i,e3(2,:),':',i,e3(3,:),i,e3(4,:),':',i,e3(5,:),i,e3(6,:),':'); %画出a1，a2，a3,b1，b2,b3的各次辨识结果的收敛情况。

title('Identification Precision'); %----图形标题(辨识精度)
Theta(:,N-2)
4.2实验结果
图3 M序列
Fig.3 M Series
图4 利用残差的方差辨识模型的阶次
Fig.4 Order identification with residual sum of squares
表2 残差平方和的值
Fig.5 Parameter identification with recursive least squares method
表3 递推参数与实际参数的比较
Fig.6 Identification precision
5实验结果分析
不同阶次下残差平方和的值如下：
V1 =
0.7520
V2 =
0.3895
V3 =
0.0102
V4 =
0.0102
V5 =
0.0102
由图4及表2以及上述值可知：系统的阶次为3阶。

由图5及表3可知被辨识系统的参数分别为：
a1=
-1.0006
a2=
0.5091
a3=
- 0.1955
b1=
0.8930;
b2=-
0. 7007
b3=
1.2054
将上述值与真实值进行比较，可以发现递推算法建立的模型可以很好的拟合真实的系统。

参考文献
[1] 杨承志等. 系统辨识与自适应控制[M]. 重庆：重庆大学出版社.2003：62-64,106-107
[2] 王琳马平系统辨识方法综述[J].电力情报，2001（4）：63-66。