最小二乘法参数辨识
201403027
摘 要:系统辨识在工程中的应用非常广泛, 系统辨识的方法有很多种, 最小
二乘法是一种应用极其广泛的系统辨识方法. 阐述了动态系统模型的建立及其最小二乘法在系统辨识中的应用, 并通过实例分析说明了最小二乘法应用于系统辨识中的重要意义.
关键词:最小二乘法; 系统辨识; 动态系统
Abstract : System identification in engineering is widely used, system identification methods there are many ways， least squares method is a very wide range of application of system identification method and the least squares method elaborated establish a dynamic system models in System Identification applications and examples analyzed by the least squares method is applied to illustrate the importance of system identification.
Keywords : Least Squares; system identification; dynamic system
引言
随着科学技术的不断发展, 人们认识自然、利用自然的能力越来越强, 对于未知对象的探索也越来越深入. 我们所研究的对象, 可以依据对其了解的程度分为三种类型:白箱、灰箱和黑箱. 如果我们对于研究对象的内部结构、内部机制了解很深入的话, 这样的研究对象通常称之为“白箱”; 而有的研究对象, 我们对于其内部结构、机制只了解一部分, 对于其内部运行规律并不十分清楚, 这样的研究对象通常称之为“灰箱”; 如果我们对于研究对象的内部结构、内部机制及运行规律均一无所知的话, 则把这样的研究对象称之为“黑箱”. 研究灰箱和黑箱时, 将研究的对象看作是一个系统, 通过建立该系统的模型, 对模型参数进行辨识来确定该系统的运行规律. 对于动态系统辨识的方法有很多, 但其中应用最广泛, 辨识 效果良好的就是最小二乘辨识方法, 研究最小二乘法在系统辨识中的应用具有现实的、广泛的意义.
1．1 系统辨识简介
系统辨识是根据系统的输入输出时间函数来确定描述系统行为的数学模型。现代控制理论中的一个分支。通过辨识建立数学模型的目的是估计表征系统行为的重要参数，建立一个能模仿真实系统行为的模型，用当前可测量的系统的输入和输出预测系统输出的未来演变，以及设计控制器。对系统进行分析的主要问题是根据输入时间函数和系统的特性来确定输出信号。对系统进行控制的主要问题是根据系统的特性设计控制输入，使输出满足预先规定的要求。而系统辨识所研究的问题恰好是这些问题的逆问题。通常，预先给定一个模型类μ={M}（即给定一类已知结构的模型），一类输入信号u 和等价准则J=L(y，yM)(一般情况下，J 是误差函数，是过程输出y 和模型输出yM 的一个泛函) ；然后选择使误差函数J 达到最小的模型，作为辨识所要求的结果。系统辨识包括两个方面：结构辨识和参数估计。在实际的辨识过程中，随着使用的方法不同，结构辨识和参数估计这两个方面并不是截然分开的，而是可以交织在一起进行的。
1.2系统辨识的目的
在提出和解决一个辨识问题时，明确最终使用模型的目的是至关重要的。它对模型类（模型结构）、输入信号和等价准则的选择都有很大的影响。通过辨识建立数学模型通常有四个目的。
①估计具有特定物理意义的参数 有些表征系统行为的重要参数是难以直接测量的，例如在生理、生态、环境、经济等系统中就常有这种情况。这就需要通过能观测到的输入输出数据，用辨识的方法去估计那些参数。
②仿真 仿真的核心是要建立一个能模仿真实系统行为的模型。用于系统分析的仿真模型要求能真实反映系统的特性。用于系统设计的仿真，则强调设计参数能正确地符合它本身的物理意义。
③预测 这是辨识的一个重要应用方面，其目的是用迄今为止系统的可测量的输入和输出去预测系统输出的未来的演变。例如最常见的气象预报，洪水预报，其他如太阳黑子预报，市场价格的预测，河流污染物含量的预测等。预测模型辨识的等价准则主要是使预测误差平方和最小。只要预测误差小就是好的预测
模型，对模型的结构及参数则很少再有其他要求。这时辨识的准则和模型应用的目的是一致的，因此可以得到较好的预测模型。
④控制 为了设计控制系统就需要知道描述系统动态特性的数学模型，建立这些模型的目的在于设计控制器。建立什么样的模型合适，取决于设计的方法和准备采用的控制策略。
2最小二乘方法
2.1.1 系统辨识最小二乘方法简介
对工程实践中测得的数据进行理论分析，用恰当的函数去模拟数据原型是一类十分重要的问题，最常用的逼近原则是让实测数据和估计数据之间的距离平方和最小，这即是最小二乘法。最小二乘法是一种经典的数据处理方法。在系统辨识领域中 ,最小二乘法是一种得到广泛应用的估计方法 ,可用于动态系统 ,静态系统 , 线性系统 ,非线性系统。可用于离线估计，也可用于在线估计。这种辨识方法主要用于在线辨识。在随机的环境下，利用最小二乘法时，并不要求观测数据提供其概率统计方面的信息，而其估计结果，却有相当好的统计特性。
MATLAB 是一套高性能数字计算和可视化软件 , 它集成概念设计 , 算法开发 , 建模仿真 , 实时实现于一体 , 构成了一个使用方便、界面友好的用户环境 , 其强大的扩展功能为各领域的应用提供了基础。对于比较复杂的生产过程 , 由于过程的输入输出信号一般总是可以测量的 , 而且过程的动态特性必然表现在这些输入输出数据中 , 那么就可以利用输入输出数据所提供的信息来建立过程的数学模型。这种建模方法就称为系统辨识。把辨识建模称作“黑箱建模”。
2.1.2 最小二乘法系统辨识结构：
本文把待辨识的过程看作“黑箱”。只考虑过程的输入输出特性，而不强调过程的内部机理。
图中，输入u(k)和输出z(k)是可以观测的；G 是系统模型，用来描述系统的输入输出特性；N 是噪声模型，v(k)是白噪声，e(k)是有色噪声，根据表示定理： 可以表示为
e(k) =N v(k)
G (z
-1
)
B (z -1) =
A (z -1)
N (z
-1
)
D (z -1) = -1
C (z )
-n a -1-1-2
⎧A (z ) =1+a z +a z ++a z 12n a ⎪⎨-n b -1-1-2B (z ) =b z +b z ++b z ⎪12n b ⎩
-n a -1-1-2⎧C (z ) =1+c z +c z + +c z ⎪12n a ⎨-n b -1-1-2D (z ) =d z +d z + +d z ⎪12n b ⎩
2.1.3 准则函数
设一个随机序列{z (k ), k ∈(1, 2, , L ) }的均值是参数θ的线性函数： E {z (k ) }=h T (k ) θ，其中h (k ) 是可测的数据向量，那么利用随机序列的一个实现，使准则函数：
J (θ) =∑[z (k ) -h T (k ) θ]
k =1
∧
L
2
(式2-2)
达到极小的参数估计值θ称作θ的最小二乘估计。
t
z (k ) =h (k ) θ+e (k ) ，θ为模型参数向量，e (k )为零均值随机 最小二乘格式：
2.2 广义最小二乘法
2.2.1 广义最小二乘数学模型
A (z -1) z (k ) =B (z -1) u (k ) +
1
v (k ) -1
C (z )
式中，u(k)和z (k ) 表示系统的输入输出；v(k)是均值为零的不相关的随机序列；且
⎧A (z -1) =1+a 1z -1+a 2z -2+ +a n a z -n a ⎪-n b -1-1-2B (z ) =b z +b z + +b z 12n ⎨b ⎪-n c -1-1-2C (z ) =1+c z +c z + +c z 12n c ⎩
2.2.2 广义最小二乘递推算法如下
ˆ(k ) =θˆ(k -1) +K (k )[z (k ) -h (k ) θˆ(k -1)]⎧θ⎪
⎪K (k ) =P (k -1) h (k ) [h (k ) P (k -1) h (k ) +1]⎪⎪P (k ) =[I -K (k ) h (k )]P (k -1) ⎨ˆˆˆ(k -1)]ˆ(k ) -h (k ) θ⎪θ(k ) =θ(k -1) +K (k )[e
⎪K (k ) =P (k -1) h (k ) [h (k ) P (k -1) h (k ) +1]⎪⎪⎩P (k ) =[I -K (k ) h (k )]P (k -1)
τ
f
f
f
τ
-1
f f f f f f
τ
f f f f
τ
e e e e e
τ
-1
e e e e e e
τ
e e e e
式中
⎧h (k ) =[-z (k -1), , -z (k -n ), u (k -1), , u (k -n )]⎪h (k ) =[-e ˆ(k -1), , -e ˆ(k -n )]⎨⎪e
⎩ˆ(k ) =z (k ) -h (k ) θˆ(k )
f
f
f
a τ
f
f
b
e
c
τ
τ
2.2.3 广义最小二乘递推算法的计算步骤：
ˆ(0) =ε(充分小的实向量) ⎧θ⎪
⎪P (0) =a I (a 为充分大的数)
1. 给定初始条件 ⎨
ˆ(0) =0⎪θ
⎪P (0) =I ⎩
z (k ) =C (z) z (k )
2利用式，计算z (k ) 和u (k ) ；
u (k ) =C (z) z (k )
2
f e e -1
f
-1
f f
f
⎧θ=[a , , a , b , , b ]
3利用式⎨，构造h (k ) ；
⎩h (k ) =[-z (k -1), , -z (k -n ), u (k -1), , u (k -n )]
τ
1
n a
1
n b
τf
f f f a f f b
ˆ(k ) =θˆ(k -1) +K (k )[z (k ) -h (k ) θˆ(k -1)]⎧θ⎪⎪ˆ(k ) ； 4利用式⎨K (k ) =P (k -1) h (k ) [h (k ) P (k -1) h (k ) +1]递推计算θ⎪P (k ) =[I -K (k ) h (k )]P (k -1) ⎪⎩
τ
f
f
f
τ
-1
f f f f f f
τ
f f f f
ˆ(k ) 和 ˆ(k ) =z (k ) -h (k ) θ5利用e
τ
ˆ(k ) ； h (k ) =[-z (k -1), , -z (k -n ), u (k -1), , u (k -n )]计算e
τ
a
b
ˆ(k -1), , -e ˆ(k -n c )]来构造h e (k ) ； 6根据h e (k ) =[-e
τ
ˆ(k ) =θˆ(k -1) +K (k )[e ˆ(k -1)]⎧θˆ(k ) -h (k ) θ⎪
7利用⎨K (k ) =P (k -1) h (k ) [h (k ) P (k -1) h (k ) +1]
⎪P (k ) =[I -K (k ) h (k )]P (k -1) ⎩
τ
e
e
e
e
e
τ
-1
e e e e e e
τ
e e e e
返回第2步进行迭代计算，直至获得满意的辨识结果。
3 工程实例 3.1典型系统建模
以某微循环流体系统模型的参数辨识为例. 我们已经得到该系统模型的差分方程形式, 取特定点的压力波作为模型的输入, 以另一点的压力波作为模型的输出. 由于我们采集的数据是实时的, 因此采用在线辨识方法. 由于建立的微循环流体系统模型是一个单输入、单输出的模型, 为使参数估计的结果很好地跟踪参数真值的变化, 我们采用渐消记忆的最小二乘法对系统模型参数进行辨识, 即强调新数据的作用, 贬低老数据的作用. 抽象出的SISO 系统的差分方程为：
（1-1） z (k ) +a 1z (k -1) +a 2z (k -2) =b 1u (k -1) +b 2u (k -2) +υ(k ) 式
T
θ 0.483 0.57 0.42]，利用MATLAB 的M 语言辨识参数取真值为：=[1.376
系统中的未知参数a 1、a 2、b 1、b 2。要求：用参数的真值利用差分方程求出z (k )
作为测量值，υ(k ) 是均值为0，方差为0.1、0.5和0.01的不相关随机序列。使用最小二乘算法辨识。
3.2 广义最小二乘递推算法的MATLAB 仿真（程序源代码见附录） 考虑仿真对象
z(k)= -1.376z(k-1)-0.483z(k-2)+0.57u(k-1)+0.42u(k-2)+v(k)
式中，v(k) 是均值为0，方差为0.01、0.1和0.5的不相关随机序列。输入信号采用4阶M 序列，幅度为1。 选择如下形式的辨识模型
其中取c1=0,c2=0.
4结果分析及算法优化
由于辨识算法中输入或噪声信号为不相关随机序列，所以每次辨识结果都不完全相同。但是，在相同输入、相同的噪声、相同的步长条件下，精度大体相同。
算法优化方案：（1）使用M 序列（具有近似白噪声的性质）为输入信号；
（2）增加数据长度去L ；
（3）减小噪声信号v(k)的方差。
4.1 广义最小二乘递推算法的的MATLAB 仿真结果及分析
(1)、输入选用题目给出的30个随机数，即数据长度去L=30，噪声选用均值 零，方差分别为0.5、0.1和0.01的随机序列，辨识结果如表表4-1。
表中给出了三种情况下辨识参数结果即表中的估计值，估计值与真值的相对误差
（2）、输入均采用M 序列，噪声选择均值为零，方差为0.5、0.1和0.01的随机序列，辨识步长均为300步，辨识结果如表4-2。
表中给出了三种情况下辨识参数结果即表中的估计值，估计值与真值的相对误差.
（3）数据结果分析：输入采用M 序列比采用随机序列得到的辨识效果更好。噪声均值相等时，方差越大，辨识效果越差，反之，方差越小辨识效果越好。可以通过增加步长的方法提高辨识精度。
下面给出以M 序列作为输入，噪声均值为零，方差为0.01的随机序列，数据长度取L=30，得到的变化曲线图：
下面给出以M 序列作为输入，噪声均值为零，方差为0.01的不相关随机序列，数据长度去L=300，得到的变化曲线图：
参考文献
[1] 李言俊，张科，系统辨识理论及应用，国防工业出版社，2006年 [2]方崇智，萧德云，过程辨识，清华大学出版社，2002年
[3]贾秋玲，袁冬莉，栾云凤，基于MATLAB7.x/Simulink/Stateflow系统仿真、分析及设计，西北工业大学出版社，2006年
[4] 侯媛彬, 汪梅, 王立琦, 系统辨识及其MATLAB 仿真, 科学出版社,2004年
附录
广义最小二乘递推算法的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.376*z(k-1)-0.483*z(k-2)+57*u(k-1)+0.42*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('误差曲线') %图形标题
考虑仿真对象z(k)+1.5z(k-1)+0.7z(k-2)=u(k-1)+0.5u(k-2)+v(k)，
其中v(k)是服从正态分布的白噪声N(0,1)。输入信号采用4阶M 序列。
选择如下形式的辨识模型
z(k)+a1z(k-1)+a2z(k-2)=b1u(k-1)+b2u(k-2)+v(k)，
试利用递推最小二乘法辨识参数a1、a2、b1、b2。
Np=15;r=4;
X1=1;X2=1;X3=1;X4=1;
m_length = r*Np;
a=1;
for i=1:1:m_length
Y4=X4;Y3=X3;Y2=X2;Y1=X1;
X4=Y3;X3=Y2;X2=Y1;
X1=xor(Y3,Y4);
if Y4==0
M(i)=-a;
else
M(i)=a;
end
end
figure;
i=1:1:m_length;
plot(i,M);
% 白噪声
noise = zeros(1,m_length);
for i=1:1:m_length
temp = noise + 0.5*rands(1,m_length);
noise = temp;
end
noise = noise/12;
%noise = temp;
% parameter of system
n=2;d=1;a1=-1;a2=0.5;b1=1;b2=0.5;
S_U0=0.2;S_Y0=0.2;
% generate u,y
u0=ones(1,m_length)*S_U0;
U = M + u0 + noise;
figure;
i=1:1:m_length;
plot(i,U);
%formulation
y(1)=0;y(2)=0;y(3)=0;
Y(1)=S_Y0+y(1)+noise(1);Y(2)=S_Y0+y(2)+noise(2);Y(3)=S_Y0+y(3) +noise(3); for k=4:m_length
y(k) = b1*U(k-1-d)+b2*U(k-2-d)-a1*y(k-1)-a2*y(k-2);
Y(k)=S_Y0+y(k)+noise(k);
end
figure;
i=1:1:m_length;
plot(i,Y);
%辨识
% premitive value
c=2;
P = (c^3)*eye(3*n+d);
sita(:,3) = [0,0,0,0,0,0,0]';
alf = 0.995;
% M
%compute U0,Y0
sum_U = 0;sum_Y = 0;
for k=1:1:m_length
sum_U = sum_U + U(k);
sum_Y = sum_Y + Y(k);
end
t_U0 = sum_U/m_length;t_Y0 = sum_Y/m_length;
for k=1:1:m_length
t_u(k) = U(k) - t_U0;
t_y(k) = Y(k) - t_Y0;
end
figure;
i=1:1:m_length;
plot(i,t_u);
figure;
i=1:1:m_length;
plot(i,t_y);
v(1)=0;v(2)=0;v(3)=0;
for k=4:1:m_length
fai = [-t_y(k-1),-t_y(k-2),t_u(k-1-d),t_u(k-2-d),v(k-1),v(k-2),v(k-3)]';
v(k) = t_y(k) - fai'*sita(:,k-1);
W = P*fai;
dan = 1/(alf + fai'*W);
sita(:,k) = sita(:,k-1) + dan*W*v(k);
P = ( P - dan * W*W')/alf;
end
figure;
i = 3:1:m_length;
plot(i,sita(1,3:m_length),'r',i,sita(2,3:m_length),'g',i,sita(3,3:m_length),'b',i,sita(4,3:m_length),'y',i,sita(5,3:m_length),'m',i,sita(6,3:m_length),'c',i,sita(7,3:m_length),'k');
最小二乘法参数辨识
201403027
摘 要:系统辨识在工程中的应用非常广泛, 系统辨识的方法有很多种, 最小
二乘法是一种应用极其广泛的系统辨识方法. 阐述了动态系统模型的建立及其最小二乘法在系统辨识中的应用, 并通过实例分析说明了最小二乘法应用于系统辨识中的重要意义.
关键词:最小二乘法; 系统辨识; 动态系统
Abstract : System identification in engineering is widely used, system identification methods there are many ways， least squares method is a very wide range of application of system identification method and the least squares method elaborated establish a dynamic system models in System Identification applications and examples analyzed by the least squares method is applied to illustrate the importance of system identification.
Keywords : Least Squares; system identification; dynamic system
引言
随着科学技术的不断发展, 人们认识自然、利用自然的能力越来越强, 对于未知对象的探索也越来越深入. 我们所研究的对象, 可以依据对其了解的程度分为三种类型:白箱、灰箱和黑箱. 如果我们对于研究对象的内部结构、内部机制了解很深入的话, 这样的研究对象通常称之为“白箱”; 而有的研究对象, 我们对于其内部结构、机制只了解一部分, 对于其内部运行规律并不十分清楚, 这样的研究对象通常称之为“灰箱”; 如果我们对于研究对象的内部结构、内部机制及运行规律均一无所知的话, 则把这样的研究对象称之为“黑箱”. 研究灰箱和黑箱时, 将研究的对象看作是一个系统, 通过建立该系统的模型, 对模型参数进行辨识来确定该系统的运行规律. 对于动态系统辨识的方法有很多, 但其中应用最广泛, 辨识 效果良好的就是最小二乘辨识方法, 研究最小二乘法在系统辨识中的应用具有现实的、广泛的意义.
1．1 系统辨识简介
系统辨识是根据系统的输入输出时间函数来确定描述系统行为的数学模型。现代控制理论中的一个分支。通过辨识建立数学模型的目的是估计表征系统行为的重要参数，建立一个能模仿真实系统行为的模型，用当前可测量的系统的输入和输出预测系统输出的未来演变，以及设计控制器。对系统进行分析的主要问题是根据输入时间函数和系统的特性来确定输出信号。对系统进行控制的主要问题是根据系统的特性设计控制输入，使输出满足预先规定的要求。而系统辨识所研究的问题恰好是这些问题的逆问题。通常，预先给定一个模型类μ={M}（即给定一类已知结构的模型），一类输入信号u 和等价准则J=L(y，yM)(一般情况下，J 是误差函数，是过程输出y 和模型输出yM 的一个泛函) ；然后选择使误差函数J 达到最小的模型，作为辨识所要求的结果。系统辨识包括两个方面：结构辨识和参数估计。在实际的辨识过程中，随着使用的方法不同，结构辨识和参数估计这两个方面并不是截然分开的，而是可以交织在一起进行的。
1.2系统辨识的目的
在提出和解决一个辨识问题时，明确最终使用模型的目的是至关重要的。它对模型类（模型结构）、输入信号和等价准则的选择都有很大的影响。通过辨识建立数学模型通常有四个目的。
①估计具有特定物理意义的参数 有些表征系统行为的重要参数是难以直接测量的，例如在生理、生态、环境、经济等系统中就常有这种情况。这就需要通过能观测到的输入输出数据，用辨识的方法去估计那些参数。
②仿真 仿真的核心是要建立一个能模仿真实系统行为的模型。用于系统分析的仿真模型要求能真实反映系统的特性。用于系统设计的仿真，则强调设计参数能正确地符合它本身的物理意义。
③预测 这是辨识的一个重要应用方面，其目的是用迄今为止系统的可测量的输入和输出去预测系统输出的未来的演变。例如最常见的气象预报，洪水预报，其他如太阳黑子预报，市场价格的预测，河流污染物含量的预测等。预测模型辨识的等价准则主要是使预测误差平方和最小。只要预测误差小就是好的预测
模型，对模型的结构及参数则很少再有其他要求。这时辨识的准则和模型应用的目的是一致的，因此可以得到较好的预测模型。
④控制 为了设计控制系统就需要知道描述系统动态特性的数学模型，建立这些模型的目的在于设计控制器。建立什么样的模型合适，取决于设计的方法和准备采用的控制策略。
2最小二乘方法
2.1.1 系统辨识最小二乘方法简介
对工程实践中测得的数据进行理论分析，用恰当的函数去模拟数据原型是一类十分重要的问题，最常用的逼近原则是让实测数据和估计数据之间的距离平方和最小，这即是最小二乘法。最小二乘法是一种经典的数据处理方法。在系统辨识领域中 ,最小二乘法是一种得到广泛应用的估计方法 ,可用于动态系统 ,静态系统 , 线性系统 ,非线性系统。可用于离线估计，也可用于在线估计。这种辨识方法主要用于在线辨识。在随机的环境下，利用最小二乘法时，并不要求观测数据提供其概率统计方面的信息，而其估计结果，却有相当好的统计特性。
MATLAB 是一套高性能数字计算和可视化软件 , 它集成概念设计 , 算法开发 , 建模仿真 , 实时实现于一体 , 构成了一个使用方便、界面友好的用户环境 , 其强大的扩展功能为各领域的应用提供了基础。对于比较复杂的生产过程 , 由于过程的输入输出信号一般总是可以测量的 , 而且过程的动态特性必然表现在这些输入输出数据中 , 那么就可以利用输入输出数据所提供的信息来建立过程的数学模型。这种建模方法就称为系统辨识。把辨识建模称作“黑箱建模”。
2.1.2 最小二乘法系统辨识结构：
本文把待辨识的过程看作“黑箱”。只考虑过程的输入输出特性，而不强调过程的内部机理。
图中，输入u(k)和输出z(k)是可以观测的；G 是系统模型，用来描述系统的输入输出特性；N 是噪声模型，v(k)是白噪声，e(k)是有色噪声，根据表示定理： 可以表示为
e(k) =N v(k)
G (z
-1
)
B (z -1) =
A (z -1)
N (z
-1
)
D (z -1) = -1
C (z )
-n a -1-1-2
⎧A (z ) =1+a z +a z ++a z 12n a ⎪⎨-n b -1-1-2B (z ) =b z +b z ++b z ⎪12n b ⎩
-n a -1-1-2⎧C (z ) =1+c z +c z + +c z ⎪12n a ⎨-n b -1-1-2D (z ) =d z +d z + +d z ⎪12n b ⎩
2.1.3 准则函数
设一个随机序列{z (k ), k ∈(1, 2, , L ) }的均值是参数θ的线性函数： E {z (k ) }=h T (k ) θ，其中h (k ) 是可测的数据向量，那么利用随机序列的一个实现，使准则函数：
J (θ) =∑[z (k ) -h T (k ) θ]
k =1
∧
L
2
(式2-2)
达到极小的参数估计值θ称作θ的最小二乘估计。
t
z (k ) =h (k ) θ+e (k ) ，θ为模型参数向量，e (k )为零均值随机 最小二乘格式：
2.2 广义最小二乘法
2.2.1 广义最小二乘数学模型
A (z -1) z (k ) =B (z -1) u (k ) +
1
v (k ) -1
C (z )
式中，u(k)和z (k ) 表示系统的输入输出；v(k)是均值为零的不相关的随机序列；且
⎧A (z -1) =1+a 1z -1+a 2z -2+ +a n a z -n a ⎪-n b -1-1-2B (z ) =b z +b z + +b z 12n ⎨b ⎪-n c -1-1-2C (z ) =1+c z +c z + +c z 12n c ⎩
2.2.2 广义最小二乘递推算法如下
ˆ(k ) =θˆ(k -1) +K (k )[z (k ) -h (k ) θˆ(k -1)]⎧θ⎪
⎪K (k ) =P (k -1) h (k ) [h (k ) P (k -1) h (k ) +1]⎪⎪P (k ) =[I -K (k ) h (k )]P (k -1) ⎨ˆˆˆ(k -1)]ˆ(k ) -h (k ) θ⎪θ(k ) =θ(k -1) +K (k )[e
⎪K (k ) =P (k -1) h (k ) [h (k ) P (k -1) h (k ) +1]⎪⎪⎩P (k ) =[I -K (k ) h (k )]P (k -1)
τ
f
f
f
τ
-1
f f f f f f
τ
f f f f
τ
e e e e e
τ
-1
e e e e e e
τ
e e e e
式中
⎧h (k ) =[-z (k -1), , -z (k -n ), u (k -1), , u (k -n )]⎪h (k ) =[-e ˆ(k -1), , -e ˆ(k -n )]⎨⎪e
⎩ˆ(k ) =z (k ) -h (k ) θˆ(k )
f
f
f
a τ
f
f
b
e
c
τ
τ
2.2.3 广义最小二乘递推算法的计算步骤：
ˆ(0) =ε(充分小的实向量) ⎧θ⎪
⎪P (0) =a I (a 为充分大的数)
1. 给定初始条件 ⎨
ˆ(0) =0⎪θ
⎪P (0) =I ⎩
z (k ) =C (z) z (k )
2利用式，计算z (k ) 和u (k ) ；
u (k ) =C (z) z (k )
2
f e e -1
f
-1
f f
f
⎧θ=[a , , a , b , , b ]
3利用式⎨，构造h (k ) ；
⎩h (k ) =[-z (k -1), , -z (k -n ), u (k -1), , u (k -n )]
τ
1
n a
1
n b
τf
f f f a f f b
ˆ(k ) =θˆ(k -1) +K (k )[z (k ) -h (k ) θˆ(k -1)]⎧θ⎪⎪ˆ(k ) ； 4利用式⎨K (k ) =P (k -1) h (k ) [h (k ) P (k -1) h (k ) +1]递推计算θ⎪P (k ) =[I -K (k ) h (k )]P (k -1) ⎪⎩
τ
f
f
f
τ
-1
f f f f f f
τ
f f f f
ˆ(k ) 和 ˆ(k ) =z (k ) -h (k ) θ5利用e
τ
ˆ(k ) ； h (k ) =[-z (k -1), , -z (k -n ), u (k -1), , u (k -n )]计算e
τ
a
b
ˆ(k -1), , -e ˆ(k -n c )]来构造h e (k ) ； 6根据h e (k ) =[-e
τ
ˆ(k ) =θˆ(k -1) +K (k )[e ˆ(k -1)]⎧θˆ(k ) -h (k ) θ⎪
7利用⎨K (k ) =P (k -1) h (k ) [h (k ) P (k -1) h (k ) +1]
⎪P (k ) =[I -K (k ) h (k )]P (k -1) ⎩
τ
e
e
e
e
e
τ
-1
e e e e e e
τ
e e e e
返回第2步进行迭代计算，直至获得满意的辨识结果。
3 工程实例 3.1典型系统建模
以某微循环流体系统模型的参数辨识为例. 我们已经得到该系统模型的差分方程形式, 取特定点的压力波作为模型的输入, 以另一点的压力波作为模型的输出. 由于我们采集的数据是实时的, 因此采用在线辨识方法. 由于建立的微循环流体系统模型是一个单输入、单输出的模型, 为使参数估计的结果很好地跟踪参数真值的变化, 我们采用渐消记忆的最小二乘法对系统模型参数进行辨识, 即强调新数据的作用, 贬低老数据的作用. 抽象出的SISO 系统的差分方程为：
（1-1） z (k ) +a 1z (k -1) +a 2z (k -2) =b 1u (k -1) +b 2u (k -2) +υ(k ) 式
T
θ 0.483 0.57 0.42]，利用MATLAB 的M 语言辨识参数取真值为：=[1.376
系统中的未知参数a 1、a 2、b 1、b 2。要求：用参数的真值利用差分方程求出z (k )
作为测量值，υ(k ) 是均值为0，方差为0.1、0.5和0.01的不相关随机序列。使用最小二乘算法辨识。
3.2 广义最小二乘递推算法的MATLAB 仿真（程序源代码见附录） 考虑仿真对象
z(k)= -1.376z(k-1)-0.483z(k-2)+0.57u(k-1)+0.42u(k-2)+v(k)
式中，v(k) 是均值为0，方差为0.01、0.1和0.5的不相关随机序列。输入信号采用4阶M 序列，幅度为1。 选择如下形式的辨识模型
其中取c1=0,c2=0.
4结果分析及算法优化
由于辨识算法中输入或噪声信号为不相关随机序列，所以每次辨识结果都不完全相同。但是，在相同输入、相同的噪声、相同的步长条件下，精度大体相同。
算法优化方案：（1）使用M 序列（具有近似白噪声的性质）为输入信号；
（2）增加数据长度去L ；
（3）减小噪声信号v(k)的方差。
4.1 广义最小二乘递推算法的的MATLAB 仿真结果及分析
(1)、输入选用题目给出的30个随机数，即数据长度去L=30，噪声选用均值 零，方差分别为0.5、0.1和0.01的随机序列，辨识结果如表表4-1。
表中给出了三种情况下辨识参数结果即表中的估计值，估计值与真值的相对误差
（2）、输入均采用M 序列，噪声选择均值为零，方差为0.5、0.1和0.01的随机序列，辨识步长均为300步，辨识结果如表4-2。
表中给出了三种情况下辨识参数结果即表中的估计值，估计值与真值的相对误差.
（3）数据结果分析：输入采用M 序列比采用随机序列得到的辨识效果更好。噪声均值相等时，方差越大，辨识效果越差，反之，方差越小辨识效果越好。可以通过增加步长的方法提高辨识精度。
下面给出以M 序列作为输入，噪声均值为零，方差为0.01的随机序列，数据长度取L=30，得到的变化曲线图：
下面给出以M 序列作为输入，噪声均值为零，方差为0.01的不相关随机序列，数据长度去L=300，得到的变化曲线图：
参考文献
[1] 李言俊，张科，系统辨识理论及应用，国防工业出版社，2006年 [2]方崇智，萧德云，过程辨识，清华大学出版社，2002年
[3]贾秋玲，袁冬莉，栾云凤，基于MATLAB7.x/Simulink/Stateflow系统仿真、分析及设计，西北工业大学出版社，2006年
[4] 侯媛彬, 汪梅, 王立琦, 系统辨识及其MATLAB 仿真, 科学出版社,2004年
附录
广义最小二乘递推算法的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.376*z(k-1)-0.483*z(k-2)+57*u(k-1)+0.42*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('误差曲线') %图形标题
考虑仿真对象z(k)+1.5z(k-1)+0.7z(k-2)=u(k-1)+0.5u(k-2)+v(k)，
其中v(k)是服从正态分布的白噪声N(0,1)。输入信号采用4阶M 序列。
选择如下形式的辨识模型
z(k)+a1z(k-1)+a2z(k-2)=b1u(k-1)+b2u(k-2)+v(k)，
试利用递推最小二乘法辨识参数a1、a2、b1、b2。
Np=15;r=4;
X1=1;X2=1;X3=1;X4=1;
m_length = r*Np;
a=1;
for i=1:1:m_length
Y4=X4;Y3=X3;Y2=X2;Y1=X1;
X4=Y3;X3=Y2;X2=Y1;
X1=xor(Y3,Y4);
if Y4==0
M(i)=-a;
else
M(i)=a;
end
end
figure;
i=1:1:m_length;
plot(i,M);
% 白噪声
noise = zeros(1,m_length);
for i=1:1:m_length
temp = noise + 0.5*rands(1,m_length);
noise = temp;
end
noise = noise/12;
%noise = temp;
% parameter of system
n=2;d=1;a1=-1;a2=0.5;b1=1;b2=0.5;
S_U0=0.2;S_Y0=0.2;
% generate u,y
u0=ones(1,m_length)*S_U0;
U = M + u0 + noise;
figure;
i=1:1:m_length;
plot(i,U);
%formulation
y(1)=0;y(2)=0;y(3)=0;
Y(1)=S_Y0+y(1)+noise(1);Y(2)=S_Y0+y(2)+noise(2);Y(3)=S_Y0+y(3) +noise(3); for k=4:m_length
y(k) = b1*U(k-1-d)+b2*U(k-2-d)-a1*y(k-1)-a2*y(k-2);
Y(k)=S_Y0+y(k)+noise(k);
end
figure;
i=1:1:m_length;
plot(i,Y);
%辨识
% premitive value
c=2;
P = (c^3)*eye(3*n+d);
sita(:,3) = [0,0,0,0,0,0,0]';
alf = 0.995;
% M
%compute U0,Y0
sum_U = 0;sum_Y = 0;
for k=1:1:m_length
sum_U = sum_U + U(k);
sum_Y = sum_Y + Y(k);
end
t_U0 = sum_U/m_length;t_Y0 = sum_Y/m_length;
for k=1:1:m_length
t_u(k) = U(k) - t_U0;
t_y(k) = Y(k) - t_Y0;
end
figure;
i=1:1:m_length;
plot(i,t_u);
figure;
i=1:1:m_length;
plot(i,t_y);
v(1)=0;v(2)=0;v(3)=0;
for k=4:1:m_length
fai = [-t_y(k-1),-t_y(k-2),t_u(k-1-d),t_u(k-2-d),v(k-1),v(k-2),v(k-3)]';
v(k) = t_y(k) - fai'*sita(:,k-1);
W = P*fai;
dan = 1/(alf + fai'*W);
sita(:,k) = sita(:,k-1) + dan*W*v(k);
P = ( P - dan * W*W')/alf;
end
figure;
i = 3:1:m_length;
plot(i,sita(1,3:m_length),'r',i,sita(2,3:m_length),'g',i,sita(3,3:m_length),'b',i,sita(4,3:m_length),'y',i,sita(5,3:m_length),'m',i,sita(6,3:m_length),'c',i,sita(7,3:m_length),'k');