系统辨识最小二乘法—课设报告.doc

课 程 设 计 报 告

学 院：

专业名称：

学生姓名：

指导教师：

时 间：

课程设计任务书
一、设计内容
SISO系统的差分方程为：

参数取真值为：，利用MATLAB 的M语言辨识系统中的未知参数、、、。
二、主要技术要求
用参数的真值及差分方程求出作为测量值，是均值为0，方差为0.1、0.5和0.01的不相关随机序列。选取一种最小二乘算法利用MATLAB的M语言辨识参数。
三、进度要求
2周（6月28日-7月11日）完成设计任务，撰写设计报告3000字以上，应包含设计过程、 计算结果、 图表等内容。具体进度安排：
u 6月28日，选好题目，查阅系统辨识相关最小二乘法原理的资料。
u 6月29日，掌握最小二乘原理，用MATLAB编程实现最小二乘一次完成算法。
u 6月30日，掌握以最小二乘算法为基础的广义最小二乘递推算法。
u 7月1日，用MATLAB编程实现广义最小二乘递推算法。
u 7月2日，针对题目要求进行参数辨识，并记录观察相关数据。
u 7月3日-7月5日，对参数辨识结果进行分析，找出存在的问题，提出改进方案，验证改进优化结果。
u 7月6日-7月7日，撰写课程设计报告。
u 7月8日，对课程设计报告进行校对。
u 7月9日，打印出报告上交。
学 生
王景
指导教师
邢小军

本科课程设计报告

1. 设计内容
设SISO系统的差分方程为：
式（1-1）
参数取真值为：，利用MATLAB 的M语言辨识系统中的未知参数、、、。
2. 设计过程
2.1 问题重述。
设SISO系统的差分方程为：
式（2-1）
参数取真值为：，利用MATLAB 的M语言辨识系统中的未知参数、、、。
要求：用参数的真值利用差分方程求出作为测量值，是均值为0，方差为0.1、0.5和0.01的不相关随机序列。选取一种最小二乘算法辨识。
2.2 最小二乘参数辨识
2.2.1、 最小二乘法的概念与应用
对工程实践中测得的数据进行理论分析，用恰当的函数去模拟数据原型是一类十分重要的问题，最常用的逼近原则是让实测数据和估计数据之间的距离平方和最小，这即是最小二乘法。最小二乘法是一种经典的数据处理方法。在系统辨识领域中 ,最小二乘法是一种得到广泛应用的估计方法 ,可用于动态系统 ,静态系统 , 线性系统 ,非线性系统。可用于离线估计，也可用于在线估计。这种辨识方法主要用于在线辨识。在随机的环境下，利用最小二乘法时，并不要求观测数据提供其概率统计方面的信息，而其估计结果，却有相当好的统计特性。
MATLAB是一套高性能数字计算和可视化软件 ,它集成概念设计 ,算法开发 ,建模仿真 ,实时实现于一体 ,构成了一个使用方便、界面友好的用户环境 ,其强大的扩展功能为各领域的应用提供了基础。对于比较复杂的生产过程 ,由于过程的输入输出信号一般总是可以测量的 ,而且过程的动态特性必然表现在这些输入输出数据中 ,那么就可以利用输入输出数据所提供的信息来建立过程的数学模型。这种建模方法就称为系统辨识。把辨识建模称作“黑箱建模”。
2.2.2、 最小二乘法系统辨识结构：
本文把待辨识的过程看作“黑箱”。只考虑过程的输入输出特性，而不强调过程的内部机理。
+
+
e(k)
图1 SISO系统辨识“黑箱”结构图
y(k)
u(k)
z(k)
v(k)

图中，输入u(k)和输出z(k)是可以观测的；G（）是系统模型，用来描述系统的输入输出特性；N（）是噪声模型，v(k)是白噪声，e(k)是有色噪声，根据表示定理：
可以表示为
e(k) =N（）v(k)

