应用最小二乘一次完成法和递推最小二乘法算法的系统辨识讲解
各类最小二乘算法
β N −1 H* = N 0
β N −2
β 2( N −1) WN = 0
β 2( N −2)
0 ⋱ 1
三、递推算法 ∵
k θ(k ) = ∑ β i =1
∧
2(k −i) h (i )h T (i )
2随着采样次数的增多数据量不断增加ls估计有可能出现所谓的数据饱和现象导致递推算法不能接近参数真二关于数据饱和现象的分析所谓数据饱和现象就是随着时间的推移采集到的数据越来越多新数据所提供的信息被淹没在老数据的海洋之中
Ⅴ 各种最小二乘类参数辨识算法 §1 概 述
最小二乘类参数辨识算法（一次完成算法、递推算法） 最小二乘类参数辨识算法 （一次完成算法 、 递推算法 ） 是一种 最基本和常用的参数估计方法。但研究和应用表明， 最基本和常用的参数估计方法。 但研究和应用表明， 这一算 法仍存在明显的不足。 法仍存在明显的不足。 一、LS 算法的主要不足之处 1、当模型噪声为有色噪声时，LS 估计不再是无偏估计、一致 、当模型噪声为有色噪声时， 估计不再是无偏估计、 估计。 估计。 2、随着采样次数的增多，数据量不断增加，LS估计有可能出 、随着采样次数的增多，数据量不断增加， 估计有可能出 现所谓的“数据饱和”现象， 现所谓的“数据饱和”现象，导致递推算法不能接近参数真 值。
于是有： 于是有：
α P ( k ) P − 1 ( k − 1) = I − P ( k ) h ( k ) h T ( k )
则：
ˆ θ ( k ) = P ( k ) H * T Z * = P ( k ) α H * −1T Z * −1 + h ( k ) z ( k ) k k k k
应用最小二乘一次完成法和递推最小二乘法算法的系统辨识讲解
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 ay 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 Φ为正定矩阵，或者说是非奇异的。
最小二乘法在系统辨识中的应用(包含相关的三种算法)
但是，数据向量)(kfh中的变量均需要按照（3.2.2）式计算，然而噪声模型)(1zC并不知道。
为此需要用迭代的方法来估计)(1zC。
令)()(1z)(1kvCke（3.2.5）置,2[k],)()]3(),2(),1([)(31ccckekekekeபைடு நூலகம்h（3.2.6）就把噪声模型（3.2.5）也化成最小二乘格式)()()(kvkkeeeh由于上式的噪声已是白噪声，所以再次利用最小二乘法可获得噪声模型参数e的无偏估计。
通过极小化（1.1.4）式来计算的方法称作最小二乘法，未知模型参数最可能的值是在实际观测值与计算值之累次误差的平方和达到最小处所得到的，这种模型输出能最好地接近实际过程的输出。
2、辨识原理考虑模型（1.1.2）式的辨识问题，其中)(kz和)(kh都是可观测的数据，是待估计参数，准则函数取（1.1.4）根据（1.1.3）的定义，准则函数)(J可写成二次型的形式)()()(HzHzllllJ（1.2.1）显然上式中的Hl代表模型的输出，或者说是过程的输出预报值。
)()() 1()()()()()(1111kkkkkiikkihhPhhhhP（2.2.3）令：
] ) 1(),2 (z),1 (z[1kzkz则：
] )i()()[1()() 1(1111111ikkkkkzikkhPzHHH于是有i111)()() 1() 1(kizikkhP（2.2.4）令] )k(),2 (z),1 (z[zkz利用（2.2.3）和（2.2.4）式，可得)]1()()()[()() 1()}()() 1()]()()(){[()]()() 1() 1()[(] )i()()[()()(1111ikkkzkkkkzkkkkkkkzkkkkzikkkkkkkhhphhhPPhPPhPzHHH（2.2.5）引进增益矩阵)(kK，定义为)()()(kkkhPK（2.2.6）则（2.2.5）式写成)]1()()()[() 1()(kkkzkkkhK（2.2.7）进一步把（2.2.3）式写成11)]()() 1([)(kkkkhhPP（2.2.8）为了避免矩阵求逆运算，利用矩阵反演公式可将（2.2.8）式演变成) 1()]()([)]()() 1([)(11kkkkkkkPhKIhhPP（2.2.9）将（2.2.9）式代入（2.2.6）式，整理后有1] 1)() 1()()[() 1()(kkkkkkhPhhPK（2.2.10）综合（2.2.7）、（2.2.9）、（2.2.10）式便得到最小二乘参数估计递推算法。
系统辨识之最小二乘法
系统辨识之最小二乘法方法一、最小二乘一次性算法：首先对最小二乘法的一次性辨识算法做简要介绍如下：过程的黑箱模型如图所示：其中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 文件，实现上述算法。
最小二乘一次完成法
20世纪60年代
态系统辨识。
A str o m
..
等人将此方法用于动
一次完成法：理论研究 递推算法： 在线辨识
问题的提法及一次完成 最小二乘估计
单输入单输出系统，如图所示
即可
差分方程的形式
y (k ) ai y (k i ) biu (k i )
i 1 i 1 n n
(3.2)
其中， y(k ) y(kT ), u(k ) u(kT )
k 1, 2,
实际的测量中，数据中是有误差的。
如测量噪声、模型误差等。因此，实际 测得的输入输出关系应修正为：
1 1
~
其中 n 是噪声协方差阵， 0 为模型参数真值
一致性
● 一致性描述参数估计值的收敛情况。 ● 定理：如果模型噪声向量是零均值白 噪声，最小二乘参数估计是一致收敛的 ，即有
1 lim LS W.P. 0 L
ˆ
有效性
噪声方差估计性质
其输入输出关系为
b1 z 1 bn z n Y ( z) H ( z) U ( z) 1 a1 z 1 an z n
其中 z e
，
sT
mn ，T为采样周期，
不失一般性，不妨取 m=n。若m<n
b b b 0 m 1 m 2 n 只要取
参数估计偏差的协方性质
参数估计偏差的协方差阵是用来评价参数估计 精度的一个依据。 ● 定理：如果模型噪声向量是零均值且与数据 矩阵统计独立的随机向量，则加权最小二乘参 ˆ ~ 数估计偏差 的协方差阵可写 WLS 0 WLS 成
第五章 最小二乘法辨识
服从正态分
❖ 4）有效性
❖ 定理4：假设 (k) 是均值为零，方差为 2I 的正态
白噪声，则最小二乘参数估计量
^
是有效估计
量，即参数估计误差的协方差达到Cramer-Rao不
等式的下界
E (^
^
)(
)T
2E
(
T N
N
) 1
M 1
❖ 其中M为Fisher信息矩阵。
4、适应算法
❖ 随着更多观测数据的处理，递推最小二乘法对线性 定常系统的参数估计并非越来越精确，有时会发现
❖ 现举例说明最小二乘法的估计精度 ❖ 例5.1：设单输入－单输出系统的差分方程为
y(k) a1y(k 1) a2 y(k 2) b1u(k 1) b2u(k 2) (k)
❖ 设 u(k)是幅值为1的伪随机二位式序列，噪声 (k)是 一个方差 2可调的正态分布 N(0, 2 )随机序列。
❖ 为了克服数据饱和现象，可以用降低旧数据影响的 办法来修正算法。而对于时变系统，估计k时刻的 参数最好用k时刻附近的数据估计较准确。否则新 数据所带来的信息将被就数据所淹没。
❖ 几种算法：渐消记忆法，限定记忆法与振荡记忆法
❖ 矩阵求逆引理：设A为 n n 矩阵，B和C为 n m 矩阵，
并且A， A和 BCT I CT都A是1B 非奇异矩阵，则有矩
阵恒等式
A BCT 1 A1 A1B(I CT A1B)1CT A1
❖
令
A
PN1
，B
N 1
，C
T N 1
，根据引理有
PN1
T N 1 N 1
1
❖ 算法中，^ N 为2n+1个存贮单元（ai ,bi ,i 1,2, , n）， 而 PN 是 (2n 1) (2n 1)维矩阵，显然，将 N 换成 PN 后，存贮量大为减少（因为n为模型的阶数，一般 远远小于N）
应用最小二乘一次完成法和递推最小二乘法算法地系统辨识
1最小二乘法的理论基础1.1最小二乘法设单输入单输出线性定长系统的差分方程表示为：其中δ(k)为服从N(0,1)的随机噪声,现分别测出n+N 个输出输入值y(1),y(2),…,y(n+N),u(1),u(2),…,u(n+N),则可写出N 个方程，写成向量－矩阵形式()()()()()()()()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 ξ=-------++-++-+L L ()()()()()()()()()()()()()()()()()()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 ξξξ+--+⎡⎤⎡⎤⎢⎥⎢⎥+-+-+⎢⎥⎢⎥=⨯⎢⎥⎢⎥⎢⎥⎢⎥+-+--+⎢⎥⎢⎥⎣⎦⎣⎦⎡⎤⎢⎥+⎡⎤⎢⎥⎢⎥⎢⎥+⎢⎥+⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥+⎢⎥⎣⎦⎢⎥⎢⎥⎣⎦LLL L M M MMLL M M M(4.1.1)则式(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 ay n a y n y b y n N b θ⎡⎤⎢⎥+⎡⎤⎢⎥⎢⎥⎢⎥+⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥+⎢⎥⎣⎦⎢⎥⎢⎥⎣⎦M M M ()()()()()()()()()()()()()()()()()()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 ξξθξξφ⎡⎤⎢⎥++⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥++⎢⎥⎢⎥===⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥++⎢⎥⎢⎥⎣⎦⎣⎦⎢⎥⎢⎥⎣⎦--+⎡⎤⎢⎥-+-+⎢⎥=⎢⎥⎢⎥-+--+⎢⎥⎣⎦M M M M L LL L M MMLL y φθξ=+ˆθˆˆyθ=Φ设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 Φ为正定矩阵，或者说是非奇异的。
系统辨识方法之最小二乘法
目录一、系统辨识的定义.................................................................................................................. - 2 -二、最小二乘法的引出.............................................................................................................. - 2 -三、最小二乘法的原理.............................................................................................................. - 3 -3.1 最小二乘法一次完成推导[1]........................................................................................ - 3 -3.2最小二乘法的缺陷[ 5].................................................................................................... - 5 -四、其他系统辨识方法.............................................................................................................. - 5 -4.1 基于BP神经网络的系统辨识方法特点[3]................................................................. - 5 -4.2 基于遗传算法的系统辨识算法................................................................................... - 6 -五、结论...................................................................................................................................... - 7 -六、参考文献.............................................................................................................................. - 7 -系统辨识方法简介摘要：在研究一个控制系统过程中，建立系统的模型十分必要。
最小二乘法的原理及其应用说课讲解
最小二乘法的原理及其应用一、研究背景在科学研究中，为了揭示某些相关量之间的关系，找出其规律，往往需要做数据拟合，其常用方法一般有传统的插值法、最佳一致逼近多项式、最佳平方逼近、最小二乘拟合、三角函数逼近、帕德（Pade）逼近等，以及现代的神经网络逼近、模糊逼近、支持向量机函数逼近、小波理论等。
其中，最小二乘法是一种最基本、最重要的计算技巧与方法。
它在建模中有着广泛的应用，用这一理论解决讨论问题简明、清晰，特别在大量数据分析的研究中具有十分重要的作用和地位。
随着最小二乘理论不断的完善，其基本理论与应用已经成为一个不容忽视的研究课题。
本文着重讨论最小二乘法在化学生产以及系统识别中的应用。
二、最小二乘法的原理人们对由某一变量t或多个变量t1…..tn 构成的相关变量y感兴趣。
如弹簧的形变与所用的力相关，一个企业的盈利与其营业额，投资收益和原始资本有关。
为了得到这些变量同y之间的关系，便用不相关变量去构建y，使用如下函数模型,q个相关变量或p个附加的相关变量去拟和。
通常人们将一个可能的、对不相关变量t的构成都无困难的函数类型充作函数模型（如抛物线函数或指数函数）。
参数x是为了使所选择的函数模型同观测值y相匹配。
（如在测量弹簧形变时，必须将所用的力与弹簧的膨胀系数联系起来）。
其目标是合适地选择参数，使函数模型最好的拟合观测值。
一般情况下，观测值远多于所选择的参数。
其次的问题是怎样判断不同拟合的质量。
高斯和勒让德的方法是，假设测量误差的平均值为0。
令每一个测量误差对应一个变量并与其它测量误差不相关（随机无关）。
人们假设，在测量误差中绝对不含系统误差，它们应该是纯偶然误差，围绕真值波动。
除此之外，测量误差符合正态分布，这保证了偏差值在最后的结果y上忽略不计。
确定拟合的标准应该被重视，并小心选择，较大误差的测量值应被赋予较小的权。
并建立如下规则：被选择的参数，应该使算出的函数曲线与观测值之差的平方和最小。
用函数表示为：用欧几里得度量表达为：最小化问题的精度，依赖于所选择的函数模型。
系统辨识各类最小二乘法汇总
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 =
使用最小二乘法法进行系统辨识的两种方法
递推最小二乘法辨识与仿真现在有如下的辨识仿真对象：图中， )(k v 是服从N )1,0(分布的不相关随机噪声。
且)(1-zG )()(11--=z A z B ，)(1-z N )()(11--=zC zD ， （1）⎪⎪⎩⎪⎪⎨⎧=+==+-=--------1)(5.00.1)()(7.05.11)(121112111z D z zz B z C z z a z A选择上图所示的辨识模型。
仿真对象选择如下的模型结构：)()2()1()2()1()(2121k w k u b k u b k y a k y a k y +-+-=-+-+可得系统模型为：)()2(5.0)1()2(7.0)1(5.1)(k w k u k u k y k y k y +-+-=-+-- 递推最小二乘法的推导公式如下：)]1(ˆ)()()[()1(ˆ)(ˆ--+-=k k k z k k k θh K θθτ )1()]()([)]()()1([)(11--=+-=--k k k k k k k τP h K I h h PP τ1]1)()1()()[()1()(-+--=k k k k k k h P h h P K τ相关程序如下：% exp053 %%递推最小二乘法程序%clear%清理工作间变量L=55;% M序列的周期y1=1;y2=1;y3=1;y4=0;%四个移位寄存器的输出初始值for i=1:L;%开始循环，长度为Lx1=xor(y3,y4);%第一个移位积存器的输入是第3个与第4个移位积存器的输出的“或”x2=y1;%第二个移位积存器的输入是第3个移位积存器的输出x3=y2;%第三个移位积存器的输入是第2个移位积存器的输出x4=y3;%第四个移位积存器的输入是第3个移位积存器的输出y(i)=y4;%取出第四个移位积存器幅值为"0"和"1"的输出信号，if y(i)>0.5,u(i)=-0.03;%如果M序列的值为"1"时,辨识的输入信号取“-0.03”else u(i)=0.03;%当M序列的值为"0"时,辨识的输入信号取“0.03”end%小循环结束y1=x1;y2=x2;y3=x3;y4=x4;%为下一次的输入信号做准备end%大循环结束，产生输入信号uw=normrnd(0, sqrt(0.1), 1, 55);%加入白噪声figure(1);%第1个图形，伪随机序列stem(u),grid on%以径的形式显示出输入信号并给图形加上网格z(2)=0;z(1)=0;%取z的前两个初始值为零for k=3:55;%循环变量从3到55z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2)+w(k);%给出理想的辨识输出采样信号endc0=[0.001 0.001 0.001 0.001,0.001]';%直接给出被辨识参数的初始值,即一个充分小的实向量p0=10^6*eye(5,5);%直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.000000005;%相对误差E=0.000000005c=[c0,zeros(5,54)];%被辨识参数矩阵的初始值及大小e=zeros(5,55);%相对误差的初始值及大小for k=3:55; %开始求Kh1=[-z(k-1),-z(k-2),u(k-1),u(k-2),w(k)]'; 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<=E break;%若参数收敛满足要求，终止计算end%小循环结束end%大循环结束c;%显示被辨识参数e;%显示辨识结果的收敛情况%分离参数a1=c(1,:); a2=c(2,:); b1=c(3,:); b2=c(4,:);d1=c(5,:);ea1=e(1,:); ea2=e(2,:); eb1=e(3,:); eb2=e(4,:);figure(2);%第2个图形i=1:55;%横坐标从1到55plot(i,a1,'r',i,a2,':',i,b1,'g',i,b2,':',i,b1,'k') %画出a1，a2，b1，b2的各次辨识结果title('系统辨识结果')%图形标题如果系统中的参数为：⎪⎪⎩⎪⎪⎨⎧+-=+==+-=----------211211121112.01)(5.00.1)()(7.05.11)(z z z D z zz B z C z z a z A 那么系统模型机构为：)2(2.0)1()()2(5.0)1()2(7.0)1(5.1)(-+--+-+-=-+--k w k w k w k u k u k y k y k y 相关程序如下： % exp054.m % %递推最小二乘法编程%clear%清理工作间变量 L=55;% M 序列的周期y1=1;y2=1;y3=1;y4=0;%四个移位寄存器的输出初始值 for i=1:L;%开始循环，长度为Lx1=xor(y3,y4);%第一个移位积存器的输入是第3个与第4个移位积存器的输出的“或” x2=y1;%第二个移位积存器的输入是第3个移位积存器的输出 x3=y2;%第三个移位积存器的输入是第2个移位积存器的输出 x4=y3;%第四个移位积存器的输入是第3个移位积存器的输出 y(i)=y4;%取出第四个移位积存器幅值为"0"和"1"的输出信号，if y(i)>0.5,u(i)=-0.03;%如果M序列的值为"1"时,辨识的输入信号取“-0.03”else u(i)=0.03;%当M序列的值为"0"时,辨识的输入信号取“0.03”end%小循环结束y1=x1;y2=x2;y3=x3;y4=x4;%为下一次的输入信号做准备end%大循环结束，产生输入信号uw=normrnd(0, sqrt(0.1), 1, 55);figure(1);%第1个图形，伪随机序列stem(u),grid on%以径的形式显示出输入信号并给图形加上网格z(2)=0;z(1)=0;%取z的前两个初始值为零for k=3:55;%循环变量从3到55z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2)+w(k)-w(k-1)+0.2*w(k-2);%给出理想的辨识输出采样信号endc0=[0.001 0.001 0.001 0.001,0.001]';%直接给出被辨识参数的初始值,即一个充分小的实向量p0=10^6*eye(5,5);%直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.000000005;%相对误差E=0.000000005c=[c0,zeros(5,54)];%被辨识参数矩阵的初始值及大小e=zeros(5,55);%相对误差的初始值及大小for k=3:55; %开始求Kh1=[-z(k-1),-z(k-2),u(k-1),u(k-2),w(k),w(k-1),w(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<=E break;%若参数收敛满足要求，终止计算end%小循环结束end%大循环结束c;%显示被辨识参数e;%显示辨识结果的收敛情况%分离参数a1=c(1,:); a2=c(2,:); b1=c(3,:); b2=c(4,:);d1=c(5,:);d2=c(6,:);d3=c(7,:); figure(2);%第2个图形i=1:55;%横坐标从1到55plot(i,a1,'r',i,a2,':',i,b1,'g',i,b2,':',i,d1,'k',i,d2,':',i,d3,'*') %画出a1，a2，b1，b2的各次辨识结果title('系统辨识结果')%图形标题。
系统辨识之最小二乘
写成矩阵的形式为：
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
应用最小二乘一次完成法和递推最小二乘法算法的系统辨识.
目录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)应用最小二乘一次完成法和递推最小二乘法算法的系统辨识摘要：本题针对一个单输入单输出系统的便是问题，辨识的输入信号采用的是伪随机二位式序列（M序列），系统噪声为独立同分布高斯随机向量序列（白噪声），辨识的算法是递推最小二乘法和广义最小二乘法，本文简单描述应用经典辨识方法的辨识方案，详细描述了输入信号、噪声的产生方法及matlab程序，阐述了用两种不同算法的辨识原理并对它们的推导过程及辨识程序编制思路做了详细的描述。
最后结合真值与估计值对不同辨识算法的优劣进行了比较。
关键词：系统辨识M序列最小二乘法1引言1.1概述系统辨识是现代控制理论中的一个分支，它是根据系统的输入输出时间函数来确定描述系统行为的数学模型。
递推最小二乘辨识讲解
（2）每增加一次观测量，都必须重新计算φ，（）T -1。
（3）如果出现φ列相关，既不满秩的情况， T 为病态矩阵，则不能得到最 小二乘估计值。 *递推最小二乘参数辨识，就是当被辨识的系统在运行时，每取得一次新的观
测数据后，就在前一次估计结果的基础上，利用新引入的观测数据对前次 估计的结果，根据递推算法进行修正，从而递推地得出新的参数估计值。 这样，随着新的观测数据的逐次引入，一次接一次的进行参数计算，直到 参数估计值达到满意的精确程度为止。
P(k -1)(k -1) (k -1)
P
(k
)
I
-
1
(k
- 1) P
(k
-1)
(k
-1)
P
(k
-1)
其中的计算顺序为先计算P(k),然后再计算 θˆ(k) .
(7) (8)
有时,为计算方便并便于理解,上述RLS估计算法又可表示 为
常规最小二乘辨识的递推算法
主要内容
1.思想及原理 2.实例与MATLAB仿真 3.应用
1.递推最小二乘法的思想及原理
1.1递推最小二乘法的引入 *最小二乘法的缺陷 （1）数据量越多，系统参数估计的精度就越高，为了获得满意的辨识结果，
矩阵的阶数 T 常常取得相当大。这样矩阵求逆的计算量很大，存储量
1.2递推算法的思想 * 递推辨识算法的思想可以概括成
新的参数估计值=旧的参数估计值+修正项 即新的递推参数估计值是在旧的递推估计值的基础上而成，
这就是递推的概念. 递推算法不仅可减少计算量和存储量,而且能实现在线
实时辨识. * 递推算法是依时间顺序,每获得一次新的观测数据就修
正一次参数估计值,随着时间的推移,便能获得满意的辨 识结果. RLS法即为成批型LS算法的递推化，即将成批型LS算法
系统辨识第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
系统辩识实验报告
一、实验目的1. 理解系统辨识的基本概念和原理。
2. 掌握递推最小二乘算法在系统辨识中的应用。
3. 通过实验，验证算法的有效性，并分析参数估计误差。
二、实验原理系统辨识是利用系统输入输出数据，对系统模型进行估计和识别的过程。
在本实验中，我们采用递推最小二乘算法对系统进行辨识。
递推最小二乘算法是一种参数估计方法，其基本思想是利用当前观测值对系统参数进行修正，使参数估计值与实际值之间的误差最小。
递推最小二乘算法具有计算简单、收敛速度快等优点。
三、实验设备1. 电脑一台，装有MATLAB软件。
2. 系统辨识实验模块。
四、实验步骤1. 打开MATLAB软件，运行系统辨识实验模块。
2. 在模块中输入已知的系数a1、a2、b1、b2。
3. 生成输入序列u（t）和噪声序列v（t）。
4. 将输入序列u（t）和噪声序列v（t）加入系统，产生输出序列y（t）。
5. 利用递推最小二乘算法对系统参数进行辨识。
6. 将得到的参数估计值代入公式计算参数估计误差。
7. 仿真出参数估计误差随时间的变化曲线。
五、实验结果与分析1. 实验结果根据实验步骤，我们得到了参数估计值和参数估计误差随时间的变化曲线。
2. 结果分析（1）参数估计值：通过递推最小二乘算法，我们得到了系统参数的估计值。
这些估计值与实际参数存在一定的误差，这是由于噪声和系统模型的不确定性所导致的。
（2）参数估计误差：从参数估计误差随时间的变化曲线可以看出，递推最小二乘算法在短时间内就能使参数估计误差达到较低水平。
这说明递推最小二乘算法具有较好的收敛性能。
（3）参数估计误差曲线：在实验过程中，我们发现参数估计误差曲线在初期变化较快，随后逐渐趋于平稳。
这表明系统辨识过程在初期具有较高的灵敏度，但随着时间的推移，参数估计误差逐渐减小，系统辨识过程逐渐稳定。
六、实验结论1. 递推最小二乘算法在系统辨识中具有较好的收敛性能，能够快速、准确地估计系统参数。
2. 实验结果表明，递推最小二乘算法能够有效减小参数估计误差，提高系统辨识精度。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
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 ay 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θ=Φ()()ˆˆTT J e e y y θθ==-Φ-Φˆθ()ˆ20ˆT J y θθ∂=-Φ-Φ=∂ˆT T y θΦΦ=Φ()1ˆT T y θ-=ΦΦΦ220ˆTJ θ∂=ΦΦ>∂1n N yyyu T+-ΦΦ⎡⎤当N 足够大时有（4.1.8）如果矩阵ΦT Φ正定，则Ru 是是对称矩阵，并且各阶主子式的行列式为正。
当N 足够大时，矩阵Ru 才是是对称的。
由此引出矩阵ΦT Φ正定的必要条件是u(k)为持续激励信号。
如果序列{u(k)}的n+1阶方阵Ru 是正定的，则称{u(k)}的n+1阶持续激励信号。
下列随机信号都能满足Ru 正定的要求 1. 有色随机信号 2. 伪随机信号 3.白噪声序列1.1.2最小二乘估计的概率性质 最小二乘估计的概率性质 1） 无偏性由于输出值y 是随机的，所以θ是随机的，但注意θ不是随机值。
如果{}{}ˆE E θθθ==，则称ˆθ是θ无偏估计 2）一致性估计误差θ的方差矩阵为（4.1.9）上式表明当N →∞时，θ以概率1趋近于θ。
当ξ(k )为不相关随机序列时，最小二乘估计具有无偏性和一致性 3)有效性如果式（1.1.2）中的ξ的均值为0且服从正态分布的白噪声向量，则最小二乘参数估计值 为有效值，参数估计偏差的方差达到Cramer-Rao 不等式的下界，即（4.1.10）式中M 为Fisher 矩阵，且（4.1.11）4）渐近正态性如果式（4.1.2）中的ξ是均值为0且服从正态分布的白噪声向量，则最小二乘参数估计值ˆθ服从正态分布，即..11y uy W P T yu u R R R R R N ⎡⎤ΦΦ−−−→=⎢⎥⎣⎦121ˆT Var E NN σθ-⎧⎫⎪⎪⎛⎫=ΦΦ⎨⎬⎪⎝⎭⎪⎪⎩⎭21ˆlim lim 0,..1N N Var R w p Nσθ-→∞→∞==ˆθ(){}121ˆT Var E M θσ--=ΦΦ=()()ˆˆln /ln /ˆˆTp y p y M E θθθθ⎧⎫⎡⎤⎡⎤∂∂⎪⎪⎢⎥⎢⎥=⎨⎬⎢⎥⎢⎥∂∂⎪⎪⎣⎦⎣⎦⎩⎭12ˆT -1.2递推最小二乘法为了实现实时控制，必须采用递推算法，这种辨识方法主要用于在线辨识 设已获得的观测数据长度为N ，则估计方差矩阵为 式中 于是如果再获得1组新的观测值，则又增加1个方程得新的参数估计值 式中应用矩阵求逆引理，可得1N P +和N P 的递推关系式矩阵求逆引理 设A 为n ×n 矩阵，B 和C 为n ×m 矩阵，并且A ，A ＋BCT 和I+CTA-1B 都是非奇异矩阵，则有矩阵恒等式（4.2.2）得到递推关系式由于1TN N N P ψψ+是标量，因而上式可以写成（4.2.3）最后，得最小二乘法辨识公式（4.2.4）有2种给出初值的办法（1）设N0（N0>n ）为N 的初始值，可算出（4.2.5）N N NY θξ=Φ+()1ˆNT T N N N NY θ-=ΦΦΦ()122T N NN Var P θσσ-=ΦΦ=()1T N N N P -=ΦΦˆNT N N N P Yθ=Φ111T N N N y ψθξ+++=+()1111ˆT N N N N N N P Y y θψ++++=Φ+()11111111TN N TT T N N N N N N P P ψψψψ---+++++⎧⎫ΦΦ⎡⎤⎡⎤⎪⎪==+⎨⎬⎢⎥⎢⎥⎣⎦⎣⎦⎪⎪⎩⎭()()111111T T T A BC A A B I C A B C A ------+=-+()11111T TN N N N N N N N NP P P I P P ψψψψ-++++=-+()111111T T N N N N N N N N NP P P P P ψψψψ-++++=-+()10000000ˆ,TT N N N N N N N P P Y θ-=ΦΦ=Φ()1111ˆˆˆT N N N N N NK y θθψθ++++=+-()11111T N N N N N N K P P ψψψ-+++=+()1111111T T N N N N N N N N NP P P P P ψψψψ-+++++=-+（4.2.1）（2）假定200ˆ0,,P c I θ==C 是充分大的常数，I 为(2n+1) ×(2n+1)单位矩阵，则经过若干次递推之后能得到较好的参数估计。
[1][2][4]2两种算法的实现方案2.1最小二乘法一次完成算法实现如果把式（1.1.6）中的ˆθ取好足够的输入—输出数据以后一次计算出来，那么这种算法就式最小二乘法的一次完成法。
2.1.1最小二乘一次完成算法程序框图2.1.2一次完成法程序 具体程序参见附录42.1.3一次完成算法程序运行结果 ans =1.5000 0.7000 1.0000 0.5000 a1 = 1.5000 a2 = 0.7000 b1 =1.00002.1.4辨识数据比较2.1.5程序运行曲线-101输入u (k )0102030405060-505输出z (k)0102030405060-505离散输出z (k )2.2递推最小二乘法的实现2.2.1递推算法实现步骤1）输入M 序列的产生程序，可以参见3.2当中M 序列产生方法（具体程序见附录。
） 2）初始化。
一种简单的方法是选择ˆ0nθ=和2n P C I =，其中C 为尽量大的数，然后以它们为起始值进行计算(参照式2.2.3)。
可以证明，当C →∞时，按照递推公式算得的将与上面用批处理算法求得的结果相同，当N 很大时，C 的选择便无关紧要了。
3）递推算法的停机标准：()()()1ˆˆ1maxˆi i iik k k θθεθ-∀--<，ε是适当小的数。
4）最小二乘估计的递推算法系统参数的步骤如下：（1）产生系统输入信号M 序列，输入系统阶次，计算输入和输出数据u （i ），y （i ），i=1,2,…,n;（2）置N=n ，10ˆ0,10N NP I θ==（I …单位矩阵）； （3）形成行向量[]1()(1)()(1N y n y N n u N u N n ψ+=--+-+-（4）计算()()()111111;TTN N N K P N N P N ψψψ-+⎡⎤=++++⎣⎦（5）输入新测量数据y ( N+1 )和u （N+1）；（6）计算()11ˆˆˆ(1)1;N N N NK y N N θθψθ++⎡⎤=++-+⎣⎦（7）计算[]111;N N N N P I K P ψ+++=- （8）输出1ˆN θ+和1N P +； （9）N ←N+1,转（3）； 2.2.2程序编制思路：（1）产生M 序列，通过计算，依据算出输出值，设定递推初始值，采用4.2.4方法2设定辨识参数初值，0ˆθ为一个很小的量，P0为一个很大值的4×4矩阵，设定误差量，作为停机标准。
（2）依据设定，计算[]1(2)(1)(1)(2)Ty y u u ψ=--，1011T P ψψ⨯⨯+可以算出为一标量，依据式4.2.5算出K1,K1为4×1矩阵。
（3）依据式4.2.5计算出1ˆθ，和1P ，加入估计参数矩阵； （4）计算误差大小，求出相对变化率，加入误差矩阵第一列。
（5）再以1ˆθ和1P 为已知参数，重复（2）－（3）计算出参数2ˆθ，并加入估计参数矩阵。
（6）如此循环，计算出参数估计ˆnθ。
（7）判断误差变化率满足要求否，如果满足，则退出循坏，不再进行计算。
（8）分离估计参数，绘出图形。
2.2.3递推最小二乘法程序框图如下所示：具体程序参见附录2.2.4程序运行曲线-0.5-0.4-0.3-0.2-0.100.10.20.30.40.5图6M 序列信号的波形102030405060708090100-0.4-0.200.20.40.60.811.21.41.6Parameter Identification with Recursive Least Squares Method0102030405060708090100-400-2000200400600800100012001400Identification Precision2.2.5测试结果 见附录62.2.6地退数据表差的分布趋势。