摘要系统辨识是描述各种各样系统运动规律的一种方法论，是研究系统的一种有效工具。

利用这个工具可以对我们要研究的系统进行定量描述。

随着现代控制理论的迅速发展，系统辨识得到迅速而蓬勃发展，并已经成功运用与多种工程应用领域。

但针对有色噪声干扰系统，传统的辨识方法不能得到良好的参数估计，而工程上大多系统都为有色噪声干扰系统。

有色噪声干扰系统的一类系统为广义输出误差模型，本文对白噪声和有色噪声两种噪声干扰下运用最小二乘法，递推最小二乘法分别比较噪声的不同，以及进一步得出最小二乘法的适用性。

进一步研究广义最小二乘法对有色噪声系统辨识的改进。

关键词：系统辨识；白噪声；有色噪声；最小二乘；递推最小二乘；广义最小二乘ABSTRACTthe law of motion of the system identification is to describe the various system and an effective tool for the system. It also can be used to study quantitative description of the system. With the rapid development of modern control theory,the development of system identification has been rapid and vigorous , and has been successfully applied to a variety of engineering applications. However, the traditional identification methods can't get good parameter estimation in engineering systems for colored noise . Colored noise jamming system is for generalized output error model.under two kinds of noise of white noise and color noise ,by using the least square method and the recursive least square method ,respectively to compare the different noise, and to furtherconcluded that the applicability of the least-square method. To further studied the improvement of generalized least squares identification of colored noise system.KEY WORDS: system identification,white noise, colored noise.,least squares;,recursive least squares,the generalized least squares目录摘要 (1)第一章系统辨识 (2)1.1辨识的定义 (2)1.2辨识的内容和步骤 (3)1.3辨识目的 (3)1.4先验知识 (3)1.5数据预处理 (4)1.6模型结构辨识 (4)1.7模型参数辨识 (4)1.8模型检验 (4)第二章题目详解 (5)2.1.1最小二乘法 (5)2.1.2递推法 (8)2.2.1最小二乘法 (12)2.2.2递推法 (13)2.3广义最小二乘法 (15)附录 (16)第一章系统辨识1.1辨识的定义Zadeh对辨识的定义（1962年）辨识就是在输入和输出数据的基础上，从一组给定的模型类中，确定一个与所测系统等价的模型。

Liung 对辨识的的定义（1978年）系统辩识有三个要素——数据、模型类和（等价）准则。

系统辩识是按照一个准则，在模型类中选择一个与数据拟合得最好的模型。

Liung 认为，实际系统的复杂性很难找到一个适用的模型与之等价。

因此，系统辩识的任务只是要求从输入输出数据出发，找到一个与实际系统相逼近的模型。

该定义体现了逼近的观点。

系统辨识是在输入和输出观测的基础上，从指定的一类模型中确定一个与系统等价的模型。

1.2辨识的内容和步骤简单地说，辨识就是一种从观测到的含有噪声的输入输出数据中提取数学模型的方法。

根据现场的情况，辨识可以离线进行，也可以在线进行；辨识的内容主要包括四个方面：实验设计，模型结构辨识，模型参数辨识，模型检验。

1.3辨识目的明确模型应用的最终目的很重要，因为它将决定模型的类型，精度要求和采用哪种辨识方法等问题。

1.4先验知识对于给定的过程进行辨识之前，要通过一些手段对过程有一些了解粗略的掌握一些过程和数据，这些先验知识对实验设计将其指导作用。

1.5数据预处理输入输出数据通常含有直流或交流成分，任何方法都无法消除对系统精度的影响。

此外，高频成分也是不利的。

因此需对输入数据进行零均值化和剔除高频成分预处理。

1.6模型结构辨识模型结构辨识包括模型验前结构和模型结构参数两部分。

1.7模型参数辨识当模型结构确定之后，就需要进行模型参数辨识。

方法很多，最基本，应用最广的是最小二乘法。

但是最小二乘法也有一些重大缺陷，比如过程辨识受到有色噪声污染时，他几乎比能适应。

1.8模型检验模型检验是辨识过程不可或缺的步骤之一。

但是，它没有一般方法可循。

它和模型结构问题密切相关。

如果模型结构不合适，模型检验是不能通过的。

