系统辨识大作业

一、 问题描述考虑仿真对象：()0.9(1)0.15(2)0.02(3)0.7(1)0.15(2)()z k z k z k z k u k u k e k +-+-+-=---+ e() 1.0e(1)0.41e(2)(),~(0,1)k k k v k v N λ+-+-=式中，u(k)和z(k)是输入输出数据，v(k)是零均值、方差为1的不相关的随机噪声；u(k)采用与e(k)不相关的随机序列。

1. 设计实验，产生输入输出数据；2. 使用基本最小二乘法估计参数；3. 考虑其他适用于有色噪声的辨识方法估计参数；4. 模型验证。

二、 问题分析对于单输入单输出系统（Single Input Single Output, SISO ），如图 1所示，将待辨识的系统看作“灰箱”，它只考虑系统的输入输出特性，而不强调系统的内部机理。

图 1中，输入u(k)和输出z(k)是可以测量的，1()G z -是系统模型，用来描述系统的输入输出特性，y(k)是系统的实际输出。

1()N z -是噪声模型，v(k)是均值为零的不相关随机噪声，e(k)是有色噪声。

图 1 SISO 系统的“灰箱”结构对于SISO 随机系统，被辨识模型()G z 为：12121212()()()1n n n n b z b z b z y z G z u z a z a z a z------+++==++++L L 其相应的差分方程为11()()()n ni i i i y k a y k i b u k i ===--+-∑∑若考虑被辨识系统或观测信息中含有噪声，被辨识模型可改写为11()()()()n ni i i i z k a y k i b u k i v k ===--+-+∑∑式中，z(k)为系统输出量的第k 次观测值；y(k)为系统输出量的第k 次真值，y(k-1)为系统输出量的第k-1次真值，以此类推；u(k)为系统的第k 个输入值，u(k-1)为系统的第k-1个输入值；v(k)为均值为0的不相关随机噪声。

1. 数据生成本部分需要生成系统的输入输出数据以及噪声数据。

1）白噪声的生成辨识数据通常包含有噪声，如果该噪声相关性较弱或者强度很小，可近似看作白噪声。

本次实验问题中()(0,1)v t N :，即服从标准正态分布，可以将噪声看作为服从正态分布的白噪声过程，在Matlab 中可以由randn 函数生成。

2）输入数据的生成伪随机二进制序列（Pseudo Random Binary Sequence, PRBS ）是广泛应用的一种伪随机序列，所谓“二进制”是指序列中每个随机变量只有“0”或“1”两种逻辑状态。

伪随机二进制序列可由多级线性反馈移位寄存器组成的随机信号发生器产生，其中具有最长循环周期的线性移位寄存器序列是伪随机二进制序列最常见的一种形式，简称M 序列（Maximal Length Sequence ）。

M 序列由于具有近似白噪声的性质，而且工程上易于实现，能够保证较好的系统辨识精度，是普遍采用的一种辨识用输入信号。

图 2 线性反馈移位寄存器产生伪随机二进制序列结构图以一个4级线性反馈移位寄存器产生伪随机二进制序列为例，如图 2所示。

假设4个移位寄存器0123,,,a a a a 输出的初态非全零，移位寄存器的工作原理是：一个移位脉冲来到后，每级移位寄存器的输出移到下一级移位寄存器作为输入，最末一级移位寄存器的输出即为伪随机二进制序列。

3）输出数据的生成根据给定的SISO 系统，可以求出z(k)的表达式：()0.9(1)0.15(2)0.02(3)0.7(1)0.15(2)v() 1.0(1)0.41(2)z k z k z k z k u k u k k e k e k λ=------+---+---- 其理想系数值为12312120.9,0.15,0.02,0.7, 1.5, 1.0,0.41a a a b b c c =====-==.可以根据生成的白噪声序列和输入序列，以及必要的0初始值，带入表达式即可得到采样输出数据。

2. 差分模型阶检验在实际场景中，辨识模型的阶数和纯时延往往是未知的，在很多情况下仅仅依靠猜测。

在模型的阶数和纯时延不确定时，设系统模型为11()()()()n n i j i j y t a y t i b y t i t ζ===--+-+∑∑其中n 为模型的阶数，1()()()t C z e t ζ-=模型的阶估计可以采用多种方法，本实验采用比较简单易行的损失函数检验法。

定义预报误差（噪声方差的估计值）的平方和为损失函数，即11ˆ()N N i J t N ε==∑当n 从小增大时，N J 应随之减小，当n 增大到某一值时，N J 应近似白噪声过程。

采用以下的检验原则：在n-1这一点，N J 最后一次出现陡峭的下降，此后就近似地保持不变或只有微小的下降，则取ˆnn =。

3. 参数辨识模型在系统辨识和参数估计领域中，最广泛的估计方法时最小二乘法和极大似然估计法。

最小二乘法作为一种最基本的估计方法应用极为广泛，其他的大多数估计算法都与最小二乘法有关。

它既可用于动态系统，也可用于静态系统；可用于线性系统，也可用于非线性系统；可用于离线估计，也可用于在线估计。

