系统辨识最小二乘法大作业

系统辨识大作业最小二乘法及其相关估值方法应用学院：自动化学院专业：信息工程学号：**********姓名：***日期：2010.11.14基于最小二乘法的多种系统辨识方法研究1.最小二乘法的引出在系统辨识中用得最广泛的估计方法是最小二乘法(LS)。

设单输入-单输出线性定长系统的差分方程为x(k)+a1x(k−1)+⋯+a n(k−n)=b0u(k)+⋯+b n u(k−n)，k=1，2，3，⋯(5.1.1)式中: u(k)为随机干扰；x(k)为理论上的输出值。

x(k)只有通过观测才能得到，在观测过程中往往附加有随机干扰。

x(k)的观测值y(k)可表示为y(k)=x(k)+n(k)(5.1.2)式中: n(k)为随机干扰。

由式(5.1.2)得x(k)=y(k)−n(k)(5.1.3)将式(5.1.3)带入式(5.1.1)得y(k)+a1y(k−1)+⋯+a n y(k−n)n=b0u(k)+b1u(k−1)+⋯+b n u(k−n)+n(k)+∑a i(k−i)i=1(5.1.4)我们可能不知道n(k)的统计特性，在这种情况下，往往把n(k)看做均值为0的白噪声。

设nξ(k)=n(k)+∑a i(k−i)i=1(5.1.5)则式(5.1.4)可写成y(k)=−a1y(k−1)−a2y(k−2)−⋯−a n y(k−n)+b0u(k)+b1u(k−1)+⋯+b n u(k−n)+ξ(k)(5.1.6)在观测u(k)时也有测量误差，系统内部也可能有噪声，应当考虑它们的影响。

因此假定ξ(k)不仅包含了x(k)的测量误差，而且包含了u(k)的测量误差和系统内部噪声。

假定ξ(k)是不相关随机序列(实际上ξ(k)是相关随机序列)。

现分别测出n+N个随机输入值y(1)，y(2)，⋯，y(n+N)，u(1)，u(2)，⋯，u(n+ N)，则可写成N个方程，即y(n+1)=−a1y(n)−a2y(n−1)−⋯−a n y(1)+b0u(n+1)+b1u(n)+⋯+b n u(1) +ξ(n+1)y(n+2)=−a1y(n+1)−a2y(n)−⋯−a n y(2)+b0u(n+2)+b1u(n+1)+⋯+b n u(2) +ξ(n+2)⋮y(n+N)=−a1y(n+N−1)−a2y(n+N−2)−⋯−a n y(N)+b0u(n+N)+b1u(n+N−1)+⋯+b n u(N)+ξ(n+N)上述N个方程可写成向量-矩阵形式[y(n+1) y(n+2)⋮y(n+N)]=[−y(n)−y(n+1)⋯⋯−y(1)u(n+1)⋯u(1)−y(2)u(n+2)⋯u(2)⋮⋮ ⋮ ⋮⋮−y(n+N−1)⋯−y(N)u(n+N)⋯u(N)]×[a1⋮a nb0⋮b n]+[ξ(n+1)ξ(n+2)⋮ξ(n+3)](5.1.7)设y=[y(n+1) y(n+2)⋮y(n+N)]，θ=[a1⋮a nb0⋮b n]，ξ=[ξ(n+1)ξ(n+2)⋮ξ(n+N)]Φ=[−y(n)−y(n+1)⋯⋯−y(1)u(n+1)⋯u(1)−y(2)u(n+2)⋯u(2)⋮⋮ ⋮ ⋮⋮−y(n+N−1)⋯−y(N)u(n+N)⋯u(N)]则式(5.1.7)可写为y=Φθ+ξ(5.1.8)式中：y为N维输出向量；ξ为N维噪声向量；θ为(2n+1)维参数向量；Φ为N×(2n+1)测量矩阵。

因此式(5.1.8)是一个含有(2n+1)个未知参数，由N个方程组成的联立方程组。

如果N<2n+1，方程数少于未知数数目，则方程组的解是不定的，不能唯一地确定参数向量。

