系统辨识实验报告

系统辨识实验报告1一、相关分析法(1)实验原理)1)(1()s(11sTsTKG相关分析法u(k)Z(k))(ˆkg)(0kg)(~kg图1实验原理图本实验的原理图如图1。过程传递函数()Gs中12120,8.3,6.2KTSecTSec；输入变量()uk，输出变量()zk，噪声服从2(0,)vN，0()gk为过程的脉冲响应理论值，ˆ()gk为过程脉冲响应估计值，()gk为过程脉冲响应估计误差。过程输入()uk采用M序列，其输出数据加白噪声()vk得到输出数据()zk。利用相关分析法估计出过程的脉冲响应值ˆ()gk，并与过程脉冲响应理论值0()gk比较，得到过程脉冲响应估计误差值()gk。M序列阶次选择说明：首先粗略估计系统的过渡过程时间TS(通过简单阶跃响应)、截止频率fM(给系统施加不同周期的正弦信号或方波信号，观察输出)。本次为验证试验，已知系统模型，经计算HzTTfM14.0121，sTS30。根据式Mft3.0及式STtN)1(，则t取值为1，此时31N，由于t与N选择时要求完全覆盖，则选择六阶M移位寄存器，即N=63。系统辨识实验报告2(2)编程说明图2程序流程图(3)分步说明①生成M序列：M序列的循环周期63126N，时钟节拍1tSec，幅度1a，移位寄存器中第5、6位的内容按“模二相加”，反馈到第一位作为输入。其中初始数据设为{1,0,1,0,0,0}。程序如下：过程仿真得到理论输出数据()zk计算脉冲响应估计值计算互相关函数，得到脉冲响应估计值计算脉冲响应估计误差计算脉冲响应理论值，得到脉冲响应估计误差人机对话噪声标准差：sigma；生成数据周期数：r生成数据生成M序列()uk；生成白噪声序列()vk系统辨识实验报告3②生成白噪声序列：程序如下：③过程仿真得到输出数据：如图2所示的过程传递函数串联，可以写成形如121211()1/1/KGsTTsTsT，其中112KKTT。图2过程仿真方框图程序如下：④计算脉冲响应估计值：互相关函数采用公式)()(1)(10kiyixNrkRNrixy，互相关函数所用的数据是从第二个周期开始的，其中r为周期数，取1-3之间。则脉冲响应估计值为：])([1)(ˆckRkkgxy，taNNk2)1(。补偿量)1(NRcxy。程序如下：系统辨识实验报告4⑤计算脉冲响应估计值：脉冲响应的理论值由式12//012()[]ktTktTKgkeeTT可计算得到。这时可得到过程脉冲相应估计误差0ˆ()()()gkgkgk。脉冲响应估计误差为：NkNkgkgkg12012))(())(~(程序如下：(4)数据记录①当噪声标准差sigma=0.1，生成数据周期r为2时：脉冲响应估计误差为0.0121。脉冲响应估计曲线为图3所示。②当噪声标准差sigma=0.5，生成数据周期r为1时：脉冲响应估计误差为0.0347。脉冲响应估计曲线为图4所示。③当噪声标准差sigma=0.5，生成数据周期r为3时：脉冲响应估计误差为0.0258。脉冲响应估计曲线为图5所示。④当噪声标准差sigma=1，生成数据周期r为3时：脉冲响应估计误差为0.0279。脉冲响应估计曲线为图6所示。系统辨识实验报告5图3sigma=0.1，r=2时脉冲响应估计曲线图4sigma=0.5，r=1时脉冲响应估计曲线图5sigma=0.5，r=3时脉冲响应估计曲线图6sigma=1，r=3时脉冲响应估计曲线(5)结果分析实验中可以看到脉冲响应估计的曲线与理论曲线的重合度还是比较高的，脉冲响应估计误差也比较小，实验证明相关分析法的估计效果还是不错的。同时，经过实验可以得出结论：固定数据周期r，给定不同的噪声标准差sigma可以发现，噪声的方差越大，也就是信噪比越大，估计的效果越不好；固定噪声标准差sigma，选择不同的数据生成周期r可以发现，数据周期越大，估计的周期越多，估计的效果越好。系统辨识实验报告6二、最小二乘法1.基本最小二乘（离线辨识）baUYY][残差为：)(ˆ)()(kykyke最小二乘目标：残差平方和最小（一阶导为0，二阶导0）。正定，且TTTLSY1ˆ从上式看出，逆存在才有解，满足条件的u(k)：(1)伪随机；(2)白噪声；(3)有色随机信号。程序如下：结果如下：result1=[-0.8287;0.1275;-0.0024;1.9884;-1.2723]2.递推最小二乘（在线辨识）RLS)1()]()([)()()1()(1)()1()()]1(ˆ)()()[()1(ˆ)(ˆkPkkKIkPkkPkkkPkKkkkykKkkTTT修正项老估值新估值NNˆ)(ˆ)(1为了启动RLS，需给初值：IcP20,0ˆ。计算框图见书P66。程序如下：系统辨识实验报告7结果如下：result2=[-0.8284;0.1274;-0.0024;1.9883;-1.2717]图7递推最小二乘法参数过渡过程数据饱和：（1）原因：01NK，不再起修正作用，引起误差变大。（2）为了克服数据饱和现象，可以用降低老的数据影响的方法：①渐消记忆法（遗忘因子法）②限定记忆法（固定窗法）当)(k为不相关序列，最小二乘有一致性与无偏性，但)(k往往为相关序列，为克服最小二乘有偏估计的缺点，引入辅助变量法和广义最小二乘法，增广最小二乘法等。系统辨识实验报告83.渐消记忆法（遗忘因子法）)1()]()([)()()1()()()1()()]1(ˆ)()()[()1(ˆ)(ˆkPkkKIkPkkPkkkPkKkkkykKkkTTT一般99.0~95.0程序如下：结果如下：result3=[-0.6862;0.1114;0.0640;2.0356;-0.9429]图8渐消记忆法参数过渡过程4.限定记忆法（固定窗法）程序如下：系统辨识实验报告9结果如下：（较前三种方法偏差较大）result4=[-1.0190;0.1725;-0.4531;2.1092;-1.3724]5.辅助变量法（IV）(1).辅助变量ZYZZZZYZZZZYZZQQZEZETTTTTTTTTTT111)(ˆ)()()(噪声不相关Z,0)(其中，辅助变量估计值由强烈相关与为非奇异阵，，与(2).计算步骤：①先根据实测数据最小二乘求粗略ˆ（为有偏估计）②)(ˆ)()(ˆ)(ˆ)(ˆ11kykuzbkyza，得到辅助变量代入③n)()(ˆ)(),(Zkukykuky构造，用构造用④nTnTYZZ(3).递推：RIV)1()]()([)()()1()()()1()()]1(ˆ)()()[()1(ˆ)(ˆkPkzkKIkPkzkPkkzkPkKkkkykKkkTTT系统辨识实验报告10初始条件：IcP20,0ˆ缺点:P0的选择非常敏感，一个改进方法是，用递推最小二乘辨识算法作为启动方法，然后转换到辅助变量法。程序如下：结果如下：result5=[-0.9369;0.1207;-0.0254;1.9781;-1.5331]图9辅助变量法参数过渡过程系统辨识实验报告116.广义最小二乘法（GLS）(1)广义最小二乘法的基本思想：由于)(k在n+k个采样周期的时差范围内具有自相关性，从而使的最小二乘估计为有偏的，所以引入一个所谓成形滤波器（白化滤波器），把相关噪声)(k转化成白噪声)(k。如果知道有色噪声序列的相关性：)()()()(11kzdkzc令)()()(111zdzczf，有)()(1)(1kzfk广义最小二乘法(GLS)是建立在最小二乘法(LS)的基础上的。基本最小二乘法只是广义最小二乘法在1)(1zf时的特例。(2)广义最小二乘法计算步骤：广义最小二乘法的关键问题是如何用比较简单的方法找到成形滤波器的系数。其计算是逐次逼近法。①应用输入输出数据按最初模型求出的最小二乘估计。这个估计值是不精确的，它只是被估参数的一次近似。②计算残差e(k)，并拟合成形滤波器的模型：)(ˆ)1(ˆ)(ˆ)(ˆ)1(ˆ)(y)(101nkybkubkubnkyakyakkenn得到TmTTfffef]ˆˆˆ[][ˆ211其中)()1()1()(NmneNnemnene③应用所得的成形滤波器，对输入输出数据滤波：)(ˆ)1(ˆ)()()(ˆ)1(ˆ)()(11mkufkufkukumkyfkyfkykymm其中，m为噪声模型的阶，一般事先不知道，实际经验表明指定m为2或3可以得到比较满意的输出。④按新的输入、输出模型求出参数的第二次估计值ˆ。结果如下：Result6=[-0.6538;0.2926;0.0454;1.8776;-1.6042]7.广义递推最小二乘法（GLS）广义最小二乘法的递推计算过程可分成两个部分：(1)按递推最小二乘法(RLS)，随着N的增大，不断计算Nˆ(逐步接近于无偏)系统辨识实验报告12和Nfˆ(逐步使噪声白化)；(2)在递推过程中，Nˆ和Nfˆ是时变的，则过滤信号)(),(kyku及残差)(ke是由时变系统产生，要不断计算。因而，递推广义最小二乘法由两组普通的递推最小二乘法组成，它们是通过滤波算法联系起来的：1)1(1)1(1)1(11)1(1)1()1(11)1(11)1(1)1(111)1(11][11)ˆ(ˆˆTNNNTNNTNNNNNNNTNNNNNTNNNNNPPPPPPPPKyK1)2(1)2(1)2(11)2(1)2()2(11)2(11)2(1)2(111)2(11][11)ˆ(ˆˆTNNNTNNTNNNNNNNTNNNNNTNNNNNPPPPPPPPKfeKff结果如下：result7=[-1.0164;0.1754;-0.0159;2.0056;-1.6438]噪声传递系数的估计结果：[-0.0307;0.0900]图10广义递推最小二乘参数过渡过程8.增广矩阵法(ELS/RELS)(增广最小二乘法)增广矩阵法是把观测矩阵适当增大，使得有偏估计的程度得到一定改善。])()1()()()()1([NNnNuNnuNyNnyTN这一方程结构适用于递推最小二乘法，但向量TN中)(k是未知的。解决这个矛盾的一个方法是用估计值)(ˆk代替)(k，即：])(ˆ)1(ˆ)()()()1([ˆNNnNuNnuNyNnyTN其中，1ˆˆ)()(ˆNTNNnyNn，)(ˆNn初始值可以取为零。按照与前面类似的递推方法，得到增广矩阵法参数估计。增广矩阵法在实际中有广泛的应用，收敛情况也比较好。系统辨识实验报告13结果如下：result8=[-0.8636;0.0747;-0.0162;1.9729;-1.1420]图10增广最小二乘参数过渡过程