2.2.3、 准则函数
设一个随机序列的均值是参数的线性函数： ，其中是可测的数据向量，那么利用随机序列的一个实现，使准则函数：
(式2-2)
达到极小的参数估计值称作的最小二乘估计。
最小二乘格式：，为模型参数向量，为零均值随机噪声。
2.3 最小二乘一次完成算法
2.3.1、 最小二乘问题的解
考虑系统模型:
(式3-1)
准则函数可写成：
（式3-2）
极小化准则函数得到：
（式3-3）
通过极小化式3-2计算的方法称作加权最小二乘法，为加权最小二乘估计，若取，则退化为一般最小二乘估计值，对应方法叫最小二乘法：
= （式3-4）
当获得一批数据后，利用式3-3或3-4可一次求得相应参数估计值，对应的称为最小二乘估计值，这样处理问题的方法就称作一次完成算法。
输入信号:
(1) 随机序列（如白噪声）；
(2) 伪随机序列（如M序列或逆M序列）；
(3) 离散序列
2.3.2、 最小二乘一次完成算法的MATLAB仿真（程序源代码见附录）
考虑仿真对象
z(k)= -1.642z(k-1)-0.715z(k-2)+0.39u(k-1)+0.35u(k-2)+v(k)
式中，v(k) 是均值为0，方差为0.01、0.1和0.5的不相关随机序列。输入信号采用4阶M序列，幅度为1。选择如下形式的辨识模型
，
构造ZL和HL，数据长度去L=30；加权阵取I，利用式（式3-4）计算参数估计值
, ，
程序框图如图3.2所示。Matlab6.0仿真程序如下：

赋输入信号初值u
定义输出观测值的长度并计算系统的输出值
显示输入和输出观测值
给样本矩阵HL和zL赋值
根据式(3-4)计算参数
从中分离出并显示出被辨识参数a1, a2, b1, b2
停机
图3.2 最小二乘一次完成算法程序框图

2.4 广义最小二乘法
2.4.1、 广义最小二乘数学模型

式中，u(k)和表示系统的输入输出；v(k)是均值为零的不相关的随机序列；且

2.4.2、 广义最小二乘递推算法如下

式中

2.4.3、 广义最小二乘递推算法的计算步骤：
1.给定初始条件
2利用式，计算和；
3利用式，构造；
4利用式递推计算；
5利用和
计算；
6根据来构造；
7利用
返回第2步进行迭代计算，直至获得满意的辨识结果。

2.4.4、 广义最小二乘递推算法的MATLAB仿真（程序源代码见附录）
考虑仿真对象
z(k)= -1.642z(k-1)-0.715z(k-2)+0.39u(k-1)+0.35u(k-2)+v(k)
式中，v(k) 是均值为0，方差为0.01、0.1和0.5的不相关随机序列。输入信号采用4阶M序列，幅度为1。

选择如下形式的辨识模型

+
y(k)
u(k)
e(k)
z(k)
+
v(k)

图2广义最小二乘法辨识实例结构图

其中取c1=0,c2=0.

3. 结果分析及算法优化
由于辨识算法中输入或噪声信号为不相关随机序列，所以每次辨识结果都不完全相同。但是，在相同输入、相同的噪声、相同的步长条件下，精度大体相同。
算法优化方案：（1）使用M序列（具有近似白噪声的性质）为输入信号；
（2）增加数据长度去L；
（3）减小噪声信号v(k)的方差。
3.1、 最小二乘一次完成算法仿真结果及分析：
采用给定的30组随机数据作为输入, 即数据长度去L=30，噪声v(k)选用均值为零，方差分别为0.5、0.1和0.01辨识结果如表3-1-1：
表3-1-1

