应用最小二乘一次完成法和递推最小二乘法算法的系统辨识
目录
1引言 (2)
1.1概述 (2)
1.2辨识的基本步骤 (2)
2系统辨识输入信号的产生方法和理论依据 (3)
2.1白噪声序列 (3)
2.1.1白噪声序列的产生方法 (3)
2.2 M序列的产生 (4)
2.2..1 伪随机噪声 (4)
2.2.2 M序列的产生方法 (4)
3应用经典辨识方法的辨识方案。 (6)
3.1经典辨识方法概述 (6)
3.2经典辨识方法的实现 (6)
4最小二乘法的理论基础 (7)
4.1最小二乘法 (7)
4.1.1最小二乘法估计中的输入信号 (9)
4.1.2最小二乘估计的概率性质 (9)
4.2递推最小二乘法 (10)
5两种算法的实现方案 (11)
5.1最小二乘法一次完成算法实现 (11)
5.1.1最小二乘一次完成算法程序框图 (11)
5.1.2一次完成法程序 (11)
5.1.3一次完成算法程序运行结果 (11)
5.1.4辨识数据比较 (12)
5.1.5程序运行曲线 (12)
5.2递推最小二乘法的实现 (12)
5.2.1递推算法实现步骤 (12)
5.2.2程序编制思路： (13)
5.2.3递推最小二乘法程序框图 (14)
5.2.4程序运行曲线 (15)
5.2.5测试结果 (16)
5.2.6递推数据表 (16)
6结论 (16)
7参考文献 (17)
8附录 (17)
1
应用最小二乘一次完成法和递推最小二乘法算法的系统辨识
摘要：本题针对一个单输入单输出系统的便是问题，辨识的输入信号采用的是伪随机二位式序列（M序列），系统噪声为独立同分布高斯随机向量序列（白噪声），辨识的算法是递推最小二乘法和广义最小二乘法，本文简单描述应用经典辨识方法的辨识方案，详细描述了输入信号、噪声的产生方法及matlab程序，阐述了用两种不同算法的辨识原理并对它们的推导过程及辨识程序编制思路做了详细的描述。最后结合真值与估计值对不同辨识算法的优劣进行了比较。
关键词：系统辨识M序列最小二乘法
1引言
1.1概述
系统辨识是现代控制理论中的一个分支，它是根据系统的输入输出时间函数来确定描述系统行为的数学模型。通过辨识建立数学模型的目的是估计表征系统行为的重要参数，建立一个能模仿真实系统行为的模型，用当前可测量的系统的输入和输出预测系统输出的未来演变，以及设计控制器。
对系统进行分析的主要问题是根据输入时间函数和系统的特性来确定输出信号。对系统进行控制的主要问题是根据系统的特性设计控制输入，使输出满足预先规定的要求。而系统辨识所研究的问题恰好是这些问题的逆问题。通常，预先给定一个模型类μ＝｛M｝（即给定一类已知结构的模型），一类输入信号u和等价准则J ＝L(y，yM)(一般情况下，J是误差函数，是过程输出y和模型输出yM的一个泛函)；然后选择使误差函数J 达到最小的模型，作为辨识所要求的结果。系统辨识包括两个方面：结构辨识和参数估计。在实际的辨识过程中，随着使用的方法不同，结构辨识和参数估计这两个方面并不是截然分开的，而是可以交织在一起进行的。
1.2辨识的基本步骤
①先验知识和建模目的的依据。先验知识指关于系统运动规律、数据以及其他方面的已有知识。这些知识对选择模型结构、设计实验和决定辨识方法等都有重要作用。用于不同目的的模型可能会有很大差别。
②实验设计。辨识是从实验数据中提取有关系统信息的过程，设计实验的目标之一是要使所得到的数据能包含系统更多的信息。主要包括输入信号设计，采样区间设计，预采样滤波器设计等。
③结构辨识。即选择模型类中的数学模型M的具体表达形式。除线性系统的结构可通过输入输出数据进行辨识外，一般的模型结构主要通过先验知识获得。
④参数估计。知道模型的结构后，用输入输出数据确定模型中的未知参数。实际测量都是有误差的，所以参数估计以统计方法为主。
⑤模型适用性检验。造成模型不适用主要有三方面原因：模型结构选择不当；实验数据误差过大或数据代表性太差；辨识算法存在问题。检验方法主要有利用先验知识检验和利用数据检验两类。[1]
1.3设待辨识系统如图1所示。
2
3
二阶系统为：
1121211212112
12()1()()1a z a z a z b z b z b z f z f z f z ---------=++=+=++
参数真值为：12120121.5;0.7;1;0.50;1;1;0.2a a b b f f f =-=====-=-
1．设1()1,()()f z k k ξε-==，()k ξ为白色噪声，简单描述应用经典辨识方法的辨识方案。
2．()k ξ为有色噪声，11212()1f z f z f z ---=++，121.00;0.2f f =-=，()k ε为独立同分布的高斯序列，
0.4σ=，
2系统辨识输入信号的产生方法和理论依据
2.1白噪声序列
如果随机序列 均值为0，并且两两不相的关的，对应自相关函数为
式中 则称这种随机序列为白噪声序列 2.1.1白噪声序列的产生方法
下面主要介绍（0，1）均匀分布和正态分布随机数的产生方法
在计算机上产生（0，1）均匀分布随机数的方法很多，其中最简单、最方便的是数学方法。产生伪随机数的数学方法很多，其中最常用的是乘同余法和混合同余法。 （1）乘同余法
这种方法先用递推同余式产生正整数序列{Xi=Axi-1(modM),i=1,2,3…式中：M 为2的方幂，k 为大于2的整数；A ≡3(mod8)或A ≡5(mod8),且A 不能太小；初值x0取正奇数，例如取x0=1.
再令 则 是伪随机序列，循环周期可达
。 （2）混合同余法
混合同余法产生伪随机数的递推同余式为 式中 K 为大于2的整数；A ≡1(mod4),即21n
A =+，其中n 为满足关系式2≦n ≦34
的整数。初
{()}w t ()2,0,1,2,w l R l l σδ==±±
{
1,0
0,0
l l l δ=≠=
,1,2,...
i
i x
i M
ξ=={}i
ξ22k -1
(mod )i i x Ax c M -=+2,k M =
4
值x0为非负整数。令,i
i x M
ξ=
则 ()ξi 是循环周期为 2k 的伪随机数序列 2.2 M 序列的产生
在进行系统辨识时，选用白噪声作为辨识输入信号可以保证获得较好的便是效果，但工程上难以实现。M 序列是一种很好的辨识输入信号，它具有近似白燥声的性质，不仅可以保证有较好的辨识效果，而且工程上易于实现。
M 序列是伪随机二位式序列的一种形式。在介绍M 序列之前，先介绍一下伪随机噪声的概念。 2.2..1 伪随机噪声
对白噪声的一个样本函数w(t)截取[0,T]时间内一段，对其它时间段[T,2T],[2T,3T],…,以 周期 T 延拖下去，这样获得的函数w(t)是周期T 的函数，在[0,T]时间内是白噪声，在此时间之外是重复的白噪声，它的自相关函数 的周期也是T 。由于在[0,T]时间内自相关函数 就是白噪声的自相关函数，它具有周期性，称为w(t)为伪随机噪声。 2.2.2 M 序列的产生方法
M 序列是一种离散二位式：随机序列，所谓“二位式”是指每个随机变量只有2种状态。可用多级线性反馈移位寄存器产生M 序列。M 序列是最长线性反馈移存器序列的简称，是由带线性反馈的移存器产生的周期最长的一种序列。具有较强的抗干扰能力和较低的截获概率，而且长的M 序列更容易在一定的强噪声中被提取，这样就能够充分保证数据的正常通信。
通常产生伪随机序列的电路为反馈移存器.一般说来,
由n 级移位寄存器产生的周期为N=2?-1的M 序列，在
一个循环周期内，“0”出现的次数为1/2N -,“1”出现的次数为1/2N +。现在我们引入M 序列的本原多项式的概念。若一个n 次多项式()f x 满足以下条件 （1）()f x 为既约的。
（2）()f x 可整除(1)xm +,21m n =+。
（3）()f x 除不尽(1)xq +,q m <则()f x 为本原多项式。
一个4级M 序列可以通过线性反馈移位寄存器产生，如下图所示：
图2 周期为15的伪随机序列产生器图
()()()w R E w t w t ττ=+????()w R τ
5
每级移位寄存器由双稳态触发器和门电路组成，称为1位，分别以0和1来表示2种状态。当移位脉冲到来时，每位的内容移至下一位，最后1位移出的内容即为输出。为了保持连续工作，将最后2级寄存器的内容经过适当的逻辑运算后反馈到第1级寄存器作为输入。当一个移位脉冲到来后，第1级寄存器的内容送到第2级，第2级寄存器的内容送到第3级，第3级寄存器的内容送到第4级，而第3级和第4级寄存器的内容作模和2相加后再反馈到第1级寄存器。产生伪随机序列时要求寄存器的初始状态不全为0，因为全0初始状态将导致各级寄存器输出永远是0。如果寄存器的初始内容都是1，第1个移位脉冲来到后，4级寄存器的内容变以0111，一个周期的变化规律为:1111 0111 0011 0001 1000 0100 0010 1001 1100 0110 1011 0101 1010 1101 1110 1111一个周期结束后，产生的15种不同的状态。
计算机模拟产生M 序列非常方便，先定义输出序列长度和一个数组，数组个数等于移位寄存器的个数，通过使用异或指令，再利用for 循环指令，即可完成任意长度和级数M 序列的产生。
图3M 序列产生的流程图
图4MATLAB 仿真M 序列图形
如图3是产生M 级长度为L 的M 序列程序流程图，图4所示是Matlab 仿真4级移位寄存器，长度15的M 序列输出图形。[2]
6
3应用经典辨识方法的辨识方案。
3.1经典辨识方法概述
用正弦信号、单位阶跃信号和M 序列辨识系统的方法都称为经典辨识方法。在经典辨识方法中，用得最多的是脉冲响应。用M 序列作为输入信号，再用相关法处理测试结果，可以很方便地得到系统的脉冲响应。用白噪声作为输入信号确定系统脉冲响应的方法，称为相关法。相关法，就是根据维纳－霍夫积分方程，利用输入信号的自相关函数和输入与输出的互相关函数确定系统脉冲响应的方法。采用周期性的伪随机信号作为输入信号可以使计算变得简单，所以用M 序列辨识系统脉冲响应是用得最广泛的一种方法。
3.2经典辨识方法的实现
由系统辨识图可得如下关系统：
(k)(k)y =x +(k )g(k )(k )
ξξ=+ （3.2.1）
其中：()k ξ为输入白噪声信号，(k)x 为实际输入信号，(k)y 为实际输出信号，g(k )为输入信号为u(k )时的
脉冲响应。
根据维纳-霍夫积分方程可得离散的维纳霍夫积分方程：
()1
()()-==???-?∑N ux x k R g k R k τμ （3.2.2）
令式中=?τμ、=?μμ、=?k k 则上式可写成以下形式：
()1
()()-==?-∑N ux x k R g k R k μμ （3.2.3）
再设：
()()()()012¨¨??
??
????=??
??
????
g g g g g n (0)(1)(2)
(1)????????=??????-??ux ux ux ux ux R R R R R N (0)(1)(1)(1)(0)(2)(1)(2)(0)--+??
??-+??=????
--??
u u u u u u u u u R R R N R R R N R R N R N R
再根据上式（3.2.3）可得：
=?ux R Rg 1
1-=
?
ux g R R （3.2.4） 由二电平M 序列的自相关函数的计算公式可以求得
7
2
111111
111??-
-
??????--??=??
??
??
--????N N R a N
N N
N
()12
211121111
2-??????=??+????
N R a N (3.2.5) 根据递推公式求m 次观测的脉冲响应g(m)
对系统进行m 次观测，m μ≥由m 次观测得到的ux R ()μ用ux R (,m )μ来表示，则有如下递推式：
[]m ux ux
ux k=011
R (m,m)=y(x)x(k-m)=R (m,m-1)+y(m)x(m-m)-R (m,m-1)m+1m+1
∑ （3.2.6） 综合以上公式就可以得到以下递推公式：
()1m m-1m 221
1x(m)12
1x(m-1)1N g =g +y(m )g m+1a N+1D 112x(m-N+1)-?????
????????
????-???
???????
?????????
(3.2.7)
根据递推式，可从1m g -及新的观测数据得到m g ，随着观测数据的增加，m g 的精确度不断的提高。[1][3]
4最小二乘法的理论基础
4.1最小二乘法
设单输入单输出线性定长系统的差分方程表示为：
其中 n(k)均值为0的白噪声,现分别测出n+N 个输出输入值y(1),y(2),…,y(n+N),u(1),u(2),…,u(n+N),则可写出N 个方程，写成向量－矩阵形式
(4.1.1)
()()()()()()()()
1201121n n y k a y k a y k a y k n b u k b u k b u k n k ξ=-------++-++-+ ()()()1
n
i i k n k a n k i ξ==+-∑
()()()()
()()()()()()()()()()()(
)()()10111121222112n n y n y n y u n u y n y n y u n u y n N y n N y N u n N u N a n a n b n N b ξξξ+--+????????+-+-+????=?????????+-+--+????????
??
??+??
??????+??+???????
???+?????????
?
8
则式(4.1.1)可写为 （4.1.2） 式中：y 为N 维输出向量；ξ为N 为维噪声向量；θ为（2n+1）维参数向量；Φ为N ×（2n+1）测量矩阵。因此，式（4.1.1）是一个含有（2n+1）个未知参数，由N 个方程组成的联立方程组。
11y θφφξ--=-
在给定输出向量y 和测量矩阵Φ的条件下求参数θ的估计，这就是系统辨识问题。
设 表示 θ 的估计值，?表示y 的最优估计，则有 （4.1.3） 式中：
()()()10??1??2??,???n n a
y n a y n y b y n N b θ????
+????????+????==????????+??????
????
设e(k)=y(k)- ?(k), e(k)称为残差，则有e=y- ?=y-Φθ 最小二乘估计要求残差的平方和最小，即按照指数函数
（4.1.4）
求Ｊ对 的偏导数并令其等于0可得：
（4.1.5）
由式（4.1.5）可得的 θ 最小二乘估计：
（4.1.6）
J 为极小值的充分条件是：
即矩阵ΦT Φ为正定矩阵，或者说是非奇异的。
()()()()()()()
()()()()()()()()()()()101122,,11112221n n a y n n y n a n y b y n N n N b y n y u n u y n y u n u y n N y N u n N u N ξξθξξφ????++????????????++?
???===??????????????++??????????
????
--+????-+-+??=????-+--+????
y φθ
ξ=+?θ??y θ=Φ()()
??T
T
J e e y y θ
θ==-Φ-Φ?θ()
?20?T J y θ
θ
?=-Φ-Φ=??T T y θ
ΦΦ=Φ()1?T T y
θ-=ΦΦΦ220?T
J θ
?=ΦΦ>?
9
4.1.1最小二乘法估计中的输入信号
当矩阵ΦT Φ的逆阵存在是，式（4.1.6）才有解。一般地，如果u(k)是随机序列或伪随机二位式序列，则矩阵ΦT Φ是非奇异的，即（ΦT Φ）－1存在，式（4.1.6）有解。 现在从ΦT Φ必须正定出发，讨论对u(k)的要求。
（4.1.7）
当N 足够大时有
（4.1.8）
如果矩阵ΦT Φ正定，则Ru 是是对称矩阵，并且各阶主子式的行列式为正。当N 足够大时，矩阵Ru 才是是对称的。
由此引出矩阵ΦT Φ正定的必要条件是u(k)为持续激励信号。如果序列{u(k)}的n+1阶方阵Ru 是正定的，则称{u(k)}的n+1阶持续激励信号。 下列随机信号都能满足Ru 正定的要求 1. 有色随机信号 2. 伪随机信号 3.
白噪声序列
4.1.2最小二乘估计的概率性质 最小二乘估计的概率性质 1） 无偏性
由于输出值y 是随机的，所以θ
是随机的，但注意θ不是随机值。
如果{}
{}?E E θ
θθ==，则称?θ是θ无偏估计 2）一致性
估计误差θ
的方差矩阵为
（4.1.9）
上式表明当N →∞时，θ
以概率1趋近于θ。
当ξ(k )为不相关随机序列时，最小二乘估计具有无偏性和一致性 3)有效性
如果式（4.1.2）中的ξ的均值为0且服从正态分布的白噪声向量，则最小二乘参数估计值 为有效值，参数估计偏差的方差达到Cramer-Rao 不等式的下界，即
（4.1.10）
式中M 为Fisher 矩阵，且
1
n N yy yu T
k n
uy uu +-=ΦΦ??ΦΦ=
??ΦΦ??
∑
..1
1y uy W P T yu
u R R R R R N ??
ΦΦ???→=????
1
21?T Var E N N
σθ-????
??=ΦΦ?? ???????
2
1?lim lim 0,..1N N Var R w p N
σθ
-→∞→∞==?θ(){
}
1
21
?T Var E M θ
σ--=ΦΦ=
10
（4.1.11）
4）渐近正态性
如果式（4.1.2）中的ξ是均值为0且服从正态分布的白噪声向量，则最小二乘参数估计值?θ服从正态分布，即
（4.1.2）
4.2递推最小二乘法
为了实现实时控制，必须采用递推算法，这种辨识方法主要用于在线辨识 设已获得的观测数据长度为N ，则
估计方差矩阵为 式中 于是
如果再获得1组新的观测值，则又增加1个方程
得新的参数估计值 式中
应用矩阵求逆引理，可得1N P +和N P 的递推关系式
矩阵求逆引理 设A 为n ×n 矩阵，B 和C 为n ×m 矩阵，并且A ，A ＋BCT 和I+CTA-1B 都是非奇异矩阵，则有矩阵恒等式
（4.2.2）
得到递推关系式
由于1T
N N N P ψψ+是标量，因而上式可以写成
（4.2.3）
最后，得最小二乘法辨识公式
()
(
)??ln /ln /??T
p y p y M E θθ
θθ
?????????????
?=?
??
????????
?
?
???
()
{
}(
)
12?,T N E θθσ-ΦΦ N N N
Y θξ=Φ+()1
?N
T T N N N N
Y θ-=ΦΦΦ()1
22T N N N
Var P θσσ-=ΦΦ= ()
1
T N N N P -=ΦΦ?N
T N N N
P Y θ=Φ111
T
N N N y ψθξ+++=+()1111
?T N N N N N N P Y y θψ++++=Φ+()1
1111111T
N N T
T T N N N N N N P P ψψψψ---+++++??ΦΦ?????
?==+????
???
???????
()()1
1
1111
T T T A BC A A B I C A B C A ------+=-+()1
1111T
T
N N N N N
N N N N
P P P I P P ψψψψ-++++=-+()1
11111T T
N N N N N N N N N
P P P P P ψψψψ-++++=-+()
1111???T N N N N N N
K y θθψθ++++=+-()
1
1111T N N N N N N K P P ψψψ-+++=+()1
111111T T N N N N N N N N N
P P P P P ψψψψ-+++++=-
+（4.2.1）
11 （4.2.4）
有2种给出初值的办法
（1）设N0（N0>n ）为N 的初始值，可算出
（4.2.5）
（2）假定200
?0,,P c I θ==C 是充分大的常数，I 为(2n+1) ×(2n+1)单位矩阵，则经过若干次递推之后能得到较好的参数估计。[1][2][4]
5两种算法的实现方案
5.1最小二乘法一次完成算法实现
如果把式（4.1.6）中的?θ取好足够的输入—输出数据以后一次计算出来，那么这种算法就式最小二乘法的一次完成法。
5.1.1最小二乘一次完成算法程序框图
5.1.2一次完成法程序 具体程序参见附录4
5.1.3一次完成算法程序运行结果 ans =
-1.5000 0.7000 1.0000 0.5000 a1 = -1.5000 a2 =
()10000000
?,T T N N N N N N N P P Y θ-=ΦΦ=Φ
12
0.7000 b1 =
1.0000
5.1.4辨识数据比较
5.1.5程序运行曲线
-10
1
输入u (k )
0102030405060
-10
010
输出z (k
)
0102030405060
-10
010离散输出z (k )
5.2递推最小二乘法的实现
5.2.1递推算法实现步骤
1）输入M 序列的产生程序，可以参见3.2当中M 序列产生方法（具体程序见附录。） 2）初始化。
一种简单的方法是选择?0n
θ=和2n P C I =，其中C 为尽量大的数，然后以它们为起始值进行计算(参照式4.2.3)。可以证明，当C →∞时，按照递推公式算得的将与上面用批处理算法求得的结果相同，当N 很大时，C 的选择便无关紧要了。 3）递推算法的停机标准：
()()()1??1max
?i i
i
i
k k k θθεθ
-?--<，ε是适当小的数。
4）最小二乘估计的递推算法系统参数的步骤如下：
13
（1）产生系统输入信号M 序列，输入系统阶次，计算输入和输出数据u （i ），y （i ），i=1,2,…,n;
（2）置N=n ，10?0,10N N
P I θ==（I …单位矩阵）； （3）形成行向量[]1()(1)()(1N y n y N n u N u N n ψ+=--+-+-
（4）计算()()()1
11111;T
T
N N N K P N N P N ψ
ψψ-+??=++++??
（5）输入新测量数据y ( N+1 )和u （N+1）；
（6）计算()11???(1)1;N N N N
K y N N θθψθ++??=++-+?
?
（7）计算[]11
1;N N N N P I K P ψ+++=- （8）输出1
?N θ+和1N P +； （9）N ←N+1,转（3）； 5.2.2程序编制思路：
（1）产生M 序列，通过计算，依据算出输出值，设定递推初始值，采用4.2.4方法2设定辨识参数初值，0
?θ为一个很小的量，P0为一个很大值的4×4矩阵，设定误差量，作为停机标准。 （2）依据设定，计算[]
1
(2)(1)(1)(2)T
y y u u ψ=--，1011T P ψ
ψ??+可以算出为一标量，依
据式4.2.5算出K1,K1为4×1矩阵。
（3）依据式4.2.5计算出1?θ，和1
P ，加入估计参数矩阵； （4）计算误差大小，求出相对变化率，加入误差矩阵第一列。
（5）再以1?θ和1
P 为已知参数，重复（2）－（3）计算出参数2?θ，并加入估计参数矩阵。 （6）如此循环，计算出参数估计?n
θ。 （7）判断误差变化率满足要求否，如果满足，则退出循坏，不再进行计算。 （8）分离估计参数，绘出图形。
5.2.3递推最小二乘法程序框图
如下所示：
14
15
具体程序参见附录 5.2.4程序运行曲线
-0.5
-0.4-0.3-0.2-0.100.10.20.30.4
0.5
图6M 序列信号的波形
0102030405060708090100
-1.5
-1
-0.5
0.5
1
1.5
Parameter Identification with Recursive Least Squares Method
16
0102030405060708090100
-500
-400-300-200-1000100200
300Identification Precision
5.2.5测试结果 见附录6
5.2.6地退数据表
差的分布趋势。（误差=
-新估计值旧估计值旧估计值
）[5]
6结论
最小二乘法是应用广泛的辨识方法之一。可用于动态系统，也可用于静态系统；可用于线性系统，也可用于非线性系统。通过本课程设计，了解和掌握了基于最小二乘法原理的一次性完成算法和递推算法的实现方
法，完成了Matlab下的仿真计算，能够精确的估计辨识参数。最小二乘一次性完成算法是离线算法，需要采集大量数据，一次性完成计算，因此，数据计算量大，当数据量很大时，数据输入不方便，但在本课程设计过程当中，考虑到了此问题，运用相应的办法，解决了矩阵输入的问题。递推算法适合于在线算法，利用原有参数估计进行下一步估计，可以做到运算量小，实时进行估计，根据仿真结果图示，可以看出计算一定量的数据后，参数估计趋于稳定，只要满足误差要求，不必计算完所有数据。两种算法不足之处，在考虑噪声干扰时，看成了白噪声，具有不相关性。如果噪声引起的方程误差是相关的，而不是白噪声，可以通过下面两种方法进行估计①先估计未受噪音污染的系统输出，然后再估计模型的参数，就是辅助变量法；②使用一种滤波器，将相关的方程的误差转换为白噪声再进行估计，就是增广矩阵法。
7参考文献
[1]李言俊张科,系统辨识理论及应用,国防工业出版社,2006.7 .
[2]刘宏才,系统辨识与参数估计,北京,冶金工业出版社,1999.
[3]黄道平 ,MATLAB与控制系统的数字仿真及CAD化学工业出版社, 2004.10
[4]侯媛彬汪梅王立奇系统辨识及其MATLAB仿真科学出版社 2004.2
[5]蔡季冰，系统辨识，北京，北京理工大学出版社，1991.
8附录
附录1：用乘同余法产生白噪声
①编程如下：
A=6; x0=1; M=255; f=2; N=100；%初始化；
x0=1; M=255;
for k=1: N %乘同余法递推100次；
x2=A*x0; %分别用x2和x0表示x i+1和x i-1；
x1=mod (x2,M); %取x2存储器的数除以M的余数放x1（x i）中；
v1=x1/256; %将x1存储器中的数除以256得到小于1的随机数放v1中；
）减去0.5再乘以存储器f中的系数，存放在矩阵存储器v的第v(:,k)=(v1-0.5 )*f; %将v1中的数（
i
k列中，v(:,k)表示行不变、列随递推循环次数变化；
x0=x1; % x i-1= x i；
v0=v1;
end %递推100次结束；
v2=v %该语句后无‘；’，实现矩阵存储器v中随机数放在v2中，且可直接显示在
MATLAB的window中；
k1=k;
%grapher %以下是绘图程序；
k=1:k1;
plot(k,v,k,v,?r?);
xlabel(…k?), ylabel(…v?);tktle(… (-1,+1)均匀分布的白噪声?)
1.2程序运行结果如图6所示。
17
18
图6 采用MA TLAB 产生的(-1,+1)均匀分布的白噪声序列
附录2 产生的(-1，1)均匀分布的白噪声序列
在程序运行结束后，产生的(-1，1)均匀分布的白噪声序列，直接从MA TLAB 的window 界面中copy 出来如下（v2中每行存6个随机数）：
v2 =
-0.9531 -0.7188 0.6875 -0.8359 -0.0156 0.9219
0.5703 0.4531 -0.2500 -0.4844 0.1016 -0.3672 0.8047 -0.1328 0.2188 0.3359 -0.9531 -0.7188 0.6875 -0.8359 -0.0156 0.9219 0.5703 0.4531 -0.2500 -0.4844 0.1016 -0.3672 0.8047 -0.1328
0.2188 0.3359 -0.9531 -0.7188 0.6875 -0.8359 -0.0156 0.9219 0.5703 0.4531 -0.2500 -0.4844
0.1016 -0.3672 0.8047 -0.1328 0.2188 0.3359 -0.9531 -0.7188 0.6875 -0.8359 -0.0156 0.9219
0.5703 0.4531 -0.2500 -0.4844 0.1016 -0.3672 0.8047 -0.1328 0.2188 0.3359 -0.9531 -0.7188 0.6875 -0.8359 -0.0156 0.9219 0.5703 0.4531 -0.2500 -0.4844 0.1016 -0.3672 0.8047 -0.1328
0.2188 0.3359 -0.9531 -0.7188 0.6875 -0.8359 -0.0156 0.9219 0.5703 0.4531 -0.2500 -0.4844
0.1016 -0.3672 0.8047 -0.1328 0.2188 0.3359 -0.9531 -0.7188 0.6875 -0.8359
附录3 M 序列产生程序
function [Out]=Mgenerator(length,plength,a)
%plength=4;length=15; %M 序列周期=2^plength-1 ,%置M 序列总长度和级数 ，a 是幅度 for n=1:plength %移位寄存器输入X(i)初T 态（1111） X(n)=1; end
for i=1:length for j=plength:-1:2
X(j)=X(j-1);
end
X(1)=xor(X(plength-1),X(plength)); %异或运算
if X(plength)==0
Out(i)=-1;
else
Out(i)=X(plength);
end
end
%save ('mdata','Out');
%绘图
i1=i
k=1:1:length;
plot(k,Out,k,Out,'rx')
xlabel('k')
ylabel('M序列')
title('移位寄存器产生的M序列')
end
附录4程序运行结果如下图所示
图7移位寄存器产生的M序列附录4 最小二乘一次完成法
19
%Oncefinish
%model is y(k)-1.5y(k-1)+0.7(k-2)=u(k-1)+0.5(k-2)+v(k);
L=50;A=1;
n=2; %output number;
u=mseries(A,L);
%u=Mserise(L,M,A); %系统辨识的输入信号为一个周期的M序列,信号长度L；幅度A；
z=zeros(1,16); %定义输出观测值的长度1*16
dd=zeros(1,4);
gg=zeros(14,1);
for k=n+1:L+1
z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); %用理想输出值作为观测值
end
subplot(3,1,1) %画三行一列图形窗口中的第一个图形
stem(u) %画出输入信号u的经线图形
ylabel('输入u(k)')
subplot(3,1,2)
i=1:1:L+1;
plot(i,z)
ylabel('输出z(k)')
subplot(3,1,3) %画三行一列图形窗口中的第三个图形
stem(z),grid on %画出输出观测值z的经线图形，并显示坐标网格
ylabel('离散输出z(k)')
%u,z; %显示输入信号和输出观测信号
%L=L-1%数据长度
for i=n+1:L+1
h=[-z(i-1) -z(i-2) u(i-1) u(i-2)];
dd(i-2,:)=h; %dd=(L-1)*2n
end
%HL=[-z(2) -z(1) u(2) u(1);-z(3) -z(2) u(3) u(2);-z(4) -z(3) u(4) u(3);-z(5) -z(4) u(5) u(4);-z(6) -z(5) u(6) u(5);-z(7) -z(6) u(7) u(6);-z(8) -z(7) u(8) u(7);-z(9) -z(8) u(9) u(8);-z(10) -z(9) u(10) u(9);-z(11) -z(10) u(11) u(10);-z(12) -z(11) u(12) u(11);-z(13) -z(12) u(13) u(12);-z(14) -z(13) u(14) u(13);-z(15) -z(14) u(15) u(14)] %给样本矩阵HL赋值
for j=n+1:L+1
zz=[z(j)];
gg(j-2,:)=zz;
end
%c=inv(dd'*dd)*dd'*gg;
20
递推最小二乘法推导(RLS)——全网最简单易懂的推导过程
递推最小二乘法推导（RLS）——全网最简单易懂的推导过程 作者：阿Q在江湖 先从一般最小二乘法开始说起 已知x和y的一系列数据，求解参数theta的估计。用矩阵的形式来表达更方便一些： 其中k代表有k组观测到的数据， 表示第i组数据的输入观测量，yi表示第i组数据的输出观测量。令： ，则最小二乘的解很简单， 等价于即参数解为：如果数据是在线的不断的过来，不停的采用最小二乘的解法来解是相当消耗资源与内存的，所
以要有一种递推的形式来保证对的在线更新。 进一步推导出递推最小二乘法（RLS） 我们的目的是从一般最小二乘法的解 推导出 的递推形式。一定要理解这里的下标k代表的意思，是说在有k组数据情况下的预测，所以k比k-1多了一组数据，所以可以用这多来的一组数据来对原本的估计进行修正，这是一个很直观的理解。下面是推导过程： 先看一般最小二乘法的解 下面分别对 和 这两部分进行推导变换，令
得到下面公式（1） 下面来变换得到公式（2） 下面再来，根据一般最小二乘法的解，我们知道下式成立，得到公式（3）（注：后续公式推导用到） 好了，有了上面最主要的三步推导，下面就简单了,将上面推导的结果依次代入公式即可：
至此，终于变成 的形式了。 通过以上推导，我们来总结一下上面RLS方程： 注：以上公式7中，左边其实是根据公式1，右边I为单位矩阵
公式（5）和（7）中，有些文献资料是用右边的方程描述，实际上是等效的，只需稍微变换即可。例如（5）式右边表达式是将公式（1）代入计算的。为简化描述，我们下面还是只讨论左边表达式为例。 上面第7个公式要计算矩阵的逆，求逆过程还是比较复杂，需要用矩阵引逆定理进一步简化。 矩阵引逆定理： 最终RLS的方程解为：
递推最小二乘法算法
题目： （递推最小二乘法） 考虑如下系统： )()4(5.0)3()2(7.0)1(5.1)(k k u k u k y k y k y ξ+-+-=-+-- 式中，)(k ξ为方差为0.1的白噪声。 取初值I P 610)0(=、00=∧ ）（θ。选择方差为1的白噪声作为输入信号)(k u ，采用PLS 法进行参数估计。 Matlab 代码如下： clear all close all L=400; %仿真长度 uk=zeros(4,1); %输入初值：uk(i)表示u(k-i) yk=zeros(2,1); %输出初值 u=randn(L,1); %输入采用白噪声序列 xi=sqrt(0.1)*randn(L,1); %方差为0.1的白噪声序列 theta=[-1.5;0.7;1.0;0.5]; %对象参数真值 thetae_1=zeros(4,1); %（）θ初值 P=10^6*eye(4); %题目要求的初值 for k=1:L phi=[-yk;uk(3:4)]; %400×4矩阵phi 第k 行对应的y(k-1),y(k-2),u(k-3), u(k-4) y(k)=phi'*theta+xi(k); %采集输出数据 %递推最小二乘法的递推公式 K=P*phi/(1+phi'*P*phi); thetae(:,k)=thetae_1+K*(y(k)-phi'*thetae_1); P=(eye(4)-K*phi')*P; %更新数据 thetae_1=thetae(:,k); for i=4:-1:2 uk(i)=uk(i-1); end uk(1)=u(k); for i=2:-1:2 yk(i)=yk(i-1);
几种最小二乘法递推算法的小结
一、 递推最小二乘法 递推最小二乘法的一般步骤： 1. 根据输入输出序列列出最小二乘法估计的观测矩阵?： ] )(u ... )1( )( ... )1([)(T b q n k k u n k y k y k ------=? 没有给出输出序列的还要先算出输出序列。 本例中， 2)]-u(k 1),-u(k 2),-1),-y(k -[-y(k )(T =k ?。 2. 给辨识参数θ和协方差阵P 赋初值。一般取0θ=0或者极小的数，取σσ,20I P =特别大，本例中取σ=100。 3. 按照下式计算增益矩阵G ： ) ()1()(1)()1()(k k P k k k P k G T ???-+-= 4. 按照下式计算要辨识的参数θ： )]1(?)()()[()1(?)(?--+-=k k k y k G k k T θ?θθ 5. 按照下式计算新的协方差阵P ： )1()()()1()(---=k P k k G k P k P T ? 6. 计算辨识参数的相对变化量，看是否满足停机准则。如满足，则不再递推；如不满足， 则从第三步开始进行下一次地推，直至满足要求为止。 停机准则：ε???<--) (?)1(?)(?max k k k i i i i 本例中由于递推次数只有三十次，故不需要停机准则。 7. 分离参数：将a 1….a na b 1….b nb 从辨识参数θ中分离出来。 8. 画出被辨识参数θ的各次递推估计值图形。 为了说明噪声对递推最小二乘法结果的影响，程序5-7-2在计算模拟观测值时不加噪 声， 辨识结果为a1 =1.6417，a2 = 0.7148，b1 = 0.3900，b2 =0.3499，与真实值a1 =1.642， a2 = 0.715， b1 = 0.3900，b2 =0.35相差无几。 程序5-7-2-1在计算模拟观测值时加入了均值为0，方差为0.1的白噪声序列，由于噪 声的影响，此时的结果为变值，但变化范围较小，现任取一组结果作为辨识结果。辨识结果为a1 =1.5371， a2 = 0.6874， b1 = 0.3756，b2 =0.3378。 程序5-7-2-2在计算模拟观测值时加入了有色噪声，有色噪声为 E(k)+1.642E(k-1)+0.715E(k-2)，E(k)是均值为0，方差为0.1的白噪声序列，由于有色噪声的影响，此时的辨识结果变动范围远比白噪声时大，任取一组结果作为辨识结果。辨识结果为a1 =1.6676， a2 = 0.7479， b1 = 0.4254，b2 =0.3965。 可以看出，基本的最小二乘法不适用于有色噪声的场合。
递推阻尼最小二乘法辨识算法公式的详细推导与说明
控制理论与控制工程 学位课程《系统辨识》考试报告 递推阻尼最小二乘法公式详细 推导 专业：控制理论与控制工程 班级：2011双控（研） 学生姓名：江南 学号：20110201016 任课教师：蔡启仲老师 2012年06月29 日
摘要 在参数辨识中，递推最小二乘法是用得最多的一种算法。但是，最小二乘法存在一些缺点，如随着协方差矩阵的减小，易产生参数爆发现象；参数向量和协方差矩阵的处置选择不当会使得辨识过程在参数收敛之前结束；在存在随机噪声的情况下，参数易产生漂移，出现不稳定等。为了防止参数爆发现象，Levenberg 提出在参数优化算法中增加一个阻尼项，以增加算法的稳定性。本文在一般的最小二乘法中增加了阻尼因子，构成了阻尼最小二乘法。又根据实时控制的要求，详细推到了递推阻尼最小二乘公式，实现在线辨识。 关键字：系统辨识，最小二乘法，递推算法 正文 1.题目的基本要求 已知单入单出系统的差分方程以及噪声，在应用最小二乘法进行辨识的时候，在性能指标中加入阻尼因子，详细推导阻尼最小二乘法的递推公式。 2.输入辨识信号和系统噪声的产生方法和理论依据 2.1系统辩识信号输入选择准则 （1）输入信号的功率或副度不宜过大，以免使系统工作在非线性区，但也不应过小，以致信噪比太小，直接影响辩识精度； （2）输入信号对系统的“净扰动”要小，即应使正负向扰动机会几乎均等； （3）工程上要便于实现，成本低。 2.2白噪声及其产生方法 （1） 白噪声过程 （2）白噪声是一种均值为0、谱密度为非0常数的平稳随机过程。 （3）白噪声过程定义：如果随机过程 () t ω的均值为0，自相关函数为 ()()2 R t t ωσδ= （2.2.1） 式中()t δ 为狄拉克（Dirac) 分布函数，即 (){ (),00,0 1t t t dt δδ∞ ∞=≠∞ ==? -且t （2.2.2） 则称该随机过程为白燥声过程。 2.3白噪声序列 (1) 定义 如果随机序列{() }w t 均值为0，并且是两两不相关的，对应的自相关函数为 ()2 ,0,1,2w l R l l σδ==±± 式中{1,0 0,0 l l l δ=≠=则称这种随机序列{()}w t 为白噪声序列。 2.4白噪声序列的产生方法 (1) （0，1）均匀分布随机数的产生 在计算机上产生（0，1）均匀分布随机数的方法很多，其中最简单、最方便的是数学方法。产生伪随机数的数学方法很多，其中最常用的是乘同余法和混合同余法。 ①乘同余法。
递推最小二乘法的应用
1递推最小二乘法在电厂模型辨识中的应用 电厂中大多数热工对象可以用一阶或二阶有迟延和非迟延的模型来表示，对这些模型中参数的辨识，递推最小二乘法是一种较好的方法。本文以火电厂部分典型一阶模型为例子，借助于某电厂现场数据，分别对以下几种环节进行辨识。 1.1 一阶惯性环节 火电厂中，来自锅炉的过热蒸汽，经高压调节汽门和导汽管道进入高压缸膨胀做功，高压缸的排汽回到锅炉再热器被重新加热，加热后的蒸汽经中压调节汽门进入中低压缸进一步膨胀做功，做功后的乏汽最终排入凝汽器变成凝结水，一般中压调节汽门的开度是高压调节汽门的3倍，即在机组负荷大于额定的30%或者滑压运行时，汽轮机的中压调门是完全开启的。因此，在简化模型中，汽机侧调速器一级压力与机组有功功率可以简化为一阶惯性环节如下： 1 11()1 K G s T s = + 将以上环节离散化，并写成差分方程的形式 11111111 ()[(1)](1)(1)()/,/y k a y k b u k v k a T T T b K T T =--+-+-=-= 其中 u 为调速器一级压力，y 为机组有功功率，()v k 为零均值方差为1的高斯白噪 声。 该论文依据递推最小二乘法原理，借助 MATLAB 工具编写程序，设定合适 的初始值和加权因子进行参数辨识，辨识结果为11?? 2.7547,0.9193a b ==-，由11??,a b 可得到134.12K =，112.36s T =进而得到系统的传递函数为： 134.12 ()12.361 G s s = + 下面运用递推最小二乘法对所得结果进行仿真：假设134.12K =，112.36s T =已知， 采样时间为1s T =，则计算可得
应用最小二乘一次完成法和递推最小二乘法算法的系统辨识讲解
1最小二乘法的理论基础 1.1最小二乘法 设单输入单输出线性定长系统的差分方程表示为： 其中δ(k)为服从N(0,1)的随机噪声,现分别测出n+N 个输出输入值y(1),y(2),…,y(n+N),u(1),u(2),…,u(n+N),则可写出N 个方程，写成向量－矩阵形式 (4.1.1) ()()()()()()()() 1201121n n y k a y k a y k a y k n b u k b u k b u k n k ξ=-------+ +-+ +-+()()()()()()101122,,n n a y n n y n a n y b y n N n N b ξξθξξ?? ??++????????????++? ???===??????????????++?????????? ???? ()()()()()()()()() () ()()()() ()( )()()10111121222112n n y n y n y u n u y n y n y u n u y n N y n N y N u n N u N a n a n b n N b ξξξ+--+???? ????+-+-+???? =?????????+-+--+???? ?? ???? ??+?? ??????+??+??????? ???+??????????
则式(1.1.1)可写为 （4.1.2） 式中：y 为N 维输出向量；ξ为N 为维噪声向量；θ为（2n+1）维参数向量；Φ为N ×（2n+1）测量矩阵。因此，式（4.1.1）是一个含有（2n+1）个未知参数，由N 个方程组成的联立方程组。 11y θφφξ--=- 在给定输出向量y 和测量矩阵Φ的条件下求参数θ的估计，这就是系统辨识问题。 设 表示 θ 的估计值，?表示y 的最优估计，则有 （4.1.3） 式中： ()()()10??1??2??,???n n a y n a y n y b y n N b θ???? +????????+????==????????+?????? ???? 设e(k)=y(k)- ?(k), e(k)称为残差，则有e=y- ?=y-Φθ 最小二乘估计要求残差的平方和最小，即按照指数函数 （4.1.4） 求Ｊ对 的偏导数并令其等于0可得： （4.1.5） 由式（4.1.5）可得的 θ 最小二乘估计： （4.1.6） J 为极小值的充分条件是： 即矩阵ΦT Φ为正定矩阵，或者说是非奇异的。 1.1.1最小二乘法估计中的输入信号 当矩阵ΦT Φ的逆阵存在是，式（1.1.6）才有解。一般地，如果u(k)是随机序列或伪随机二位式序列，则矩阵ΦT Φ是非奇异的，即（ΦT Φ）－1存在，式（1.1.6）有解。 现在从ΦT Φ必须正定出发，讨论对u(k)的要求。 y φθξ=+?θ??y θ=Φ()() ??T T J e e y y θ θ==-Φ-Φ?θ() ?20?T J y θ θ ?=-Φ-Φ=??T T y θ ΦΦ=Φ()1 ?T T y θ -=ΦΦΦ220?T J θ ?=ΦΦ>?1 n N yy yu T +-ΦΦ??
广义递推最小二乘辨识
广义递推最小二乘辨识 一、实验目的 1 通过实验掌握广义最小二乘辨识算法； 2 运用MATLAB编程，掌握算法实现方法。 二、实验原理 广义最小二乘法的基本思想是基于对数据先进行一次滤波预处理，然后利用普通最小二乘法对滤波后的数据进行辨识。如果滤波模型选择得合适，对数据进行了较好的白色化处理，那么直接利用普通最小二乘法就能获得无偏一致估计。 广义最小二乘法所用的滤波模型实际上就是一种动态模型，在整个迭代过程中不断靠偏差信息来调整这个滤波模型，使它逐渐逼近于一个较好的滤波模型，以便对数据进行较好的白色化处理，使模型参数估计称为无偏一致估计。理论上说，广义最小二乘法所用的动态模型经过几次迭代调整后，便可对数据进行较好的白化处理，但是，当过程的输出噪信比比较大或模型参数比较多时，这种数据白色化处理的可靠性就会下降。此时，准则函数可能出现多个局部收敛点，因而辨识结果可能使准则函数收敛于局部极小点上而不是全局极小点上。这样，最终的辨识结果往往也会是有偏的。 其收敛速度比较慢，需要经过多次迭代计算，才能得到较准确的参数估计值。一般情况下，经过多次迭代后，估计值便会收敛到稳态值。但在某些情况下（如噪声比较低时）存在局部极小值，估计值不一定收敛到准则函数的全局极小值上。为了防止参数估计值收敛到局部极小值，最好选定初值接近最优解，一般可以用最小二乘法的批处理估计值作为初值。如果系统是时变的，或为了克服数据饱和现象，可以在两次RLS算法中分别引进遗忘因子。 三、实验内容 <1> 数据获取：实验数据按照表9-1，为二阶线性离散系统的输入输出数据 <2> 数据处理：为了提高辨识精度，实验者必须对原始数据进行剔除坏数据、零均值化、工频滤波等处理。实验进行了白化滤波处理。 <3> 辨识算法：利用处理过的数据（取适当的数据长度），选择某种辨识方法（如RLS递推最小二乘法、RELS、RIV或RML等参数估计算法及F-检验或AIC定
用matlab实现最小二乘递推算法辨识系统参数
自动化专业综合设计报告 设计题目：最小二乘递推算法辨识系统参数所在实验室：自动化系统仿真实验室 指导教师： 学生姓名 班级计082-2 班 学号 撰写时间：2012-3-1 成绩评定：
一．设计目的 1、学会用Matlab实现最小二乘法辨识系统参数。 2、进一步熟悉Matlab的界面及基本操作； 3、了解并掌握Matlab中一些函数的作用与使用； 二．设计要求 最小二乘递推算法辨识系统参数，利用matlab编程实现，设初始参数为零。z(k)-1.5*z(k-1)+0.7*z(k-2)=1*u(k-1)+0.5*u(k-2)+v(k); 选择如下形式的辨识模型： z(k)+a1*z(k-1)+a2*z(k-2)=b1*u(k-1)+b2*u(k-2)+v(k); 三．实验程序 m= 3; N=100; uk=rand(1,N); for i=1:N uk(i)=uk(i)*(-1)^(i-1); end yk=zeros(1,N); for k=3:N yk(k)=1.5*yk(k-1)-0.7*yk(k-2)+uk(k-1)+0.5*uk(k-2); end %j=100;kn=0; %y=yk(m:j)'; %psi=[yk(m-1:j-1);yk(m-2:j-2);uk(m-1:j-1);uk(m-2:j-2)]'; %pn=inv(psi'*psi); %theta=(inv(psi'*psi)*psi'*y); theta=[0;0;0;0]; pn=10^6*eye(4); for t=3:N ps=([yk(t-1);yk(t-2);uk(t-1);uk(t-2)]); pn=pn-pn*ps*ps'*pn*(inv(1+ps'*pn*ps)); theta=theta+pn*ps*(yk(t)-ps'*theta); thet=theta'; a1=thet(1); a2=thet(2); b1=thet(3); b2=thet(4); a1t(t)=a1; a2t(t)=a2;b1t(t)=b1;b2t(t)=b2; end t=1:N; plot(t,a1t(t),t,a2t(t),t,b1t(t),t,b2t(t));
递推最小二乘法
线性方程组的最优求解方法 一．递推最小二乘法 设线性方程组 b Ax = (1) 则有 k b k =x :A ),(, (n k Λ,2,1=) (2) 其中，[]kn k k a a a k ,,,:),(21Λ=A ，[]T n x x x ,,,21Λ=x 。 设 x :A ),()(k k f = (3) 下面采用基于递推最小二乘法(RLS)的神经网络算法来训练权值向量x ，以获得线性方程组（1）的解x 。由式（3）可知，若以)(k f 为神经网络输出，以k b 为神经网络训练样本，以x 为神经网络权值向量，[]kn k k a a a k ,,,:),(21Λ=A 为神经网络输入向量，则解线性方程组的神经网络模型如同1所示。 图1 神经网络模型 采用RLS 算法训练神经网络权值向量x ，其算法如下： （1）神经网络输出： x :A ),()(k k f = (4) （2）误差函数：
)()(k f b k e k -= (5) （3）性能指标： ∑==n k k e J 1 2)(21 (6) （4）使min =J 的权值向量x ，即为所求的神经网络权值向量x ，这是一个多变量线性优化问题，为此，由 0=??x J 可得最小二乘递推法（RLS ）： ]),([1k k k k k k b x :A Q x x -+=+ (7) ),(),(1),(:A P :A :A P Q k k k T k T k k += (8) k k k k P :A Q I P )],([1-=+ (9) ()n k ,,2,1Λ= 随机产生初始权值向量)1,(0 n rand =x ，设n n ?∈=R I P α0（α是足够大的正数（一般取10610~10=α），n n ?∈R I 是单位矩阵），通过对样本数据训练，即可获得神经网络权值 向量x ，此即为线性方程组（1）的解。 二．具有遗忘因子的递推最小二乘估计公式为： ]),([1k k k k k k b x :A Q x x -+=+ (10) ),(),(),(:A P :A :A P Q k k k T k T k k +=λ (11) k k k k P :A Q I P )],([11-= +λ (12) 式中，1:)],(:),([)(-=k A k A k T W P ，W 为加权对角阵：
(完整word版)多种最小二乘算法分析+算法特点总结
第一部分：程序设计思路、辨识结果分析和算法特点总结 (2) 一：RLS遗忘因子法 (2) RLS遗忘因子法仿真思路和辨识结果 (2) 遗忘因子法的特点： (3) 二：RFF遗忘因子递推算法 (4) 仿真思路和辨识结果 (4) 遗忘因子递推算法的特点： (5) 三：RFM限定记忆法 (5) 仿真思路和辨识结果 (5) RFM限定记忆法的特点： (7) 四：RCLS偏差补偿最小二乘法 (7) 仿真思路和辨识结果 (7) RCLS偏差补偿最小二乘递推算法的特点： (9) 五：增广最小二乘法 (9) 仿真思路和辨识结果 (9) RELS增广最小二乘递推算法的特点： (11) 六：RGLS广义最小二乘法 (12) 仿真思路和辨识结果 (12) RGLS广义最小二乘法的特点： (14) 七：RIV辅助变量法 (14) 仿真思路和辨识结果 (14) RIV辅助变量法的特点： (16) 八：Cor-ls相关最小二乘法（二步法） (17) 仿真思路和辨识结果 (17) Cor-ls相关最小二乘法（二步法）特点： (18) 九：MLS多级最小二乘法 (19) 仿真思路和辨识结果 (19) MLS多级最小二乘法的特点： (22) 十：yule_walker辨识算法 (23) 仿真思路和辨识结果 (23) yule_walker辨识算法的特点： (24) 第二部分：matlab程序 (24) 一：RLS遗忘因子算法程序 (24) 二：RFF遗忘因子递推算法 (26) 三：RFM限定记忆法 (28) 四：RCLS偏差补偿最小二乘递推算法 (31) 五：RELS增广最小二乘的递推算法 (33) 六;RGLS 广义最小二乘的递推算法 (36) 七：Tally辅助变量最小二乘的递推算法 (39) 八：Cor-ls相关最小二乘法（二步法） (42) 九：MLS多级最小二乘法 (45) 十yule_walker辨识算法 (49)
递推最小二乘法的应用
1 递推最小二乘法在电厂模型辨识中的应用 电厂中大多数热工对象可以用一阶或二阶有迟延和非迟延的模型来表示，对这些模型中参数的辨识，递推最小二乘法是一种较好的方法。本文以火电厂部分典型一阶模型为例子，借助于某电厂现场数据，分别对以下几种环节进行辨识。 1.1 一阶惯性环节 火电厂中，来自锅炉的过热蒸汽，经高压调节汽门和导汽管道进入高压缸膨胀做功，高压缸的排汽回到锅炉再热器被重新加热，加热后的蒸汽经中压调节汽门进入中低压缸进一步膨胀做功，做功后的乏汽最终排入凝汽器变成凝结水，一般中压调节汽门的开度是高压调节汽门的3倍，即在机组负荷大于额定的30%或者滑压运行时，汽轮机的中压调门是完全开启的。因此，在简化模型中，汽机侧调速器一级压力与机组有功功率可以简化为一阶惯性环节如下： 1 11()1 K G s T s = + 将以上环节离散化，并写成差分方程的形式 11111111 ()[(1)](1)(1)()/,/y k a y k b u k v k a T T T b K T T =--+-+-=-= 其中 u 为调速器一级压力，y 为机组有功功率，()v k 为零均值方差为1的高斯白噪 声。 该论文依据递推最小二乘法原理，借助 MATLAB 工具编写程序，设定合适 的初始值和加权因子进行参数辨识，辨识结果为11?? 2.7547,0.9193a b ==-，由11??,a b 可得到134.12K =，112.36s T =进而得到系统的传递函数为： 134.12 ()12.361 G s s = + 下面运用递推最小二乘法对所得结果进行仿真：假设134.12K =，112.36s T =已知，采样时间为1s T =，则计算可得
参数估计带遗忘因子递推最小二乘法仿真(RLS)
参数估计带遗忘因子递推最小二乘法仿真（RLS ） 模型：y (N+1) = ?N+1T θ + ε(N+1) 其中 带遗忘因子的RLS 法递推算式： θ N+1 = θ N + K N+1(y (N+1) – ?N+1T θ N ) 式（2-3-5） 式（2-4-1） 式（2-4-2） 参考程序（BASIC ） 40 N=200： M=2： D=2 ’（ N —数据量；M —参数维数；D —滞后量 D ≧1） 50 DIM Y （N ），U （N ），A （M ），P （M ，M ）， C （M ，N ），X （M ），PX （M ），E （N ），A1（N ），B1（N ） ’ Y —输出；U —输入；A —参数估计；P —估计误差协方差阵；C —贮存参数估计结果；E —随机干扰；PX —工作单元 ； X —观测数据向量?；A1和B1—参数真值（一阶系统） 60 RANDOMIZE 77： ’ 伪随机数初始化 70 FOR I =1 TO M ：A （I ）=0 : P （I ，I ）=1000000.：NEXT I ’ 赋初值 111 1+++++=N N T N N N N P P K ??α?α??ρ??1112111???? ??+-=+++++N N T N N T N N N N N P P P P P )] (),...,1(),(),...,1([] ,...,,,,....,,[2121n k u k u n k y k y b b b a a a T k n n T k ------==?θ
80 B=1.：’赋遗忘因子赋值 90 FOR K= 5 TO N：’主循环 100 A1（K）= —0.9： B1（K）=1.0： C1=0. ：’赋真值110 IF K〉=50 THEN B1（K）=2.0：’参数时变120 US=15： IF RUN（0）〉0.8 THEN US= —1*US 130 U（K）= US：’给定输入 140 E（K）=RND（1）— 0.5 ’噪声 150 Y（K）= —A1（K）*Y（K—1）+ B1（K）*U（K—D）+E（K）+C1*E（K—1）：’过程仿真 160 X（1）= —Y（K—1）： X（2）= U（K—D）: ' 观测数据向量?赋值 170 GOSUB 220 ’调用RLS子程序 180 FOR I=1 TO M： C（I，K）=A（I）： NEXT I ’存入参数估计结果 200 NEXT K 210 END 220 ’********** 参数估计 RLS 子程序 **************** 225 ’** W 用于存放? N+1T θ N ； Z 用于存放? N+1 T P N ? N+1 ** 230 Z=0： W=0： FOR I=1 TO M ： PX（I）= 0： NEXT I 240 FOR I=1 TO M： FOR J=1 TO M 250 Z=Z+P（I，J）*X（I）*X（J）：PX（I）= PX（I）+ P（I，J）*X（J）： NEXT J 260 W=W+A（I）*X（I）： NEXT I： RE=Y（K）—W 270 FOR I=1 TO M： A（I）=A（I）+ PX（I）*RE/ （B+Z）： FOR J=1 TO M 280 P（I，J）=（P（I，J）—PX（I）*PX（J）/ （B+Z））/B 290 NEXT J： NEXT I 300 RETURN
带遗忘因子的递推最小二乘法
带遗忘因子的递推最小二乘法 %开环系统参数辨识，带遗忘因子的递推最小二乘估计法（FFRLS)，系统为单入单出的CAR(带控制量的自回归模型）模型,三阶系统 clear all clc a=[1-1.10.60.1];b=[10.7];d=4;%实际模型系数矩阵与纯迟延 L=1000;%仿真长度 na=length(a)-1;nb=length(b)-1;%na,nb为输出输入系数矩阵A,B的阶数 yk=zeros(na,1);%输出矩阵初始化 yk_m=zeros(na,1);%模型输出 uk=zeros(nb+d,1);%输入矩阵初始化 theta_e0=zeros(na+nb+1,1);%theta_e0为估计参数初值， a1,a2....an,b0,b1,...bn,共na+nb+1个 phi=zeros(na+nb+1,1);%phi为当前实际输出输入构成的矩阵 P=10^6*eye(na+nb+1);%修正系数初值 beta=0.99;%遗忘因子，在0.95到1之间 u=randn(L,1);%输入信号，方差为1的白噪声序列 omega=sqrt(0.1)*randn(L,1);%干扰信号，方差为0.1的白噪声序列 for i=1:L theta(:,i)=[a(2:na+1),b]';%系统实际参数值 phi=[-yk;uk(d:d+nb)];%系统输出输入矩阵 phi_e=[-yk_m;uk(d:d+nb)];%模型输出输入矩阵 y(i)=phi'*theta(:,i)+omega(i);%系统实际输出 y_m(i)=phi_e'*theta_e0;%模型输出 %递推公式 K=P*phi/(beta+phi'*P*phi); theta_e(:,i)=theta_e0+K*(y(i)-phi'*theta_e0); P=(eye(na+nb+1)-K*phi')*P/beta; %数据更新 theta_e0=theta_e(:,i); for j=na:-1:2 yk(j)=yk(j-1); yk_m(j)=yk_m(j-1); end yk(1)=y(i); yk_m(1)=y_m(i); for j=(nb+d):-1:2 uk(j)=uk(j-1); end
几种最小二乘法递推算法的小结
递推最小二乘法的一般步骤： 1. 根据输入输出序列列出最小二乘法估计的观测矩阵?： ] )(u ... )1( )( ... )1([)(T b q n k k u n k y k y k ------=? 没有给出输出序列的还要先算出输出序列。 本例中， 2)]-u(k 1),-u(k 2),-1),-y(k -[-y(k )(T =k ?。 2. 给辨识参数θ和协方差阵P 赋初值。一般取0θ=0或者极小的数，取σσ,20I P =特别 大，本例中取σ=100。 3. 按照下式计算增益矩阵G ： ) ()1()(1)()1()(k k P k k k P k G T ???-+-= 4. 按照下式计算要辨识的参数θ： )]1(?)()()[()1(?)(?--+-=k k k y k G k k T θ?θθ 5. 按照下式计算新的协方差阵P ： )1()()()1()(---=k P k k G k P k P T ? 6. 计算辨识参数的相对变化量，看是否满足停机准则。如满足，则不再递推；如不满足， 则从第三步开始进行下一次地推，直至满足要求为止。 停机准则：ε???<--) (?)1(?)(?max k k k i i i i 本例中由于递推次数只有三十次，故不需要停机准则。 7. 分离参数：将a 1….a na b 1….b nb 从辨识参数θ中分离出来。 8. 画出被辨识参数θ的各次递推估计值图形。 为了说明噪声对递推最小二乘法结果的影响，程序5-7-2在计算模拟观测值时不加噪声， 辨识结果为a1 =1.6417，a2 = 0.7148，b1 = 0.3900，b2 =0.3499，与真实值a1 =1.642， a2 = 0.715， b1 = 0.3900，b2 =0.35相差无几。 程序5-7-2-1在计算模拟观测值时加入了均值为0，方差为0.1的白噪声序列，由于噪声 的影响，此时的结果为变值，但变化范围较小，现任取一组结果作为辨识结果。辨识结果为a1 =1.5371， a2 = 0.6874， b1 = 0.3756，b2 =0.3378。 程序5-7-2-2在计算模拟观测值时加入了有色噪声，有色噪声为 E(k)+1.642E(k-1)+0.715E(k-2)，E(k)是均值为0，方差为0.1的白噪声序列，由于有色噪声的影响，此时的辨识结果变动范围远比白噪声时大，任取一组结果作为辨识结果。辨识结果为a1 =1.6676， a2 = 0.7479， b1 = 0.4254，b2 =0.3965。 可以看出，基本的最小二乘法不适用于有色噪声的场合。
递推增广最小二乘法
题目： （递推增广最小二乘法） 考虑如下系统： )2-(2.0)1()()4(5.0)3()2(7.0)1(5.1)(k k k k u k u k y k y k y ξξξ+--+-+-=-+-- 式中，)(k ξ为方差为0.1的白噪声。 取初值I P 610)0(=、00=∧ ）（θ。选择方差为1的白噪声作为输入信号)(k u ，采用RELS 法进行参数估计。 何为递推增广最小二乘法？原理？ 在CARMA 模型中 )()()()()()(111k z C d k u z B k y z A ξ---+-= 在式子中与前两种方法不同的是1)(1≠-z C ，在这里噪声)(k e 为有色噪声，即 )()1()()()()(11k c k c k k z C k e n ξξξξ++-+==- 也就是说我们要估计的参数又多了噪声中的参数。比原来的估计矩阵又多了参数故称作递推增广最小二乘法。 究其原理我认为仅仅是在递推最小二乘法的计算原理的基础上增加了一个方程来计算噪声： ∧ ∧ ∧-=θ?ξ)()()(k k y k T Matlab 代码如下： clear all close all L=1000; %仿真长度 uk=zeros(4,1); %输入初值：uk(i)表示u(k-i) yk=zeros(2,1); %输出初值 xik=zeros(2,1);%噪声的初值 u=randn(L,1); %输入采用白噪声序列 xiek=zeros(2,1);%噪声的给定初值 xi=sqrt(0.1)*randn(L,1); %方差为0.1的白噪声序列 theta=[-1.5;0.7;1.0;0.5;-1;0.2]; %对象参数真值,多了噪声的两个参数真值 thetae_1=zeros(6,1); %这回要辨识的值多了两个 P=10^6*eye(6); %题目要求的初值，单位阵阶数随着辨识参数的增加而增加
- 递推最小二乘法的应用
- 带遗忘因子的递推最小二乘法
- 递推最小二乘法
- 递推最小二乘法算法
- 递推最小二乘法的Matlab程序
- 参数估计带遗忘因子递推最小二乘法仿真(RLS)
- 递推最小二乘辨识讲解
- 最小二乘法递推最小二乘法C语言源程序
- 递推最小二乘法的应用
- 递推最小二乘法与广义最小二乘法matlab源程序
- 加权递推最小二乘法仿真
- 基于MATLAB的递推最小二乘法辨识与仿真
- 递推最小二乘法推导(RLS)——全网最简单易懂的推导过程
- 广义递推最小二乘辨识
- 几种最小二乘法递推算法的小结
- 递推最小二乘和增广矩阵法
- 递推增广最小二乘法
- 递推阻尼最小二乘法辨识算法公式的详细推导与说明
- 递推最小二乘法
- 递推最小二乘辨识