方法一、最小二乘一次性算法：首先对最小二乘法的一次性辨识算法做简要介绍如下：过程的黑箱模型如图所示：其中u(k)和z(k)分别是过程的输入输出，)(1-z G 描述输入输出关系的模型，成为过程模型。

过程的输入输出关系可以描述成以下最小二乘格式： )()()(k n k h k z T +=θ （1）其中z(k)为系统输出，θ是待辨识的参数，h(k)是观测数据向量，n(k)是均值为0的随机噪声。

利用数据序列{z （k ）}和{h （k ）}极小化下列准则函数：∑=-=Lk T k h k z J 12])()([)(θθ （2）使J 最小的θ的估计值^θ，成为最小二乘估计值。

具体的对于时不变SISO 动态过程的数学模型为 )()()()()(11k n k u z B k z z A +=-- （3）应该利用过程的输入、输出数据确定)(1-z A 和)(1-Z B 的系数。

对于求解θ的估计值^θ，一般对模型的阶次a n ，b n 已定，且b a n n >；其次将（3）模型写成最小二乘格式 )()()(k n k h k z T +=θ （4）式中⎪⎩⎪⎨⎧=------=T n n T b a b a b b b a a a n k u k u n k z k z k h ],,,,,,,[)](,),1(),(,),1([)(2121 θ （5）L k ,,2,1 =因此结合式（4）（5）可以得到一个线性方程组L L L n H Z +=θ （6）其中⎩⎨⎧==T L TL L n n n n L z z z z )](),2(),1([)](),2(),1([ （7）对此可以分析得出，L H 矩阵的行数为),max (b a n n L -，列数b a n n +。

在过程的输入为2n 阶次，噪声为方差为1，均值为0的随机序列，数据长度)(b a n n L +>的情况下，取加权矩阵L Λ为正定的单位矩阵I ，可以得出：L T L L T L z H H H 1^)(-=θ （8） 其次，利用在Matlab 中编写M 文件，实现上述算法。

此次算法的实现，采用6阶M 序作为过程黑箱的输入；噪声采用方差为1，均值为0的随机数序列；黑箱模型假设为：y(k)-1.5y(k-1)+0.7y(k-2)=2u(k-1)+0.5u(k-2)，则系统输出为Z(k)-1.5Z(k-1)+0.7Z(k-2)=2U(k-1)+0.5U(k-2)+n （k ）；模型的阶次2,2==b a n n ；数据长度取L=200。

程序清单如下见附录：最小二乘一次性算法Matlab 程序运行结果如下：图1 最小二乘一次性算法参数真值与估计值其中re 为真值，ans 为估计值^θ结果发现辨识出的参数与真值之间存在细微误差，这是由于系统噪声以及数据长度L 的限制引起的，最小二乘辨识法是一种无偏估计方法。

方法二、最小二乘递推算法：最小二乘一次性算法计算量大，并且浪费存储空间，不利于在线应用，由此引出最小二乘参数估计的递推算法，其基本思想是按观测次序一步一修正，即：^^^'(1)()[()()(1)]k K k z k h k k θθθ=-+--'11()()(1)()[()(1)()]k K k P k h k h k P k h k -Λ=--+'()[()()](1)P k I K k h k P k =--这里选取估计参数的初值为全0矩阵，P 矩阵设置参数为100的对角阵，增益阵K 初值设置为10。

程序清单如下见附录：最小二乘递推算法Matlab 程序。

程序运行结果如下：图2 最小二乘递推算法参数真值与估计值图3 最小二乘递推算法待估参数过渡过程附录：1. 最小二乘一次性算法Matlab 程序%最小二乘一次性算法%Z(k)=1.5*Z(k-1)-0.7*Z(k-2)+2*u(k-1)+0.5u(k-2)+v(k)% na=2,nb=2clearclctic %计时开始L=200; %数据长度%生成M序列6阶a1=1;a2=0;a3=1;a4=0;a5=1;a6=0;M=zeros(L,1);for i=1:1:LM(i)=a1;a1=a2;a2=a3;a3=a4;a4=a5;a5=a6;a6=xor(M(i),a1);end%噪声v=randn(1,L);% na=2, 设置测量输出z(1)到z(na)的值为0z=[];z(1)=0;z(2)=0;%产生输出序列z 从z(na+1)开始,即z(3)开始for i=3:Lz(i)=1.5*z(i-1)-0.7*z(i-2)+2*M(i-1)+0.5*M(i-2)+v(i-2); end%构造HL矩阵na=2;nb=2; %模型阶次，并找出最大者赋给max if na>nbmax=na;elsemax=nb;end%HL矩阵的的维数hang=L-max;lie=na+nb;%构建HL矩阵用z和M填充HL矩阵与na,nb的最大数max有关HL=zeros(hang,lie);a=na;for j=1:1:nafor i=1:1:hangHL(i,j)=-z(i+a-1);enda=a-1;endb=nb;for j=na+1:1:liefor i=1:1:hangHL(i,j)=M(i+b-1);endb=b-1;end%构建ZL矩阵ZL=zeros(hang,1);for i=1:1:hang;ZL(i)=z(i+max);end%参数估计值th=zeros(na+nb,1);th=inv((HL'*HL))*HL'*ZL;th;%比较估计误差re为真值，th为估计值re=[-1.5,0.7,2,0.5]th'toc %计时结束2. 最小二乘递推算法Matlab程序%最小二乘递推算法%Z(k)=1.5*Z(k-1)-0.7*Z(k-2)+2*u(k-1)+0.5u(k-2)+v(k)% na=2,nb=2clearclctic %计时开始L=200; %数据长度%生成M序列6阶a1=1;a2=0;a3=1;a4=0;a5=1;a6=0;M=zeros(L,1);for i=1:1:LM(i)=a1;a1=a2;a2=a3;a3=a4;a4=a5;a5=a6;a6=xor(M(i),a1);end%噪声v=randn(1,L);% na=2, 设置测量输出z(1)到z(na)的值为0z=[];z(1)=0;z(2)=0;%产生输出序列z 从z(na+1)开始,即z(3)开始for i=3:Lz(i)=1.5*z(i-1)-0.7*z(i-2)+2*M(i-1)+0.5*M(i-2)+v(i-2); endna=2;nb=2; %模型阶次，并找出最大者赋给max if na>nbmax=na;elsemax=nb;end%HL矩阵的的维数hang=L-max;lie=na+nb;%构建HL矩阵用z和M填充HL矩阵与na,nb的最大数max有关HL=zeros(hang,lie);a=na;for j=1:1:nafor i=1:1:hangHL(i,j)=-z(i+a-1);enda=a-1;endb=nb;for j=na+1:1:liefor i=1:1:hangHL(i,j)=M(i+b-1);endb=b-1;end%P是(na+nb)*(na+nb)维，K是(na+nb)*1维，h'是1*(na+nb)维%递推求解P=100*eye(na+nb);Pstore=zeros(na+nb,hang);Pstore(:,1)=[P(1,1),P(2,2),P(3,3),P(4,4)];th=zeros(na+nb,hang);% 预估参数初值矩阵K初值th(:,1)=[0.001;0.001;0.001;0.001];K=[10;10;10;10];for i=1:hangh=HL(i,:)';K=P*h*inv(h'*P*h+1);th(:,i+1)=th(:,i)+K*(z(i+na)-h'*th(:,i));P=(eye(na+nb)-K*h')*P;Pstore(:,i+1)=[P(1,1),P(2,2),P(3,3),P(4,4)];endre=[-1.5,0.7,2,0.5]th(:,hang+1)'toc %计时结束i=1:hang+1;figure(1)plot(i,th(1,:),i,th(2,:),i,th(3,:),i,th(4,:))title('待估参数过渡矩阵')figure(2)plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:)) title('估计方差变化过程')。