真值
噪声方差为0.5
噪声方差为0.1
噪声方差为0.01
估计值
相对误差
估计值
相对误差
估计值
相对误差
a1
1.642
1.6184
-0.0144
1.6428
0.0005
1.6413
0.0004
a2
0.715
0.6851
-0.0418
0.7254
0.0145
0.7396
0.0344
b1
0.39
0.4535
0.1628
0.4149
0.0638
0.3941
0.0105
b2
0.35
0.3727
0.0649
0.3778
0.0794
0.3725
0.0643

由辨识结果可知，v(k)的方差越小，辨识效果越好。但由于v(k) 噪声是均值为零，方差分别为0.5、0.1和0.01的不相关随机序列，故一次完成算法每次辨识结果都不一样，得到的参数与真值间的误差较大,不太理想。
如V(k)选用均值为零，方差为0.1的不相关随机序列。增加数据长度L=300,效果如表3-1-2：
表3-1-2

真值
噪声方差为0.1，且输入为M序列，采集数据30个
噪声方差为0.1，且输入为M序列，采集数据300个
估计值
相对误差
估计值
相对误差
a1
1.642
1.6428
0.0005
1.6439
0.0011
a2
0.715
0.7254
0.0145
0.7125
-0.0034
b1
0.39
0.4149
0.0638
0.3963
0.0162
b2
0.35
0.3778
0.0794
0.3548
0.0137
通过数据比较（见上表3-1-1和3-1-2），可知增加数据长度L效果明显，但一次完成算法计算数据量很大；需要计算矩阵的逆，当数据很多时，会出还现“数据饱和”现象。
3.2、 广义最小二乘递推算法的的MATLAB仿真结果及分析
(1)、输入选用题目给出的30个随机数，即数据长度去L=30，噪声选用均值为零，方差分别为0.5、0.1和0.01的随机序列，辨识结果如表3-2-1
表中给出了三种情况下辨识参数结果即表中的估计值，估计值与真值的相对误差
表3-2-1

真值
噪声方差为0.5
噪声方差为0.1
噪声方差为0.01
估计值
相对误差
估计值
相对误差
估计值
相对误差
a1
1.642
1.4555
-0.1135
1.5630
0.1231
1.6523
0.0063
a2
0.715
0.5884
-0.1771
0.6392
0.1061
0.7315
0.0231
b1
0.39
0.4916
0.2605
0.3190
0.1821
0.4099
0.0510
b2
0.35
0.4213
0.2014
0.3461
0.0114
0.3418
0.0234

输入M序列，30步，噪声方差0.5时：
=12.556;
输入M序列，30步，噪声方差0.1时：
=2.5822;
输入M序列，30步，噪声方差0.01时：
=0.1706;

（2）、输入均采用M序列，噪声选择均值为零，方差为0.5、0.1和0.01的随机序列，辨识步长均为300步，辨识结果如表3-2-2。
表中给出了三种情况下辨识参数结果即表中的估计值，估计值与真值的相对误差.
表3-2-2

真值
噪声方差为0.5
噪声方差为0.1
噪声方差为0.01
估计值
相对误差
估计值
相对误差
估计值
相对误差
a1
1.642
1.5960
-0.0280
1.6550
0.0079
1.6429
0.0005
a2
0.715
0.6649
-0.0701
0.7359
0.0292
0.7101
-0.0068
b1
0.39
0.3413
-0.1249
0.4127
0.0582
0.3920
0.0051
b2
0.35
0.3212
-0.0823
0.3844
0.0983
0.3483
-0.0049

输入采用M序列，噪声方差0.5时： =135.288

输入采用M序列，噪声方差0.1时： =28.917

输入采用M序列，噪声方差0.01时：=2.7374

（3）数据结果分析：输入采用M序列比采用随机序列得到的辨识效果更好。噪声均值相等时，方差越大，辨识效果越差，反之，方差越小辨识效果越好。可以通过增加步长的方法提高辨识精度。
（4）下面给出递推辨识过程中参数及绝对误差变化过程，通过使用MATLAB作图，可以很直观地观测到变化过程。

下面给出以M序列作为输入，噪声均值为零，方差为0.01的随机序列，数据长度去L=30，得到的变化曲线图：