如果N=2n+1，方程组正好与未知数数目相等，当噪声ξ=0时，就能准确地解出θ=Φ−1y(5.1.9)如果噪声ξ≠0，则θ=Φ−1y−Φ−1ξ(5.1.10)从上式可以看出噪声ξ对参数估计是有影响的，为了尽量较小噪声ξ对θ估值的影响。

在给定输出向量y和测量矩阵Φ的条件下求系统参数θ的估值，这就是系统辨识问题。

可用最小二乘法来求θ的估值，以下讨论最小二乘法估计。

2.最小二乘法估计算法设θ̂表示 θ的最优估值，ŷ表示 y的最优估值，则有ŷ=Φθ̂(5.1.11)ŷ=[ŷ(n+1) ŷ(n+2)⋮ŷ(n+N)]，θ̂=[a1̂⋮a n̂b0̂⋮b n̂]写出式(5.1.11)的某一行，则有ŷ(k)=−a1̂y(k−1)−a2̂y(k−2)−⋯−a n̂y(k−n)+b0̂u(k)+b1̂u(k−1)+⋯+b n̂u(k−n)+ξ(k)=−∑a îy(k−i)ni=1+∑b îu(k−i)ni=0，k=n+1，n+2，⋯，n+N(5.1.12) 设e(k)表示y(k)与ŷ(k)之差，即e(k)= y(k)-ŷ(k)=y(k)−ŷ(k)=y(k)—∑a îy(k−i)ni=1+∑b îu(k−i)ni=0=(1+a1̂z−1+⋯+a n̂z−n)y(k)−(b0̂+b1̂z−1+⋯+b n̂z−n)u(k)=â(z−1)y(k)−b̂(z−1)u(k)，k=n+1，n+2，⋯，n+N(5.1.13)式中â(z−1)=1+a1̂z−1+⋯+a n̂z−nb̂(z−1)=b0̂+b1̂z−1+⋯+b n̂z−ne(k)成为残差。

把k=n+1，n+2，⋯，n+N分别代入式(5.1.13)可得残差e(n+1)，e(n+2)，⋯，e(n+N)。

设e=[e(n+1)e(n+2)⋯e(n+N)]T则有e(k)= y−ŷ=y−Φθ̂(5.1.14) 最小二乘估计要求残差的平方和为最小，即按照指数函数J=e T e=(y−Φθ̂)T(y−Φθ̂)(5.1.15) 为最小来确定估值θ̂。

求J对θ̂的偏导数并令其等于0可得∂J∂θ̂=−2ΦT(y−Φθ̂)=0(5.1.16)ΦTΦθ̂=ΦT y(5.1.17) 由式(5.1.17)可得θ的最小二乘估计θ̂=(ΦTΦ)−1ΦT y(5.1.18) 3.递推最小二乘法为了实现实时控制，必须采用递推算法，这种辨识方法主要用于在线辨识。

设已获得的观测数据长度为N，将式(5.1.8)中的y和ξ分别用Y N，ΦN，ξ̅N来代替，即Y N=ΦNθ+ξ̅N(5.3.1) 用θ̂N表示 θ的最小二乘估计，则θ̂N=(ΦN TΦN)−1ΦN T Y N(5.3.2)设P N =（ΦN TΦN ）−1(5.3.5)于是 θ̂N =P N ΦN T Y N (5.3.6) 如果再获得1组新的观测值u (n +N +1)和y(n +N +1)，则又增加1个方程y N+1=ψN+1T θ+ξN+1(5.3.7)式中y N+1=y (n +N +1)，ξN+1=ξ(n +N +1)ψN+1T =[−y (n +N ) ⋯ −y (N +1) u (n +N +1) ⋯ u(N +1)]将式(5.3.1)和式(5.3.7)合并，并写成分块矩阵形式，可得[Y N y N+1]=[ΦN ψN+1T ] θ+[ξ̅N ξN+1](5.3.8)根据上式可得到新的参数估值θ̂N+1={[ΦN ψN+1T ]T [ΦNψN+1T ]}−1[ΦN ψN+1T ]T [Y N y N+1]=P N+1[ΦN ψN+1T ]T [Y N y N+1]=P N+1(ΦN T Y N +ψN+1y N+1)(5.3.9)式中P N+1={[ΦN ψN+1T ]T [ΦNψN+1T ]}−1=（ΦN TΦN +ψN+1ψN+1T）−1根据矩阵求逆引理可以求得递推最小二乘法辨识公式θ̂N+1=θ̂N +K N+1(y N+1−ψN+1T θ̂N )(5.3.19) K N+1=P N+1ψN+1(1+ψN+1T P N ψN+1)−1(5.3.20)P N+1=P N −P N ψN+1(1+ψN+1T P N ψN+1)−1ψN+1T P N(5.3.21)由于进行递推计算需要给出θ̂N 和P N 初值P 0和θ̂0，通过计算证明，可以取初值：θ̂0=0，P 0=c 2I ，c 是充分大的常数，I 为(2n +1)(2n +1)单位矩阵，则经过若干次递推之后能够得到较好的参数估计。