在随机环境下利用最小二乘法时，无须知道观测数据的概率统计信息，而这种方法获得的结果，却有相当好的统计性质。

最小二乘参数估计方法来源于数理统计的回归分析，它能提供一个在最小方差意义上与实验数据最好的你和的模型。

该估计在一定条件下有最佳的统计特性，即它们是一致的，无偏的和有效的。

最小二乘法时一个经典的方法，概念简明，适应范围广，在一些情况下，可得到与极大似然法一样好的统计效果，它能很方便地与其他辨识算法建立关系。

（1）递推最小二乘算法当获得一批数据后，可一次求得相应的参数估计值，这样处理问题的方法称为一次完成算法或批处理算法。

它在理论研究方面有许多方便之处，但当矩阵的维数增加时，矩阵求逆运算的计算量会急剧增加，将给计算机的计算速度和存储量带来负担，而且不适合在线辨识，无法跟踪参数随时间变化的情况。

为了减少计算量，减少数据在计算机中所占的存储量，也为了实时地辨识出动态系统的特性，在用最小二乘法进行参数估计时，把它转化成参数递推的估计。

参数递推估计是指被辨识的系统，每取得一次新的测量数据后，就在前一次估计结果的基础上，利用新引入的测量数据对前一次估计的结果进行修正，从而递推地得出新的参数估计值。

这样，随着新测量数据的引入，一次接一次地进行参数估计，直到估计值达到满意的精确程度为止。

最小二乘递推算法的基本思想可以概括为：当前的估计值ˆ()k θ=上次估计值ˆ(1)k θ-+修正项 即新的估计值ˆ()k θ是在旧的估计值ˆ(1)k θ-的基础上，利用新的观测数据对旧的估计值进行修正而成的。

考虑如下模型：11()()()()()A z z k B z u k v k --=+，其中()u k ，()z k 分别是系统的输入和输出；()v k 是均值为零，方差为一的不相关白噪声。

且满足：1121211212()1......()......nn m m A z a z a z a z B z b z b z b z--------⎧=++++⎨=+++⎩ 令1212()[(1),...,(),(1),...()][,,...,,,...]Tn m h k z k z k n u k u k m a a a b b b θ⎧=------⎨=⎩则使用加权最小二乘参数估计递推算法（Recursive Weighted Least Squares, RWLS ）有： 1ˆˆˆ()(1)()[()()(1)]1()(1)()[()(1)()]()()[()()](1)()[(1),...(),(1),...,()]T T T T k k K k z k h k k K k P k h k h k p k h k k P k I K k h k P k h k z k z k n u k u k m θθθ-⎧=-+--⎪⎪=--+⎪Λ⎨⎪=--⎪⎪=------⎩图 3 最小二乘递推辨识参数估计过程中信息的变换可以看出，取()1k Λ=的时候，加权最小二乘估计就退化成了最小二乘参数估计的递推算法（Recursive Least Squares, RLS ）。

加权参数1Λ可以在(0,1]范围内选择，如果选11()k =Λ，所有的采样数据都是等同加权的，如果11()k <<Λ,则表示对新近获得的数据给予充分大的加权因子，而削弱历史观测数据的影响。

实际计算时，需要首先确定初始参数ˆ(0)θ和(0)P 。

一般说来选取 2(0)(0)P I αθε⎧=⎨=⎩αε为充分大实数为充分小的向量对于这样的系统，使用最小二乘法参数估计的递推算法进行辨识可以得到无偏估计，但是如果噪声模型必须用1()()C z v k -表示时，此时就无法得到无偏估计了，因为该方法没有把噪声模型考虑进去。

使用递推最小二乘法进行系统辨识的过程如图4所示。

（2）增广最小二乘法当噪声均值为0时，最小二乘参数估计算法为无偏估计；当噪声的均值不为0时，最小二乘参数估计算法为有偏估计。

为了解决最小二乘参数估计的有偏性，将噪声模型的辨识同时考虑进去，因此称为增广最小二乘法。

该算法可以看成是对一般最小二乘参数辨识算法的简单推广或扩充，因此又称为扩充最小二乘算法。

考虑如下模型：111()()()()()()A z z k B z u k C z v k ---=+，其中()u k ，()z k 分别是系统的输入和输出；()v k 是均值为零，方差为一的不相关白噪声。

且满足：112121*********()1...()...()1...nn m m vv A z a z a z a z B z b z b z b zC z c z c z c z ------------⎧=++++⎪=+++⎨⎪=++++⎩在模型阶次,,n m v 已经确定的情况下，令121212()[(1),...,(),(1),...(),(1),...()][,,...,,,...,,,...,]TT n m v h k z k z k n u k u k m v k v k v a a a b b b c c c θ⎧=--------⎨=⎩由于()v k 是白噪声，所以利用增广最小二乘法这一形式上的变换，即可获得参数θ的无偏估计。

不过由于数据向量()h k 中包含着不可测的噪声量(1),(2),...()v k v k v k v ---，因此要用相应的估计值代替。