下面给出以M序列作为输入，噪声均值为零，方差为0.01的不相关随机序列，数据长度去L=300，得到的变化曲线图：

4. 参考文献
1) 侯媛彬, 汪梅, 王立琦,系统辨识及其MATLAB仿真,科学出版社,2004年
2) 方崇智，萧德云，过程辨识，清华大学出版社，1988年
3) 贾秋玲，袁冬莉，栾云凤，基于MATLAB7.x/Simulink/Stateflow系统仿真、分析及设计，西北工业大学出版社，2006年
4) 李言俊，张科，系统辨识理论及应用，国防工业出版社，2006年

12

附录
最小二乘一次完成算法的MATLAB仿真程序源代码：

clear %清理工作间变量
n=normrnd(0, sqrt(0.01), 1, 30)
u=[ 1.147 0.201 -0.787 -1.589 -1.502 0.866 1.152 1.573 0.626 0.433 -0.985 0.810 -0.044 0.974 -1.474 -0.719 -0.086 1.099 1.450 1.151 0.485 1.633 0.043 1.326 1.706 -0.340 0.890 0.144 1.177 -0.390];
z(2)=0;z(1)=0; %设z的前二个初始值为零
for k=3:31
z(k)=-1.642*z(k-1)-0.715*z(k-2)+0.39*u(k-1)+0.35*u(k-2)+n(k-1); %输出采样信号（测量值）
end
u,z %显示输入信号和输出观测信号
% 给样本矩阵ZL赋值
ZL=[z(3);z(4);z(5);z(6);z(7);z(8);z(9);z(10);z(11);z(12);z(13);z(14);z(15);z(16);z(17);z(18);z(19);z(20);z(21);z(22);z(23);z(24);z(25);z(26);z(27);z(28);z(29);z(30);z(31)];
%给样本矩阵HL赋值
h0=[-z(2) -z(1) u(2) u(1)]';
HLT=[h0,zeros(4,28)];
for k=3:30
h1=[-z(k) -z(k-1) u(k) u(k-1)]';
HLT(:,k-1)=h1;
end
HL=HLT'
c1=HL'*HL; c2=inv(c1); c3=HL'*ZL; c=c2*c3 %计算并显示
a1=c(1), a2=c(2), b1=c(3),b2=c(4) %从 中分离出并显示a1 、a2、 b1、 b2

广义最小二乘递推算法的MATLAB仿真程序源代码：

