(完整)系统辨识—最小二乘法汇总,推荐文档
(完整版)最小二乘法
的关系吗？
提示 对于单变量样本数据而言，平均数是样本数据的中 心，类似地，对于双变量样本数据，假设样本点为(x1，y1)， (x2，y2)，…，(xn，yn)，记－x ＝n1(x1＋x2＋…＋xn)，－y ＝n1(y1 ＋y2＋…＋yn)，则(－x ，－y )为样本点的中心，回归直线一定 过这一点．
课前探究学习
课堂讲练互动
i 1 2 3 4 5 6 7 8 9 10 xi 2 4 4 6 6 6 7 7 8 10 yi 0.9 1.4 1.6 2.0 2.1 1.9 1.8 2.1 2.2 2.3 xiyi 1.8 5.6 6.4 12 12.6 11.4 12.6 14.7 17.6 23
课前探究学习
课堂讲练互动
2．回归直线方程求解的方法步骤 根据最小二乘法的思想和公式，利用计算器或计算机，可 以方便地求出回归方程． 求线性回归方程的步骤： 第1步：列表xi，yi，xiyi；
第 2 步：计算－x ，－y ，
第3步：代入公式计算b，a的值； 第4步：写出回归方程y＝a＋bx. 利用回归直线对总体进行估计：
利用回归直线，我们可以进行预测．若回归方程为y＝bx
＋a，则x＝x0处的估计值为：y＝bx0＋a.
课前探究学习
课堂讲练互动
题型一 求线性回归的方程
【例1】 某地10户家庭的年收入和年饮食支出的统计资料如下 表：
年收入 x(万元)
2 4 4 6 6 6 7 7 8 10
年饮食支出 y(万元) 0.9 1.4 1.6 2.0 2.1 1.9 1.8 2.1 2.2 2.3
（x1－－x ）（y1－－y ）＋（x2－－x ）（y2－－y ）＋…＋（xn－－x ）（yn－－y ） （x1－－x ）2＋（x2－－x ）2＋…＋（xn－－x ）2
系统辨识相关分析最小二乘
相关分析法辨识系统单位脉冲响应1辨识原理对于下图示的单输入单输出线性系统，其输入输出的因果关系可用卷积公式描述。
公式为：0()()()y t g x t d λλλ∞=-⎰把变量t 换成t +τ，上式两边同乘以x (t )，取时间的平均值，得11lim()(+)(){lim()(+)}22TTTTT T x t y t dt g x t x t dt d TTτλτλλ∞--→∞→∞=-⎰⎰⎰即 0()()()x y x R g R d τστλλ∞=-⎰上式即为维纳-霍夫方程，其给出了输入的自相关函数，输入、输出的互相关函数及脉冲响应函数三者之间的关系。
令x (t )为白噪声信号，则其自相关函数为：()(), ()()x x R k R k τδττλδτλ=-=-代入维纳-霍夫方程得：()()()()xy x R g R d kg τλτλλτ∞=-=⎰则有：()()xy R g kττ=这样，只要记录x(t)、y(t)的值，并计算它们的互相关函数，即可求得脉冲响应函数g(τ)。
在系统有正常输入的情形下，辨识脉冲响应的原理图如下图所示。
2辨识过程2.1预备实验以二阶系统22()2G s s s ++=作为辨识对象。
在实验前首先要进行预备实验，以了解系统特性。
通过简单阶跃响应确定系统过度过程时间T s 大约为11s ，如下图所示。
给系统施加不同周期的正弦信号，系统输出为输入的0.707倍时，确定截止频率f M 大约为0.318Hz 。
2.2选择二位式伪随机序列的参数（1）选择t 和N 由0.3Mt f ∆≤，得0.94t s ∆≤。
因为系统的时间常数1nT s ζω=，根据时间常数可按照0.050.1t T ∆= （）选择t ∆。
由二位式伪随机序列周期要大于系统过渡过程时间，若t ∆选择0.94s ，则由(1)s N t T -⨯∆≥，得12.7021N ≥；若t ∆选择0.195s ，则由(1)s N t T -⨯∆≥，得57.4103N ≥。
系统辨识之最小二乘法
系统辨识之最小二乘法方法一、最小二乘一次性算法：首先对最小二乘法的一次性辨识算法做简要介绍如下：过程的黑箱模型如图所示：其中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 文件，实现上述算法。
最小二乘法辨识
T 1 T T 1 T ˆ E [ θ ] E [ θ ( Φ Φ ) Φ ξ ] θ E [( Φ Φ ) Φ ξ ]
LS无偏估计的充要条件为：
E [( Φ Φ )
T 1
Φ
T
ξ] 0
下面讨论无偏估计的充分条件。
y ( k ) a1 y ( k 1) a n y ( k n ) b 0 u ( k ) b1u ( k 1) b n u ( k n ) ( k )
最小二乘的最早思想： 未知量的最大可能的值是这样 一个数值，它是实际观测值和计算 值的差值的平方和达到最小的数值。
基本的最小二乘估计 解决问题：在模型阶次n已知的情况下，根据系 统的输入输出数据，估计出系统差分方程的各 项系数。 1.基于输入/输出数据的系统模型描述
SISO系统的差分方程为
x ( k ) a 1 x ( k 1) a n x ( k n ) b 0 u ( k ) b n u ( k n ) y (k ) x(k ) n(k )
(k ) n(k )
a n(k i)
i i 1
n
则当前输出为：
y ( k ) a1 y ( k 1) a n y ( k n ) b 0 u ( k ) b1u ( k 1) b n u ( k n ) ( k )
ˆ min J
下面我们推导θ估计值的计算方法。
J取得最小值，也即J为极值，则有：
J ˆ θ
0
ˆ T ˆ [ (Y Φ θ ) (Y Φ θ ) ] ˆ θ
0
T ˆ 2 Φ (Y Φ θ ) 0
最小二乘参数辨识方法
《系统辨识基础》第17讲要点第5章 最小二乘参数辨识方法5.9 最小二乘递推算法的逆问题辨识是在状态可测的情况下讨论模型的参数估计问题，滤波是在模型参数已知的情况下讨论状态估计问题，两者互为逆问题。
5.10 最小二乘递推算法的几种变形最小二乘递推算法有多种不同的变形，常用的有七种情况：① 基于数据所含的信息内容不同，对数据进行有选择性的加权； ② 在认为新近的数据更有价值的假设下，逐步丢弃过去的数据； ③ 只用有限长度的数据；④ 加权方式既考虑平均特性又考虑跟综能力； ⑤ 在不同的时刻，重调协方差阵P (k )； ⑥ 设法防止协方差阵P (k )趋于零； 5.10.1 选择性加权最小二乘法 把加权最小二乘递推算法改写成[]⎪⎪⎩⎪⎪⎨⎧--=+--=--+-=-)1()]()([)(1)()1()()()()1()()()]1(ˆ)()()[()1(ˆ)(ˆ1k k k k k k k k k k k k k k k z k k k P h K P h P h h P K h K τττθθθI ΛΛ算法中引进加权因子，其目的是便于考虑观测数据的可信度．选择不同的加权方式对算法的性质会有影响，下面是几种特殊的选择：① 一种有趣的情况是Λ()k 取得很大，在极限情况下，算法就退化成正交投影算法。
也就是说，当选择⎩⎨⎧=-≠-∞=0)()1()(,00)()1()(,)(k k k k k k k h P h h P h ττΛ 构成了正交投影算法⎪⎪⎩⎪⎪⎨⎧--=--=--+-=)1()]()([)()()1()()()1()()]1(ˆ)()()[()1(ˆ)(ˆk k k k k k k k k k k k k z k k k P h K P h P h h P K h K τττθθθI 算法初始值取P ()0=I 及 ()θε0=(任定值)，且当0)()1()(=-k k k h P h τ时，令K ()k =0。
系统辨识—最小二乘法_3
---------------------------------------------------------------最新资料推荐------------------------------------------------------系统辨识—最小二乘法最小二乘法参数辨识 1 引言系统辨识是根据系统的输入输出时间函数来确定描述系统行为的数学模型。
现代控制理论中的一个分支。
通过辨识建立数学模型的目的是估计表征系统行为的重要参数，建立一个能模仿真实系统行为的模型，用当前可测量的系统的输入和输出预测系统输出的未来演变，以及设计控制器。
对系统进行分析的主要问题是根据输入时间函数和系统的特性来确定输出信号。
对系统进行控制的主要问题是根据系统的特性设计控制输入，使输出满足预先规定的要求。
而系统辨识所研究的问题恰好是这些问题的逆问题。
通常，预先给定一个模型类={M}（即给定一类已知结构的模型），一类输入信号 u 和等价准则 J=L(y，yM)(一般情况下，J 是误差函数，是过程输出 y 和模型输出 yM 的一个泛函)；然后选择使误差函数J 达到最小的模型，作为辨识所要求的结果。
系统辨识包括两个方面：结构辨识和参数估计。
在实际的辨识过程中，随着使用的方法不同，结构辨识和参数估计这两个方面并不是截然分开的，而是可以交织在一起进行的。
2 系统辨识的目的在提出和解决一个辨识问题时，明确最终使1 / 17用模型的目的是至关重要的。
它对模型类（模型结构）、输入信号和等价准则的选择都有很大的影响。
通过辨识建立数学模型通常有四个目的。
①估计具有特定物理意义的参数有些表征系统行为的重要参数是难以直接测量的，例如在生理、生态、环境、经济等系统中就常有这种情况。
这就需要通过能观测到的输入输出数据，用辨识的方法去估计那些参数。
②仿真仿真的核心是要建立一个能模仿真实系统行为的模型。
用于系统分析的仿真模型要求能真实反映系统的特性。
系统辨识—最小二乘法汇总
最小二乘法参数辨识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 系统辨识简介系统辨识是根据系统的输入输出时间函数来确定描述系统行为的数学模型。
最小二乘法
y (k ) = x(k ) + v(k )
y (k ) +
x(k ) = y (k ) − v(k )
n i=0
∑
n
i =1
ai y (k − i) =
∑
biu ( k − i ) + v ( k ) +
∑
n
i =1
aiv (k − i)
ξ (k ) = v(k ) +
∑
n
i =1
a iv (k − i)
系统辨识——最小二 (n) L − y (1) u (n + 1) L y (n + 1) y (n + 2) − y (n + 1) L − y (2) u (n + 2) L = M M M M M y (n + N ) − y (n + N − 1) L − y ( N ) u (n + N ) L a1 a 2 u (1) M ξ (n + 1) u (2) an ξ (n + 2) + b0 M u ( N ) b1 ξ (n + N ) M bn
系统辨识——最小二乘算法 系统辨识——最小二乘算法
系统参数a 系统参数a
a 1.5 a1 a2
^
1
a1 = −1.5003
a2 = 0.7006
^
0.5
a1=-1.5, a2= 0.7
0 theta -0.5 -1 -1.5 -2 0
10
20
30 time steps
40
50
最小二乘法系统辨识
最小二乘法一次完成算法的MATLAB仿真一、二阶系统的最小二乘法一次完成算法的辨识程序及结果程序源代码：clear % 工作间清零u = [-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1]; % 系统辨识的输入信号为一个周期的M序列z = zeros( 1, 16); % 定义输出观测值的长度for k = 3 : 16% 用理想输出值作为观测值z(k) = 1.8*z(k-1) - 0.8*z(k-2) + 1.1*u(k-1) + 0.4*u(k-2);ZL( k, : ) = z(k); % 给样本矩阵ZL赋值Hl = [ -z(k-1), -z(k-2), u(k-1), u(k-2) ]; % 用理想输出值作为观测值HL( k, : ) = Hl;endsubplot(3,1,1) % 画三行一列图形窗口中的第一个图形stem(u) % 画出输入信号u的经线图形subplot(3,1,2) % 画三行一列图形窗口中的第二个图形i=1:1:16; % 横坐标范围是1到16，步长为1% 图形的横坐标是采样时刻i, 纵坐标是输出观测值z, 图形格式为连续曲线plot(i,z)subplot(3,1,3) % 画三行一列图形窗口中的第三个图形stem(z),grid on % 画出输出观测值z的经线图形，并显示坐标网格u,z % 显示输入信号和输出观测信号c1=HL'*HL; % 计算参数并显示c2=inv(c1);c3=HL'*ZL;c=c2*c3a1=c(1), a2=c(2), b1=c(3), b2=c(4) % 从中分离出并显示a1 、a2、b1、b2运算结果如下：u =-1 1 -1 1 1 1 1 -1 -1 -1 1 -1 -1 1 1z =Columns 1 through 90 0 0.7000 0.5600 1.1480 3.1184 6.1947 10.1558 12.6246 Columns 10 through 1613.0997 11.9798 11.7838 10.9270 8.7416 7.6933 8.3546c =-1.80000.80001.10000.4000a1 =-1.8000a2 =0.8000b1 =1.1000b2 =0.4000所画图形如下：二、实际压力系统的最小二乘辨识程序及结果程序源代码：clear % 工作间清零V = [ 54.3, 61.8, 72.4, 88.7, 118.6, 194.0]' % 赋初值V，并显示P = [ 61.2, 49.5, 37.6, 28.4, 19.2, 10.1]' % 赋初值P，并显示% P、V之间的关系:logP=-alpha*logV% +logbeita=[-logV,1][alpha,log(beita)]'=HL*sitafor i = 1 : 6; % 循环变量的取值为从1到6 Z(i) = log( P (i) ); % 赋系统的输出采样值ZL( i, : ) = Z(i); % 给样本矩阵ZL赋值Hl = [ log( V(i) ),1 ]; % 给HL赋值HL( i, : ) = Hlend % 循环结束c1=HL'*HL; % 计算参数并显示c2=inv(c1);c3=HL'*ZL;c4=c2*c3alpha = c4(1) % 为c4的第1个元素beita = exp(c4(2)) % 为以自然数为底的c4的第2个元素的指数运算结果如下：V =54.300061.800072.400088.7000118.6000194.0000P =61.200049.500037.600028.400019.200010.1000HL =3.9945 1.0000HL =3.9945 1.00004.1239 1.0000HL =3.9945 1.00004.1239 1.00004.2822 1.0000HL =3.9945 1.00004.1239 1.00004.2822 1.00004.4853 1.0000HL =3.9945 1.00004.1239 1.00004.2822 1.00004.4853 1.00004.7758 1.0000HL =3.9945 1.00004.1239 1.00004.2822 1.00004.4853 1.00004.7758 1.00005.2679 1.0000 c4 =-1.40429.6786alpha =-1.4042beita =1.5972e+004。
系统辨识各类最小二乘法汇总
yk(k)=1.5*yk(k-1)-0.7*yk(k-2)+uk(k-1)+0.5*uk(k-2)+y1(k); end figure(3); plot(yk); title('对应输出曲线');
theta=[0;0;0;0]; p=10^6*eye(4);
9
for t=3:N h=([-yk(t-1);-yk(t-2);uk(t-1);uk(t-2)]); x=1+h'*p*h; p=(p-p*h*1/x*h'*p); theta=theta+p*h*(yk(t)-h'*theta);
12
p=(p-p*h*1/x*h'*p); theta=theta+p*h*(yk(t)-h'*theta);
a1t(t)=theta(1); a2t(t)=theta(2); b1t(t)=theta(3); b2t(t)=theta(4); d1t(t)=theta(5); d2t(t)=theta(6);
end 5、RGLS 试验程序（部分） for t=3:N
he=([-e(t-1);-e(t-2)]); xe=1+he'*pe*he; pe=(pe-pe*he*1/xe*he'*pe); thete=thete+pe*he*(e(t)-he'*thete);
c1t(t)=thete(1); c2t(t)=thete(2);
7
RELS： 当噪声模型： e k = D Z −1 ∗ v k ( v(k) 为白噪声 )时，我们采用增广最 小二乘方法。能辨识出参数（包括噪声参数）的无偏估计。 RGLS： 当噪声模型： e k =
(完整word版)最小二乘法(word文档良心出品)
最小二乘法基本原理：成对等精度测得一组数据，试找出一条最佳的拟合曲线，使得这条曲线上的各点值与测量值的平方和在所有的曲线中最小。
我们用最小二乘法拟合三次多项式。
最小二乘法又称曲线拟合，所谓的“拟合”就是不要求曲线完全通过所有的数据点，只要求所得的曲线反映数据的基本趋势。
曲线的拟合几何解释：求一条曲线，使所有的数据均在离曲线的上下不远处。
第一节 最小二乘法的基本原理和多项式拟合 一最小二乘法的基本原理从整体上考虑近似函数)(x p 同所给数据点),(i i y x (i=0,1,…,m)误差i i i y x p r -=)((i=0,1,…,m)常用的方法有以下三种：一是误差i i i y x p r -=)((i=0,1,…,m)绝对值的最大值im i r ≤≤0max ，即误差 向量T m r r r r ),,(10 =的∞—范数；二是误差绝对值的和∑=mi ir 0，即误差向量r 的1—范数；三是误差平方和∑=mi ir02的算术平方根，即误差向量r 的2—范数；前两种方法简单、自然，但不便于微分运算 ，后一种方法相当于考虑 2—范数的平方，因此在曲线拟合中常采用误差平方和∑=mi ir02来 度量误差i r (i=0，1，…，m)的整体大小。
数据拟合的具体作法是：对给定数据 ),(i i y x (i=0,1,…，m)，在取定的函数类Φ中,求Φ∈)(x p ,使误差i i i y x p r -=)((i=0,1,…,m)的平方和最小，即∑=m i ir 02=[]∑==-mi ii y x p 02min)(从几何意义上讲，就是寻求与给定点),(i i y x (i=0,1,…,m)的距离平方和为最小的曲线)(x p y =（图6-1）。
函数)(x p 称为拟合 函数或最小二乘解，求拟合函数)(x p 的方法称为曲线拟合的最小二乘法。
Φ可有不同的选取方法.6—1二多项式拟合假设给定数据点),(i i y x (i=0,1,…,m)，Φ为所有次数不超过)(m n n ≤的多项式构成的函数类，现求一Φ∈=∑=nk k k n x a x p 0)(,使得[]min )(00202=⎪⎭⎫⎝⎛-=-=∑∑∑===mi mi n k i k i k i i n y x a y x p I (1)当拟合函数为多项式时，称为多项式拟合，满足式（1）的)(x p n 称为最小二乘拟合多项式。
系统辨识之最小二乘
写成矩阵的形式为：
Y = Fq + e
3.6
é - y(n)
- y(n -1) … - y(1) u(n) … u(1) ù
F
=
ê ê ê
-
y(n "
-1)
- y(n) "
… - y(2)
u(n +1)
…
u(2)
ú ú
…"
" # "ú
êë- y(n + N ) - y(n + N -1) ! - y(N ) u(n + N ) ! u(N )úû
q = F -1y
如果噪声x ¹ 0，则 q = F-1y - F-1x
从上式可以看出噪声 x 对参数估计有影响，为了尽量减小噪声 x 对q 估值的影响，应取 N>(2n+1), 即方程数大于未知数数目。在这种情况下，不能用解方程的办法来求q ，而要采用数理统计的办法， 以便减小噪声对q 估值的影响。在给定输出向量 y 和测量矩阵 F 的条件下求系统参数q 的估值，这 就是系统辨识问题。可用最小二乘法来就q 的估值。
在科学研究中，为了揭示某些相关量之间的关系，找出其规律，往往需要求解其函数解析式。
一种方法是采用插值逼近法，即所构造的近似函数 f (x) 在已知节点 xi 上必须满足 f(xi) = yi 要求逼 近函数 f (xi) 与被逼近函数 f (x) 在各已知点 xi 处的误差为零，即要求 f (x) 的曲线必须通过所有的
取泛函 J (q )为
N
å å J (q ) = (Y - Fq )2 = e2 (n + i) = eT • e = (Y - Fq )T (Y - Fq ) i =1
(完整word版)2003版系统辨识最小二乘法大作业
西北工业大学系统辩识大作业题目：最小二乘法系统辨识一、 问题重述：用递推最小二乘法、加权最小二乘法、遗忘因子法、增广最小二乘法、广义最小二乘法、辅助变量法辨识如下模型的参数离散化有z^4 - 3.935 z^3 + 5.806 z^2 - 3.807 z + 0.9362---------------------------------------------- =z^4 - 3.808 z^3 + 5.434 z^2 - 3.445 z + 0.8187噪声的成形滤波器离散化有4.004e-010 z^3 + 4.232e-009 z^2 + 4.066e-009 z + 3.551e-010-----------------------------------------------------------------------------=z^4 - 3.808 z^3 + 5.434 z^2 - 3.445 z + 0.8187采样时间0.01s要求：1.用Matlab 写出程序代码；2.画出实际模型和辨识得到模型的误差曲线；3.画出递推算法迭代时各辨识参数的变化曲线；最小二乘法:在系统辨识领域中 ,最小二乘法是一种得到广泛应用的估计方法 ,可用于动态 ,静态 , 线性 ,非线性系统。
在使用最小二乘法进行参数估计时 ,为了实现实时控制 ,必须优化成参数递推算法 ,即最小二乘递推算法。
这种辨识方法主要用于在线辨识。
MATLAB 是一套高性能数字计算和可视化软件 ,它集成概念设计 ,算法开发 ,建模仿真 ,实时实现于一体 ,构成了一个使用方便、界面友好的用户环境 ,其强大的扩展功能为各领域的应用提供了基础。
对于一个简单的系统 ,可以通过分析其过程的运动规律 ,应用一些已知的定理和原理,建立数学模型 ,即所谓的“白箱建模 ”。
但对于比较复杂的生产过程 ,该建模方法有很大的局限性。
由于过程的输入输出信号一般总是可以测量的 ,而且过程的动态特性必然表现在这些输入输出数据中 ,那么就可以利用输入输出数据所提供的信息来建立过程的数学模型。
系统辨识-最小二乘法
5.3 Matlab源程序：%最小二乘估计clearu=[ 1.147 0.201 -0.787 -1.589 -1.052 0.866 1.152 1.573 0.6260.433 -0.985 0.810 -0.044 0.947 -1.474 -0.719 -0.086 -1.099 1.4501.151 0.485 1.633 0.043 1.326 1.706 -0.340 0.890 0.144 1.177-0.390];n=normrnd(0, sqrt(0.1), 1, 31); %方差为0.1z=zeros(1,30);for k=3:31z(k)=-1.642*z(k-1)-0.715*z(k-2)+0.39*u(k-1)+0.35*u(k-2)+n(k)+1.642*n(k-1)+0.715 *n(k-2);endh0=[-z(2) -z(1) u(2) u(1)]';HLT=[h0,zeros(4,28)];for k=3:30h1=[-z(k) -z(k-1) u(k) u(k-1)]';HLT(:,k-1)=h1;endHL=HLT';y=[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)];%求出FAIc1=HL'*HL;c2=inv(c1);c3=HL'*y;c=c2*c3;%display('方差=0.1时，最小二乘法估计辨识参数θ如下：');a1=c(1)a2=c(2)b1=c(3)b2=c(4)%递推最小二乘法估计clearu=[ 1.147 0.201 -0.787 -1.589 -1.052 0.866 1.152 1.573 0.6260.433 -0.985 0.810 -0.044 0.947 -1.474 -0.719 -0.086 -1.099 1.4501.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;n=normrnd(0, sqrt(0.1), 1, 31); %方差为0.1for k=3:31z(k)=-1.642*z(k-1)-0.715*z(k-2)+0.39*u(k-1)+0.35*u(k-2)+n(k)+1.642*n(k-1)+0.715 *n(k-2);endc0=[0.001 0.001 0.001 0.001]'; %直接给出被辨识参数的初始值,即一个充分小的实向量p0=10^6*eye(4,4); %直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.000000005; %取相对误差E=0.000000005c=[c0,zeros(4,30)]; %被辨识参数矩阵的初始值及大小e=zeros(4,30); %相对误差的初始值及大小for k=3:30; %开始求Kh1=[-z(k-1),-z(k-2),u(k-1),u(k-2)]';x=h1'*p0*h1+1; x1=inv(x); %开始求K(k)k1=p0*h1*x1;%求出K的值d1=z(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; %如果参数收敛情况满足要求，终止计算endend%display('方差为0. 1递推最小二乘法辨识后的结果是：');a1=c(1,:);a2=c(2,:);b1=c(3,:);b2=c(4,:);%display('a1,a2,b1,b2经过递推最小二乘法辨识的结果是：');for i=3:31;if(c(1,i)==0)q1=c(1,i-1);break;endendfor i=3:31;if(c(2,i)==0)q2=c(2,i-1);break;endendfor i=3:31;if(c(3,i)==0)q3=c(3,i-1);break;endendfor i=3:31;if(c(4,i)==0)q4=c(4,i-1);break;endenda1=q1a2=q2b1=q3b2=q4%辅助变量递推最小二乘法估计clearna=2;nb=2;siitt=[1.642 0.715 0.39 0.35]';siit0=0.001*eye(na+nb,1);p=10^6*eye(na+nb);siit(:,1)=siit0;y(2)=0;y(1)=0;x(1)=0;x(2)=0;j=0;u=[ 1.147 0.201 -0.787 -1.589 -1.052 0.866 1.152 1.573 0.6260.433 -0.985 0.810 -0.044 0.947 -1.474 -0.719 -0.086 -1.099 1.4501.151 0.485 1.633 0.043 1.326 1.706 -0.340 0.890 0.144 1.177-0.390];n=normrnd(0, sqrt(0.1), 1, 31); %方差为0.1for k=3:31;h=[-y(k-1),-y(k-2),u(k-1),u(k-2)]';y(k)=h'*siitt+n(k)+1.642*n(k-1)+0.715*n(k-2);hx=[-x(k-1),-x(k-2),u(k-1),u(k-2)]';kk=p*hx/(h'*p*hx+1);p=[eye(na+nb)-kk*h']*p;siit(:,k-1)=siit0+kk*[y(k)-h'*siit0];x(k)=hx'*siit(:,k-1);j=j+(y(k)-h'*[1.642 0.715 0.39 0.35]')^2;e=max(abs((siit(:,k-1)-siit0)./siit0));ess(:,k-2)=siit(:,k-1)-siitt;siit0=siit(:,k-1);enda1=siit0(1)a2=siit0(2)b1=siit0(3)b2=siit0(4)clear%广义最小二乘估计clear;nn = normrnd(0,sqrt(0.1),1,31)'; %方差为0.1uk=[1.147 0.201 -0.787 -1.589 -1.052 0.866 1.152 1.573 0.626 0.433 -0.958 0.810 -0.044 0.947 -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 1.144 1.177 -0.390];yk(1)=0;yk(2)=0;for i=1:29;yk(i+2)=-1.642*yk(i+1)-0.715*yk(i)+0.39*uk(i+1)+0.35*uk(i)+nn(i+2)+1.642*nn(i+1 )+0.715*nn(i);end;for i=1:29;A(i,:)=[-yk(i+1) -yk(i) uk(i+1) uk(i)];endsiit=inv(A'*A)*A'*(yk(3:31)+nn(2:30)')';e(1)=yk(1);e(2)=yk(2)+siit(1)*yk(1)-siit(3)*uk(1);for i=3:31;e(i)=yk(i)+siit(1)*yk(i-1)+siit(2)*yk(i-2)-siit(3)*uk(i-1)-siit(4)*uk(i-2); endfor i=1:29;fai(i,:)=[-e(i+1) -e(i)];endf=inv(fai'*fai)*fai'*e(3:31)';for i=3:31;yk(i)=yk(i)+f(1)*yk(i-1)+f(2)*yk(i-2);endyk(2)=yk(2)+f(1)*yk(1);for i=3:30;uk(i)=uk(i)+f(1)*uk(i-1)+f(2)*uk(i-2);enduk(2)=uk(2)+f(1)*uk(1);for j=1:30for i=1:29;A(i,:)=[-yk(i+1) -yk(i) uk(i+1) uk(i)];endsiit=inv(A'*A)*A'*yk(3:31)';e(1)=yk(1);e(2)=yk(2)+siit(1)*(yk(1))-siit(3)*uk(1);for i=3:31;e(i)=yk(i)+siit(1)*(yk(i-1))+siit(2)*(yk(i-2))-siit(3)*uk(i-1)-siit(4)*uk(i-2); endfor i=1:29;fai(i,:)=[-e(i+1) -e(i)];endf=inv(fai'*fai)*fai'*e(3:31)';k1(j)=f(1);k2(j)=f(2);for i=3:31;yk(i)=yk(i)+f(1)*(yk(i-1))+f(2)*(yk(i-2));endyk(2)=yk(2)+f(1)*yk(1);for i=3:30uk(i)=uk(i)+f(1)*uk(i-1)+f(2)*uk(i-2);enduk(2)=uk(2)+f(1)*uk(1);endsiit'%增广矩阵估计参数clearu=[ 1.147 0.201 -0.787 -1.589 -1.052 0.866 1.152 1.573 0.6260.433 -0.985 0.810 -0.044 0.947 -1.474 -0.719 -0.086 -1.099 1.4501.151 0.485 1.633 0.043 1.326 1.706 -0.340 0.890 0.144 1.177-0.390];n=normrnd(0, sqrt(0.1), 1, 31); %方差为0.1z=zeros(7,30);zs=zeros(7,30);zm=zeros(7,30);zmd=zeros(7,30);%输出采样、不考虑噪声时系统输出、不考虑噪声时模型输出、模型输出矩阵的大小z(2)=0;z(1)=0;zs(2)=0;zs(1)=0;zm(2)=0;zm(1)=0;zmd(2)=0;zmd(1)=0;%给输出采样、不考虑噪声时系统输出、不考虑噪声时模型输出、模型输出赋初值for k=3:31z(k)=-1.642*z(k-1)-0.715*z(k-2)+0.39*u(k-1)+0.35*u(k-2)+n(k)+1.642*n(k-1)+0.715 *n(k-2);endc0=[0.001 0.001 0.001 0.001]'; %直接给出被辨识参数的初始值,即一个充分小的实向量p0=10^6*eye(4,4); %直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.000000005; %取相对误差E=0.000000005c=[c0,zeros(4,30)]; %被辨识参数矩阵的初始值及大小e=zeros(4,30); %相对误差的初始值及大小for k=3:30; %开始求Kh1=[-z(k-1),-z(k-2),u(k-1),u(k-2)]';x=h1'*p0*h1+1; x1=inv(x); %开始求K(k)k1=p0*h1*x1;%求出K的值d1=z(k)-h1'*c0; c1=c0+k1*d1; %求被辨识参数czs(k)=-1.642*z(k-1)-0.715*z(k-2)+0.39*u(k-1)+0.35*u(k-2);%系统在M序列的输入下不考虑扰动时的输出响应zm(k)=[-z(k-1),-z(k-2),u(k-1),u(k-2)]*[c1(1);c1(2);c1(3);c1(4)];%模型在输入下不考虑扰动时的的输出响应zmd(k)=h1'*c1;%模型在输入下的的输出响应e1=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; %如果参数收敛情况满足要求，终止计算endend%display('方差为0. 1递推最小二乘法辨识后的结果是：');a1=c(1,:);a2=c(2,:);b1=c(3,:);b2=c(4,:);%display('a1,a2,b1,b2经过递推最小二乘法辨识的结果是：');for i=3:31;if(c(1,i)==0)q1=c(1,i-1);break;endendfor i=3:31;if(c(2,i)==0)q2=c(2,i-1);break;endendfor i=3:31;if(c(3,i)==0)q3=c(3,i-1);break;endendfor i=3:31;if(c(4,i)==0)q4=c(4,i-1);break;endenda1=q1a2=q2b1=q3b2=q4%用夏氏偏差修正法估计参数clear;u=[ 1.147 0.201 -0.787 -1.589 -1.052 0.866 1.152 1.573 0.6260.433 -0.985 0.810 -0.044 0.947 -1.474 -0.719 -0.086 -1.099 1.4501.151 0.485 1.633 0.043 1.326 1.706 -0.340 0.890 0.144 1.177-0.390];n=normrnd(0, sqrt(0.1), 1, 31); %方差为0.1z=zeros(1,30);for k=3:31z(k)=-1.642*z(k-1)-0.715*z(k-2)+0.39*u(k-1)+0.35*u(k-2)+n(k)+1.642*n(k-1)+0.715 *n(k-2);endh0=[-z(2) -z(1) u(2) u(1)]';HLT=[h0,zeros(4,28)];for k=3:30h1=[-z(k) -z(k-1) u(k) u(k-1)]';HLT(:,k-1)=h1;endHL=HLT';y=[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)];%求出FAIc1=HL'*HL;c2=inv(c1);c3=HL'*y;c=c2*c3;fai1_xiashi=HL;siit1_ls_guzhi=c;ka_ma=inv(fai1_xiashi'*fai1_xiashi)*fai1_xiashi';M_xiashi=eye(29)-fai1_xiashi*ka_ma;siit_N0_xshi=siit1_ls_guzhi;siit_bxs=[10,10,10,10]';for n=1:30e=zeros(29,1);e(3)=y(5)+siit_N0_xshi(1,1)*y(4)+siit_N0_xshi(2,1)*y(3)-siit_N0_xshi(3,1)*u(4)-siit_N0_xshi(4,1)*u(3);e(4)=y(6)+siit_N0_xshi(1,1)*y(5)+siit_N0_xshi(2,1)*y(4)-siit_N0_xshi(3,1)*u(5)-siit_N0_xshi(4,1)*u(4);e(5:29) = y(5:29) - fai1_xiashi(5:29,:) * siit_N0_xshi;ou_mei_ka=zeros(29,1);for i=1:29ou_mei_ka(i,:)=-e(i);endD_xiashi=ou_mei_ka'*M_xiashi*ou_mei_ka;f_xiashi=inv(D_xiashi)*ou_mei_ka'*M_xiashi*y(1:29);siit_b=ka_ma*ou_mei_ka*f_xiashi;siit_xiashi=siit1_ls_guzhi-siit_b;if norm(siit_b-siit_bxs)/norm(siit_bxs)<0.0001break ;endsiit_N0_xshi=siit_xiashi;siit_bxs=siit_b;endsiit1_xiashi=siit_xiashi;%display('利用夏氏修正法得出的辨识结果：');a1=siit1_xiashi(1)a2=siit1_xiashi(2)b1=siit1_xiashi(3)b2=siit1_xiashi(4)。
系统辨识第5章 线性动态模型参数辨识 最小二乘法
度函数
，则称uS(uk()为) “持续激励”信号。
● 定义4 一个具有谱密度 Fn (为z 1的) 平f1z稳1 信f2号z 2u(k)称fn为z nn 阶
“持续激励”Fn信(e号j )，2 S若u (对) 一0 切形如 Fn (e j ) 0
的滤波器，关系式
，意味着
。
● 定理2 设输入信号u(kR)u是(0)平稳R随u (1机) 信号，Ru (如n 果1)相关函数矩阵
式中
zL H L nL
nzHLLL[[zn(h(hh11TT)T),((,(zL12n())()22)),,,,znz(((LzLzL)(()]10]))1)
z(1 na ) z(2 na )
z(L na )
u(0) u(1)
u(L 1)
u(1 nb )
u(2
nb
)
u(L nb )
5.2 最小二乘法的基本概念
● 两种算法形式
① 批处理算法：利用一批观测数据，一次计算或经反复迭代，
以获得模型参数的估计值。
②
递推算法：在上次模型参数估计值
ˆ
(k
1)的基础上，根据当
前获得的数据提出修正，进而获得本次模型参数估计值ˆ (k )，
广泛采用的递推算法形式为
(k ) (k 1) K (k )h(k d )~z (k )
z(k ) h (k ) n(k )
式中z(k)为模型输出变量，h(k)为输入数据向量， 为模型参
数向量，n(k)为零均值随机噪声。为了求此模型的参数估计值， 可以利用上述最小二乘原理。根据观测到的已知数据序列
和{z(k)} ，{h极(k小)} 化下列准则函数
L
J ( ) [z(k ) h (k ) ]2
系统辨识之最小二乘
设
n
x (k) = n(k) + å ain(k - i) i =1
（4）
则式（3）可以写成
y(k) = -a1y(k -1) -!- an y(k - n) + b0u(k) + b1u(k -1) +!+
bnu(k - n) + x (k)
(5)
在测量 u(k) 时也有测量误差，系统内部也可能有噪声，应当考虑它们的影响。因此假定
ê ê
y(n
+
2)
ú ú
=
ê ê
- y(n + 1)
ê!úê
"
! - y(2)
u(n + 2)
u(2)
ú ú
´
"
"
"ú
ê ë
y(n
+
N
)úû
êë- y(n + N -1) ! - y(N ) u(n + N ) u(N )úû
éa1 ù
ê ê
!
êêêbb0n
ê! ê
ú ú ú ú ú ú ú
+
é x (n + 1) ù
个方程，即
y(n +1) = -a1 y(n) - a2 y(n -1) -!- an y(1) + b0u(n +1) + b1u(n) + !+ bnu(1) + x (n +1) y(n + 2) = -a1 y(n +1) - a2 y(n) -!- an y(2) +
b0u(n + 2) + b1u(n +1) +!+ bnu(2) + x (n + 2)
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
最小二乘法参数辨识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))(1-zG )()(11--=z A z B )(1-zN )()(11--=z C z D 1121211212()1()a a b b n n n n A z a z a z a z B z b z b z b z --------⎧=++++⎪⎨=+++⎪⎩⎪⎩⎪⎨⎧+++=++++=--------b b a a n n n n z d z d z d z D z c z c z c z C 2211122111)(1)(2.1.3 准则函数设一个随机序列{}),,2,1(),(L k k z ∈的均值是参数θ的线性函数： {}θ)()(k h k z E T =，其中)(k h 是可测的数据向量，那么利用随机序列的一个实现，使准则函数：21])()([)(∑=-=Lk T k h k z J θθ (式2-2)达到极小的参数估计值∧θ称作θ的最小二乘估计。
最小二乘格式：)()()(k e k h k z t+=θ，θ为模型参数向量，()k e 为零均值随机2.2 广义最小二乘法2.2.1 广义最小二乘数学模型)()(1)()()()(111k v z C k u z B k z z A ---+= 式中，u(k)和)(k z 表示系统的输入输出；v(k)是均值为零的不相关的随机序列；且⎪⎩⎪⎨⎧++++=+++=++++=------------c c b b a a n n n n n n z c z c z c z C z b z b z b z B z a z a z a z A 2211122111221111)()(1)(2.2.2 广义最小二乘递推算法如下[][]⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧--=+--=--+-=--=+--=--+-=--)1()]()([)(1)()1()()()1()()]1(ˆ)()(ˆ)[()1(ˆ)(ˆ)1()]()([)(1)()1()()()1()()]1(ˆ)()()[()1(ˆ)(ˆ11k k k k k k k k k k k k k ek k k k k k k k k k k k k k k k z k k k ee e e ee e e e e ee e e e ff f f ff f f f f fffP h K P h P h h P K h K P h K P h P h h P K h K ττττττθθθθθθI I式中⎪⎩⎪⎨⎧-=----=------=)(ˆ)()()(ˆ)](ˆ,),1(ˆ[)()](,),1(),(,),1([)(k k k z k en k e k e k n k u k u n k z k z k cebffa fffθτττh h h2.2.3 广义最小二乘递推算法的计算步骤：1.给定初始条件 ⎪⎪⎩⎪⎪⎨⎧====I )0(P )0(ˆ)(I )0(P )()0(ˆe e 2f 0θθ为充分大的数充分小的实向量a a ε2利用式)()(z )()()(z )(11k z C k u k z C k z ff--==，计算)(k z f 和)(k u f；3利用式⎩⎨⎧------==ττ)](,),1(),(,),1([)(],,,,,[f f f f f 11ba n n n k u k u n k z k z kb b a a h θ，构造)(fk h ；4利用式[]⎪⎪⎩⎪⎪⎨⎧--=+--=--+-=-)1()]()([)(1)()1()()()1()()]1(ˆ)()()[()1(ˆ)(ˆ1k k k k k k k k k k k k k z k k k ff f f ff f f f f fffP h K P h P h h P K h K τττθθθI 递推计算)(ˆk θ； 5利用)(ˆ)()()(ˆk k k z k eθτh -=和 τ)](,),1(),(,),1([)(ban k u k u n k z k z k ------= h 计算)(ˆk e ；6根据τ)](ˆ,),1(ˆ[)(c e n k e k ek ----= h 来构造)(k e h ；7利用[]⎪⎩⎪⎨⎧--=+--=--+-=-)1()]()([)(1)()1()()()1()()]1(ˆ)()(ˆ)[()1(ˆ)(ˆ1k k k k k k k k k k k k k e k k k ee e e ee e e e e eeeeeP h K P h P h h P K h K τττθθθI返回第2步进行迭代计算，直至获得满意的辨识结果。