第二章 题目详解2.1.1最小二乘法一个单输入－单输出线性定常系统的差分方程)()1()()()2()1()(1021n k u b L k u b k u b n k x a L k x a k x a k x n n -++-+=-++-+-+式中)(k u 为输入信号，)(k x 为理论上的输出值。

)(k x 通过观测得到，在观测过程中往往附加有随机干扰。

观测值用)(k y 表示：)()()(k v k x k y += ξ（k ）=v （k ）+Σ)(i k v a i - 则∑∑+-+--=)()()()(k i k u b i k y a k y i i ξ；如果 u （k ）也有测量误差，则在ξ（k ）中应包含这一测量误差。

ζφθ+=Y 则()YTTθθφφ1ˆ-=原题设为二阶系统，且θ的真值已知，输入)(k u 已知，求)(k y ，并且由输入输出序列得出参数估值，其程序已在matlab中实现。

实验结果如下图：2.1.2递推法图1.2 最小二乘递推算法辨识的Malab6.0程序流程图Matlab仿真结果：c =Columns 1 through 121.6420 0 1.6420 0 0 0 0 0 0 0 0 00.7150 0 0.7150 0 0 0 0 0 0 0 0 00.3900 0 0.3689 0 0 0 0 0 0 0 0 00.3500 0 0.2298 0 0 0 0 0 0 0 0 0Columns 13 through 240 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 00 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0Columns 25 through 300 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 0e =Columns 1 through 120 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 00 0 -0.0540 0 0 0 0 0 0 0 0 00 0 -0.3433 0 0 0 0 0 0 0 0 0Columns 13 through 240 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0Columns 25 through 300 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 0结果分析：仿真结果表明，大约递推到第五步时，参数辨识的结果基本达到稳定状态。

此时，参数的相对变化量E ≤0.000000005。

从整个辨识过程来看，精度的要求直接影响辨识的速度。

虽然最终的精度可以达到很小，但开始阶段的相对误差会很大。

2.2.1最小二乘法题设将白噪声改为有色噪声，)2()1()()(21-+-+=k a k a k k εεεξ；有色噪声可以看做是白噪声经过成型滤波器后得到的。

Matlab 仿真结果：结果分析：由于所用的输出观测值有有色噪声成分，所以辨识结果有误差，且用最小二乘法计算得到的参数估计误差较大。

2.2.2递推法结果：结果分析：第一问和第二问的区别只是噪声不同，但是从仿真结果图像来看，都是从第五步开始达到稳定状态的，但有色噪声参数辨识显然误差更大更明显。

2.3广义最小二乘法广义最小二乘法辨识的计算步骤如下：（1）应用已得到的输入和输出数据 u(k) 和y(k) (2)计算残差e(k),并用其代替)(k ξ （3）计算和)(k u -和)(k y -(4)应用得到的)(k u -和)(k y -，)()()()()()2(1)1(1k k b k a u z y z ε+=--，再用最小二乘法重新估计θ ，得θ的第2次估值 θ)2( 。

然后按步骤⑵计算残差 e(k) 。

重新估计f ，得到估值 f)2( 。

再按步骤⑶计算y)2(和u)2(，按步骤⑷求θ 的第3次估计θ)3( 。

重复上循环步骤，直到θ 的估值θ收敛为止。

附录程序%1.1 二阶系统的最小二乘一次完成算法辨识程序，在光盘中的文件名：FLch3LSeg1.m u=[1.147,0.201,-0.787,-1.589,-1.052,0.866,1.152,1.573,0.626,0.433,-0.958,0.81,-0.044,0.947,-1.4 74,-0.719,-0.086,-1.099,1.45,1.151,0.485,1.633,0.043,1.326,1.706,-0.340,0.89,1.444,1.177]; %系统辨识的输入信号为一个周期的M序列y=zeros(1,30); %定义输出观测值的长度w=0.1*randn(1,30)for k=3:30y(k)=-1.642*y(k-1)-0.715*y(k-2)+0.39*u(k-1)+0.35*u(k-2)+w(k); %用理想输出值作为观测值endsubplot(3,1,1) %画三行一列图形窗口中的第一个图形stem(u) %画输入信号u的径线图形subplot(3,1,2) %画三行一列图形窗口中的第二个图形i=1:1:30; %横坐标范围是1到16，步长为1plot(i,y) %图形的横坐标是采样时刻i, 纵坐标是输出观测值z, 图形格式为连续曲线subplot(3,1,3) %画三行一列图形窗口中的第三个图形stem(y),grid on %画出输出观测值z的径线图形，并显示坐标网格u,y %显示输入信号和输出观测信号%L=30%数据长度HL=[-y(2) -y(1) u(2) u(1);-y(3) -y(2) u(3) u(2);-y(4) -y(3) u(4) u(3);-y(5) -y(4) u(5) u(4);-y(6) -y(5) u(6) u(5);-y(7) -y(6) u(7) u(6);-y(8) -y(7) u(8) u(7);-y(9) -y(8) u(9) u(8);-y(10) -y(9) u(10) u(9);-y(11) -y(10) u(11) u(10);-y(12) -y(11) u(12) u(11);-y(13) -y(12) u(13) u(12);-y(14) -y(13) u(14) u(13);-y(15) -y(14) u(15) u(14);-y(16) -y(15) u(16) u(15);-y(17) -y(16) u(17) u(16);-y(18) -y(17) u(18) u(17);-y(19) -y(18) u(19) u(18);-y(20) -y(19) u(20) u(19);-y(21) -y(20) u(21) u(20);-y(22) -y(21) u(22) u(21);-y(23) -y(22) u(23) u(22);-y(24) -y(23) u(24) u(23);-y(25) -y(24) u(25) u(24);-y(26) -y(25) u(26) u(25);-y(27) -y(26) u(27) u(26);-y(28) -y(27) u(28) u(27);-y(29) -y(28) u(29) u(28)] %给样本矩阵HL赋值Y=[y(3);y(4);y(5);y(6);y(7);y(8);y(9);y(10);y(11);y(12);y(13);y(14);y(15);y(16);y(17);y(18);y(19) ;y(20);y(21);y(22);y(23);y(24);y(25);y(26);y(27);y(28);y(29);y(30)] % 给样本矩阵z L赋值%Calculating Parametersc1=HL'*HL; c2=inv(c1); c3=HL'*Y; c=c2*c3 %计算并显示a1=c(1), a2=c(2), b1=c(3),b2=c(4) %从中分离出并显示a1 、a2、b1、b2%1.2 RLS递推最小二乘辨识c0=[1.642 0.715 0.39 0.35]'; %被辨识参数的初始值p0=10^6*eye(4,4); %直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.000000005; %取相对误差E=0.000000005c=[c0,zeros(4,29)]; %被辨识参数矩阵的初始值及大小e=zeros(4,30); %相对误差的初始值及大小u=[1.147,0.201,-0.787,-1.589,-1.052,0.866,1.152,1.573,0.626,0.433,-0.958,0.81,-0.044,0.947,-1.4 74,-0.719,-0.086,-1.099,1.45,1.151,0.485,1.633,0.043,1.326,1.706,-0.340,0.89,1.444,1.177]; %系统辨识的输入信号为一个周期的M序列y=zeros(1,30); %定义输出观测值的长度w=0.5*randn(1,30);for k=3:30y(k)=-1.642*y(k-1)-0.715*y(k-2)+0.39*u(k-1)+0.35*u(k-2)+w(k); %用理想输出值作为观测值endfor k=3:30; %开始求Kh1=[-y(k-1),-y(k-2),u(k-1),u(k-2)]';x=h1'*p0*h1+1;x1=inv(x); %开始求K(k)k1=p0*h1*x1;%求出K的值d1=y(k)-h1'*c0;c1=c0+k1*d1; %求被辨识参数ce1=c1-c0; %求参数当前值与上一次的值的差值e2=e1./c0; %求参数的相对变化e(:,k)=e2; %把当前相对变化的列向量加入误差矩阵的最后一列c0=c1; %新获得的参数作为下一次递推的旧参数c(:,k)=c1; %把辨识参数c 列向量加入辨识参数矩阵的最后一列p1=p0-k1*k1'*(h1'*p0*h1+1); %求出p(k)的值p0=p1; %给下次用if e2<=Ebreak; %如果参数收敛情况满足要求，终止计算end %小循环结束end %大循环结束c,e %显示被辨识参数及其误差(收敛)情况%分离参数a1=c(1,:); a2=c(2,:); b1=c(3,:); b2=c(4,:); ea1=e(1,:); ea2=e(2,:); eb1=e(3,:); eb2=e(4,:);figure(2); %第二个图形i=1:30; %横坐标从1到30plot(i,a1,'r',i,a2,':',i,b1,'g',i,b2,':') %画出a1，a2，b1，b2的各次辨识结果title('Parameter Identification with Recursive Least Squares Method') %图形标题figure(3); %第三个图形i=1:30; %横坐标从1到30plot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:') %画出a1，a2，b1，b2的各次辨识结果的收敛情况title('Identification Precision') %图形标题%2.1 二阶系统的最小二乘一次完成算法辨识程序，在光盘中的文件名：FLch3LSeg1.m u=[1.147,0.201,-0.787,-1.589,-1.052,0.866,1.152,1.573,0.626,0.433,-0.958,0.81,-0.044,0.947,-1.4 74,-0.719,-0.086,-1.099,1.45,1.151,0.485,1.633,0.043,1.326,1.706,-0.340,0.89,1.444,1.177]; %系统辨识的输入信号为一个周期的M序列y=zeros(1,30); %定义输出观测值的长度w=0.1*randn(1,30);for k=3:30W(k)=w(k)+1.642*w(k-1)+0.715*w(k-2)y(k)=-1.642*y(k-1)-0.715*y(k-2)+0.39*u(k-1)+0.35*u(k-2)+W(k); %用理想输出值作为观测值endsubplot(3,1,1) %画三行一列图形窗口中的第一个图形stem(u) %画输入信号u的径线图形subplot(3,1,2) %画三行一列图形窗口中的第二个图形i=1:1:30; %横坐标范围是1到16，步长为1plot(i,y) %图形的横坐标是采样时刻i, 纵坐标是输出观测值z, 图形格式为连续曲线subplot(3,1,3) %画三行一列图形窗口中的第三个图形stem(y),grid on %画出输出观测值z的径线图形，并显示坐标网格u,y %显示输入信号和输出观测信号%L=30%数据长度HL=[-y(2) -y(1) u(2) u(1);-y(3) -y(2) u(3) u(2);-y(4) -y(3) u(4) u(3);-y(5) -y(4) u(5) u(4);-y(6) -y(5) u(6) u(5);-y(7) -y(6) u(7) u(6);-y(8) -y(7) u(8) u(7);-y(9) -y(8) u(9) u(8);-y(10) -y(9) u(10) u(9);-y(11) -y(10) u(11) u(10);-y(12) -y(11) u(12) u(11);-y(13) -y(12) u(13) u(12);-y(14) -y(13) u(14) u(13);-y(15) -y(14) u(15) u(14);-y(16) -y(15) u(16) u(15);-y(17) -y(16) u(17) u(16);-y(18) -y(17) u(18) u(17);-y(19) -y(18) u(19) u(18);-y(20) -y(19) u(20) u(19);-y(21) -y(20) u(21) u(20);-y(22) -y(21) u(22) u(21);-y(23) -y(22) u(23) u(22);-y(24) -y(23) u(24) u(23);-y(25) -y(24) u(25) u(24);-y(26) -y(25) u(26) u(25);-y(27) -y(26) u(27) u(26);-y(28) -y(27) u(28) u(27);-y(29) -y(28) u(29) u(28)] %给样本矩阵HL赋值Y=[y(3);y(4);y(5);y(6);y(7);y(8);y(9);y(10);y(11);y(12);y(13);y(14);y(15);y(16);y(17);y(18);y(19) ;y(20);y(21);y(22);y(23);y(24);y(25);y(26);y(27);y(28);y(29);y(30)] % 给样本矩阵z L赋值%Calculating Parametersc1=HL'*HL; c2=inv(c1); c3=HL'*Y; c=c2*c3 %计算并显示a1=c(1), a2=c(2), b1=c(3),b2=c(4) %从中分离出并显示a1 、a2、b1、b2%2.2 RLS递推最小二乘辨识c0=[1.642 0.715 0.39 0.35]'; %被辨识参数的初始值p0=10^6*eye(4,4); %直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.000000005; %取相对误差E=0.000000005c=[c0,zeros(4,29)]; %被辨识参数矩阵的初始值及大小e=zeros(4,30); %相对误差的初始值及大小u=[1.147,0.201,-0.787,-1.589,-1.052,0.866,1.152,1.573,0.626,0.433,-0.958,0.81,-0.044,0.947,-1.4 74,-0.719,-0.086,-1.099,1.45,1.151,0.485,1.633,0.043,1.326,1.706,-0.340,0.89,1.444,1.177]; %系统辨识的输入信号为一个周期的M序列y=zeros(1,30); %定义输出观测值的长度w=0.5*randn(1,30);for k=3:30W(k)=w(k)+1.642*w(k-1)+0.715*w(k-2)y(k)=-1.642*y(k-1)-0.715*y(k-2)+0.39*u(k-1)+0.35*u(k-2)+W(k); %用理想输出值作为观测值endfor k=3:30; %开始求Kh1=[-y(k-1),-y(k-2),u(k-1),u(k-2)]';x=h1'*p0*h1+1;x1=inv(x); %开始求K(k)k1=p0*h1*x1;%求出K的值d1=y(k)-h1'*c0;c1=c0+k1*d1; %求被辨识参数ce1=c1-c0; %求参数当前值与上一次的值的差值e2=e1./c0; %求参数的相对变化e(:,k)=e2; %把当前相对变化的列向量加入误差矩阵的最后一列c0=c1; %新获得的参数作为下一次递推的旧参数c(:,k)=c1; %把辨识参数c 列向量加入辨识参数矩阵的最后一列p1=p0-k1*k1'*(h1'*p0*h1+1); %求出p(k)的值p0=p1; %给下次用if e2<=Ebreak; %如果参数收敛情况满足要求，终止计算end %小循环结束end %大循环结束c,e %显示被辨识参数及其误差(收敛)情况%分离参数a1=c(1,:); a2=c(2,:); b1=c(3,:); b2=c(4,:); ea1=e(1,:); ea2=e(2,:); eb1=e(3,:); eb2=e(4,:);figure(2); %第二个图形i=1:30; %横坐标从1到30plot(i,a1,'r',i,a2,':',i,b1,'g',i,b2,':') %画出a1，a2，b1，b2的各次辨识结果title('Parameter Identification with Recursive Least Squares Method') %图形标题figure(3); %第三个图形i=1:30; %横坐标从1到30plot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:') %画出a1，a2，b1，b2的各次辨识结果的收敛情况title('Identification Precision') %图形标题%3 广义最小二乘递推法u=[1.147,0.201,-0.787,-1.589,-1.052,0.866,1.152,1.573,0.626,0.433,-0.958,0.81,-0.044,0.947,-1.4 74,-0.719,-0.086,-1.099,1.45,1.151,0.485,1.633,0.043,1.326,1.706,-0.340,0.89,1.444,1.177]; %系统辨识的输入信号为一个周期的M序列y=zeros(1,30); %定义输出观测值的长度w=0.1*randn(1,30);for k=3:30y(k)=-1.642*y(k-1)-0.715*y(k-2)+0.39*u(k-1)+0.35*u(k-2)+w(k); %用理想输出值作为观测值endcf0=[0.001 0.001 0.001 0.001]';pf0=10^6*eye(4);E=0.000000005;ef=zeros(4,30);nc=2; % 设定噪声模型的最高阶次ce0=zeros(nc,1);pe0=eye(nc);v1=0;for k=nc+1:30zf(k)=[1 ce0']*y(k:-1:k-nc)';uf(k)=[1 ce0']*u(k:-1:k-nc)'; % 计算zf(k)、uf(k)hf=[-zf(k-1),-zf(k-2),uf(k-1),uf(k-2)]'; % 构造hf(k)jf=hf'*pf0*hf+1; jf=inv(jf);kf=pf0*hf*jf; % 准备并求出kf的值df=zf(k)-hf'*cf0;cf1=cf0+kf*df; % 求出theta的估计值ef1=(cf1-cf0)./cf0; % 求出两次估计值之间的相对误差并保存pf1=(eye(4)-kf*hf')*pf0; % 求出pf的值% 以下为噪声模型参数估计h=[-y(k-1),-y(k-2),u(k-1),u(k-2)]';v1(k)=y(k)-h'*cf1;he=-v1(k-1:-1:k-nc)'; % 构造噪声模型参数估计时的he(k)ze=he'*pe0*he+1; ze=inv(ze);ke=pe0*he*ze; % 准备并求出ke的值de=v1(k)-he'*ce0; ce1=ce0+ke*de; % 求出theta的估计值pe1=(eye(nc)-ke*he')*pe0; % 求出pe的值ee1=(ce1-ce0)./(ce0+eps);% 保存各种参数供下次循环和作图用cf0=cf1;cf(:,k)=cf1;pf0=pf1;ef0=ef1;ef(:,k)=ef1;ce0=ce1;pe0=pe1;ee0=ee1;% 如果满足以下条件，则跳出循环if (abs(ef1)<E)&&(ef1~=0)=0if abs(ee1)<Ej1=k;break;endendend% 分离估计参数a1=cf(1,:);a2=cf(2,:);b1=cf(3,:);b2=cf(4,:);ea1=ef(1,:);ea2=ef(2,:);eb1=ef(3,:);eb2=ef(4,:);end% 绘图figurej=1:k;plot(j,a1,'r',j,a2,'b',j,b1,'m',j,b2,'k');title('NSR=0.1时有色噪声影响下广义最小二乘递推法系统参数辨识结果'); xlabel('噪声模型的最高阶次');ylabel('系统参数辨识结果');axis tightlegend('a1=-1.5','a2=0.7','b1=2','b2=0.5')grid on;figureplot(j,ea1,'r',j,ea2,'b',j,eb1,'m',j,eb2,'k');title('NSR=0.1时有色噪声影响下广义最小二乘递推法系统参数辨识误差'); xlabel('递推次数');ylabel('系统参数辨识误差');axis([1,1024,-0.02,0.02]);grid on;legend('a1','a2','b1','b2')- 21 -。