clear %清理工作间变量
L=300; % M序列的周期
y1=1;y2=1;y3=1;y4=0; %四个移位寄存器的输出初始值
for i=1:L;%开始循环，长度为L
x1=xor(y3,y4); %第一个移位寄存器的输入是第三个与第四个移位寄存器的输出的“或”
x2=y1; %第二个移位寄存器的输入是第一个移位寄存器的输出
x3=y2; %第三个移位寄存器的输入是第二个移位寄存器的输出
x4=y3; %第四个移位寄存器的输入是第三个移位寄存器的输出
y(i)=y4; %取出第四个移位寄存器的幅值为"0"和"1"的输出信号，即M序列
if y(i)>0.5,u(i)=-1; %如果M序列的值为"1", 辨识的输入信号取“-1”
else u(i)=1; %如果M序列的值为"0", 辨识的输入信号取“1”
end %小循环结束
y1=x1;y2=x2;y3=x3;y4=x4; %为下一次的输入信号做准备
end %大循环结束，产生输入信号u
figure(1); %第一个图形
stem(u),grid on %显示出输入信号M序列径线图并给图形加上网格
v=normrnd(0, sqrt(0.01), 1, 300);%均值为零的，方差为0.01或0.5或0.1不相关的随机噪声
ze(2)=0;ze(1)=0;
for k=3:301;
ze(k)=0*ze(k-1)+0*ze(k-2)+v(k-1);%C(z~1)=1,即取c1=0,c2=0
end
z(2)=0;z(1)=0; %设z的前两个初始值为零
for k=3:301; %循环变量从3到301
z(k)=-1.642*z(k-1)-0.715*z(k-2)+0.39*u(k-1)+0.351*u(k-2)+ze(k-1); %输出采样信号（测量值）
end
%RGLS广义最小二乘辨识
c0=[0.0001 0.0001 0.0001 0.0001]'; %直接给出被辨识参数的初始值,即一个充分小的实向量
pf0=10^6*eye(4,4); %直接给出初始状态P0，即一个充分大的实数单位矩阵
ce0=[0.001 0.001]';
pe0=eye(2,2);
c=[c0,zeros(4,299)]; %被辨识参数矩阵的初始值及大小
ce=[ce0,zeros(2,299)];
e=zeros(4,300); %相对误差的初始值及大小
ee=zeros(2,300);
s=0;
%广义最小二乘递推算法的计算步骤
for k=3:300;
zf(k)=z(k)+ce(1,k-2)*z(k-1)+ce(2,k-2)*z(k-2);
uf(k)=u(k)+ce(1,k-2)*u(k-1)+ce(2,k-2)*u(k-2);
hf1=[-zf(k-1),-zf(k-2),uf(k-1),uf(k-2)]';
x=hf1'*pf0*hf1+1; x1=inv(x); %开始求K(k)
k1=pf0*hf1*x1;%求出K的值
d1=zf(k)-hf1'*c0; c1=c0+k1*d1; %求被辨识参数c
e1=c1-c0; %求参数当前值与上一次的值的差值
e2=e1./c0; %求参数的相对变化
e(:,k)=e2; %把当前相对变化的列向量加入误差矩阵的最后一列
c0=c1; %新获得的参数作为下一次递推的旧参数
c(:,k)=c1; %把辨识参数c 列向量加入辨识参数矩阵的最后一列
pf1=pf0-k1*hf1'*pf0; %求出 p(k)的值
pf0=pf1; %给下次用
h1=[-z(k-1),-z(k-2),u(k-1),u(k-2)]';
s=s+(z(k)-h1'*[1.642 0.715 0.39 0.35]')^2;%求准则函数
ee(k)=z(k)-h1'*c1;
he1=[-ee(k-1),-ee(k-2)]';
x=he1'*pe0*he1+1; x1=inv(x);
k1=pe0*he1*x1;
d1=ee(k)-he1'*ce0;
ce1=ce0+k1*d1;
pe1=pe0-k1*he1'*pe0;
ce0=ce1;
ce(:,k)=ce1;
pe0=pe1;
end %大循环结束
c%辨识参数变化矩阵
%显示被辨识参数及其误差(收敛)情况
%分离参数
a1=c(1,1:300); a2=c(2,1:300);b1=c(3,1:300);b2=c(4,1:300); c1=ce(1,1:300);c2=ce(2,1:300);
ea1=e(1,1:300); ea2=e(2,1:300);eb1=e(3,1:300); eb2=e(4,1:300);
figure(2); %第二个图形
i=1:300; %横坐标从1到300
plot(i,a1,'r',i,a2,'k',i,b1,'b',i,b2,'c',i,c1,'b',i,c2,'r') %画出a1，a2,b1，b2,c1,c2的各次辨识结果
title('参数变化曲线') %图形标题
figure(3); %第三个图形
i=1:300; %横坐标从1到300
plot(i,ea1,'r',i,ea2,'k:',i,eb1,'b',i,eb2,'k:') %画出a1，a2，b1，b2,的各次辨识结果的收敛情况
title('误差曲线') %图形标题

课程设计成绩评定表
质量评价指标（在相应栏目打√）
评 价 项 目
评 价 质 量
优秀
良好
一般
及格
不及格
工作量和态度

实验、计算可靠性

文字和图表质量

总体评价

评定成绩（百分制）

评定小组成员签名

2010年 月 日

备 注

16