3. 辅助变量法辅助变量法是一种可克服最小二乘有偏估计的一种方法，对于原辨识方程y =Φθ+ξ(5.4.1) 当ξ(k)是不相关随机序列时，最小二乘法可以得到参数向量θ的一致无偏估计。

但是，在实际应用中ξ(k)往往是相关随机序列。

假定存在着一个(2n +1)×N 的矩阵Z 满足约束条件{limN→∞1NZ Tξ=0limN→∞1NZ TΦ=Q(5.4.2)式中Q是非奇异的。

用Z T乘以式(5.4.1)等号两边得Z T y=Z TΦθ+Z T ξ(5.4.3) 由上式得θ=(Z TΦ)−1Z T y−(Z TΦ)−1Z Tξ(5.4.4) 如果取θ̂IV=(Z TΦ)−1Z T y(5.4.5) 作为θ估值，则称θ̂IV为辅助变量估值，矩阵Z成为辅助变量矩阵，Z中的元素称为辅助变量。

常用的辅助变量法有递推辅助变量参数估计法，自适应滤波法，纯滞后等。

4.广义最小二乘法广义最小二乘法是能克服最小二乘法有偏估计的另一种方法，这种方法计算比较复杂但效果比较好。

下面直接介绍广义最小二乘法的计算步骤：(1)应用得到的输入和输出数据u(k)和y(k)(k=1，2，3，⋯，n+N)，按模型a(z−1)y(k)=b(z−1)u(k)+ξ(k)求出θ 的最小二乘估计θ̂(1)=[a1̂(1)⋮a n̂(1) b0̂(1)⋮b n̂(1)](2)计算残差e(1)(k)e(1)(k)=â(1)(z−1)y(k)−b̂(1)(z−1)u(k)(3)用残差e(1)(k)代替 ξ(k)，计算f̂(1)=[(Ω(1))TΩ(1)]−1(Ω(1))T e(1)(4)计算y̅(1)(k)和u̅(1)(k)y̅(1)(k)=y(k)+f̂1(1)y(k−1)+⋯+f̂m(1)y(k−m)u̅(1)(k)=u(k)+f̂1(1)u(k−1)+⋯+f̂m(1)u(k−m)(5)应用得到的y̅(1)(k)和u̅(1)(k)按模型a(z−1)y̅(1)(k)=b(z−1)u̅(1)(k)+ε(k)用最小二乘法重新估计θ，得到θ的第2次估值θ̂(2)。

然后按步骤(2)计算残差e(2)(k)，按步骤(3)重新估计f，得到估值f̂(2)。

再按照步骤(4)计算y̅(2)(k)和u̅(2)(k)，按照步骤(5)求θ的第3次估值θ̂(3)。

重复上述循环，之道θ的估值θ̂(i)收敛为止。

5.一种交替的广义最小二乘法求解技术(夏式法)这种方法是夏天长提出来的，又称夏式法。

以上讨论过的广义最小二乘法的特点在于系统的输入和输出信号反复过滤。