系统辨识大作业

《系统辨识大作业》由会员分享，可在线阅读，更多相关《系统辨识大作业（21页珍藏版）》请在装配图网上搜索。

1、系统辨识大作业递推增广最小二乘法题目一： 已知一系统为两输入单输出系统，观测数据受有色噪声污染，噪信比为N/S=0.1。系统经2000次采样，存放于文件T3T.TXT中。系统输入u1为7级M序列，u2为u1的63步移位序列。模型类可选为：A(q-1)y(k)=B1(q)u1(k)+ B2(q)u2(k)+w(k)/C(q-1)。要求编制程序，辨识出该模型的结构及参数。（注：可以将模型变形，以适合算法）作业文档要求：描述问题；选择辨识方法并简单说明所选方法中的结构辨识原理和参数估计原理；程序流程图及程序清单；说明程序中用到的一些技术，如数据标准化、UD分解、稳定性判断等；结构搜索路线及各结构下的

2、参数、残差；给出最终结果：A(q-1)=B1(q)=B2(q)=C(q-1)=并给出选择此最终结果的理由；用你的辨识结果来预报系统输出误差e(k)=y(k)-(k)，并画出e(101)-e(400)的曲线图。1.问题描述已知一系统为两输入单输出系统，观测数据受有色噪声污染，噪信比为N/S=0.1。系统经2000次采样，存放于文件T3T.TXT中。系统输入u1为7级M序列，u2为u1的63步移位序列。模型类可选为：A(q-1)y(k)=B(q-1)u1(k)+ C(q-1)u2(k)+ D(q-1)w(k)。2. 辨识原理和参数估计原理2.1递推增广最小二乘法的辨识原理 递推增广最小二乘法(简称

3、RELS方法)是用于控制系统参数估计和结构辩识的一种方法,它基于岱( 最小二乘) 方法,使被控对象数学模型在误差信号平方和最小的意义上由实验数据拟合,同时把有色噪声看成是由白噪声合成的, 从而解决了最小二乘法算法的有偏性和非一致性问题。其中n=0；n任意时是广义递推最小二乘法就是递推增广最小二乘法。2.3 结构辨识在系统辨识中，系统的阶次是一个十分重要的参数，但在实际的情况中往往无法提前获知系统的阶数。这就需要在参数辨识之前先通过结构辨识，来确定阶数。一般来说经典的结构辨识方法主要有以下3种：F检验定阶法，AIC准则法，FPE准则法。（1）F检验定阶法假设系统的阶次为（），计算出模型输出和观测

4、输出的残差平方和函数 (18)其中为观测输出，为模型输出。一般来说会随着模型阶次的增加而减少。实际进行系统辨识时，输入输出的信号采样的个数大大超过待辨识的参数个数，在这种情况下，随着模型阶次的增加值先是显著下降，当模型的阶次大于真实系统的阶次时，的显著下降现象就中止。利用这个原理可以判定模型应有的阶次。（2）AIC准则法赤池信息量准则，即Akaike information criterion、简称AIC，是衡量统计模型拟合优良性的一种标准，是由日本统计学家赤池弘次创立和发展的。赤池信息量准则建立在熵的概念基础上，可以权衡所估计模型的复杂度和此模型拟合数据的优良性。在一般的情况下，AIC可以表

5、示为： (19)其中：是参数的数量，是似然函数。 假设条件是模型的误差服从独立正态分布。设为观察数，RSS为剩余平方和，那么AIC变为：增加自由参数的数目提高了拟合的优良性，AIC鼓励数据拟合的优良性但是尽量避免出现过度拟合(Overfitting)的情况。所以优先考虑的模型应是AIC值最小的那一个。赤池信息准则的方法是寻找可以最好地解释数据但包含最少自由参数的模型。（3）FPE准则法对AR模型的一种定阶法.设是平稳AR(p)序列，是样本，是样本自协方差函数， (j=1，2，n)是拟合AR(n)参数的最小二乘估计. （20）为最终预报误差.使FPE(n)达到最小的n称为AR(p)模型阶数p的估

6、计.此方法是赤池弘次于1969年提出的所谓改进的残差方差图方法.最终预报误差准则名称的来源是对AR(p)序列观测值拟合k阶AR模型，使得一步预报均方误差达到最小，只有k取真阶p才能达到.所以定阶依据这种准则是合乎实际目的的。3 程序流程图3.1 结构辨识结果在参数辨识之前首先进行结构辨识，使用F检验法。取最后200个数据计算残差平方和函数值。具体数值如表1所示表1 F 检验法n Jt（n-1，n）26.639830.8680645.005340.2363256.636550.23450.7292系统模型为。模型阶次每增加1，待辨识参数数目就增加4。查询F分布表可知。观察表1中t的值，当阶次为4

7、时，当阶次为5时，。根据F检验法的思想，在模型阶次为5时，残差平方函数值显著下降的现象中止，所以实际系统的阶次为5。3.2 程序流程图及程序清单模型完成定阶以后，就开始参数的辨识。由于系统的阶次为5，所以模型的递推表达式可以写成 （21）仿真程序的流程图如图1所示。清零产生有色噪声信号画出有色噪声序列图加载T3t.txt文件并分别赋给z，u1，u2给辨识参数和P赋初值按照式(17)的第二式计算K(k)按照式(17)的第一式计算(k)按照式(17)的第三式计算P(k)计算被辨识参数的相对变化量参数是否满足要求？分离参数画出被辨识参数的各次递推估计式图形画出被辨识参数的相对误差的图形停机NY图1

8、递推输出误差辨识算法的程序流程图递推增广最小二乘法的辨识程序如下：%递推增广最小二乘法clearL = 2000;n = 5;size = 5+4*(n-1);v = sqrt(0.1)*randn(2000,1);M = load(T3t.txt);z(:,1) = M(1:L,1);u1(:,1) = M(:,2);u2(:,1) = M(:,3);%递推增广最小二乘法c0=0.0001*ones(size,1);%直接给出被辨识参数的初始值,即一个充分小的实向量p0=1000000*eye(size,size);%直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.00005;%取相

9、对误差Ec=c0,zeros(size,L-1);%被辨识参数矩阵的初始值及大小e=zeros(size,L);%相对误差的初始值及大小zm = zeros(L,1);ek = zeros(L,1);for k=6:L; %开始求K h1=-z(k-1,1),-z(k-2,1),-z(k-3,1),-z(k-4,1),-z(k-5,1),u1(k-1,1),u1(k-2,1),u1(k-3,1),u1(k-4,1),u1(k-5,1),u2(k-1,1),u2(k-2,1),u2(k-3,1),u2(k-4,1),u2(k-5,1),v(k),v(k-1),v(k-2),v(k-3),v(k-

10、4),v(k-5);%为求K(k)作准备 x=h1*p0*h1+1; x1=inv(x); k1=p0*h1*x1; %K d1=z(k)-h1*c0; c1=c0+k1*d1;%辨识参数c e1=c1-c0; e2=e1./c0; %求参数的相对变化 e(:,k)=e2; c0=c1;%给下一次用 c(:,k)=c1;%把辨识参数c 列向量加入辨识参数矩阵 p1=p0-k1*k1*h1*p0*h1+1; %p1=p0-k1*h1*p0; p0=p1;%给下次用 end%循环结束for k =6:Lh2=-zm(k-1,1),-zm(k-2,1),-zm(k-3,1),-zm(k-4,1),-

11、zm(k-5,1),u1(k-1,1),u1(k-2,1),u1(k-3,1),u1(k-4,1),u1(k-5,1),u2(k-1,1),u2(k-2,1),u2(k-3,1),u2(k-4,1),u2(k-5,1),v(k),v(k-1),v(k-2),v(k-3),v(k-4),v(k-5); zm(k) = h2*c(:,L); ek(k,1) = z(k,1) - zm(k);endJ = ek*ek;J,c(:,L)%显示被辨识参数及参数收敛情况figure(1);k =1:L;plot(k,z,r,k,zm,b);xlabel(k);ylabel(幅值);title(观测输出与模

12、型输出);figure(2);k=101:400;plot(k,ek(101:400,1),b);xlabel(k);ylabel(幅值);title(残差);a1 = c(1,:); a2 = c(2,:); a3 = c(3,:); a4 = c(4,:); a5 = c(5,:);figure(3);subplot(3,1,1);k = 1: L;plot(k,a1(1,1:L),b);xlabel(k);ylabel(幅值);title(a1);subplot(3,1,2);k = 1: L;plot(k,a2(1,1:L),b);xlabel(k);ylabel(幅值);title(

13、a2);subplot(3,1,3);k = 1: L;plot(k,a3(1,1:L),b);xlabel(k);ylabel(幅值);title(a3);figure(4);subplot(2,1,1);k = 1: L;plot(k,a4(1,1:L),b);xlabel(k);ylabel(幅值);title(a4);subplot(2,1,2);k = 1: L;plot(k,a5(1,1:L),b);xlabel(k);ylabel(幅值);title(a5);b1 = c(6,:); b2 = c(7,:); b3 = c(8,:); b4 = c(9,:); b5 = c(10

14、,:);figure(5);subplot(3,1,1);k = 1: L;plot(k,b1(1,1:L),b);xlabel(k);ylabel(幅值);title(b1);subplot(3,1,2);k = 1: L;plot(k,b2(1,1:L),b);xlabel(k);ylabel(幅值);title(b2);subplot(3,1,3);k = 1: L;plot(k,b3(1,1:L),b);xlabel(k);ylabel(幅值);title(b3);figure(6);subplot(2,1,1);k = 1: L;plot(k,b4(1,1:L),b);xlabel(

15、k);ylabel(幅值);title(b4);subplot(2,1,2);k = 1: L;plot(k,b5(1,1:L),b);xlabel(k);ylabel(幅值);title(b5);c1 = c(11,:); c2 = c(12,:); c3 = c(13,:); c4 = c(14,:); c5 = c(15,:);figure(7);subplot(3,1,1);k = 1: L;plot(k,c1(1,1:L),b);xlabel(k);ylabel(幅值);title(c1);subplot(3,1,2);k = 1: L;plot(k,c2(1,1:L),b);xla

16、bel(k);ylabel(幅值);title(c2);subplot(3,1,3);k = 1: L;plot(k,c3(1,1:L),b);xlabel(k);ylabel(幅值);title(c3);figure(8);subplot(2,1,1);k = 1: L;plot(k,c4(1,1:L),b);xlabel(k);ylabel(幅值);title(c4);subplot(2,1,2);k = 1: L;plot(k,c5(1,1:L),b);xlabel(k);ylabel(幅值);title(c5);d0 = c(16,:); d1 = c(17,:); d2 = c(18

17、,:); d3 = c(19,:); d4 = c(20,:);d5 = c(21,:);figure(9);subplot(3,1,1);k = 1: L;plot(k,d0(1,1:L),b);xlabel(k);ylabel(幅值);title(d0);subplot(3,1,2);k = 1: L;plot(k,d1(1,1:L),b);xlabel(k);ylabel(幅值);title(d1);subplot(3,1,3);k = 1: L;plot(k,d2(1,1:L),b);xlabel(k);ylabel(幅值);title(d2);figure(10);subplot(3

18、,1,1);k = 1: L;plot(k,d3(1,1:L),b);xlabel(k);ylabel(幅值);title(d3);subplot(3,1,2);k = 1: L;plot(k,d4(1,1:L),b);xlabel(k);ylabel(幅值);title(d4);subplot(3,1,3);k = 1: L;plot(k,d5(1,1:L),b);xlabel(k);ylabel(幅值);title(d5);4 程序中采用的技术4.1 稳定性判断递推算法失败的一个原因是模型在修正过程中变为不稳定模型。解决的办法是在每次修正后判断模型的稳定性，若模型稳定，则接受修正，否则，将

19、模型的修正量减半。Jurry-Astrom稳定判据：设线性定常离散系统的特征方程为 A(z)=a0zn+a1zn-1+.+an-1z+an=0 特征方程稳定的充分必要条件是：a0 0, b0 0, . r0 0, q0 0 以上判据很容易编制成程序。 4.2 程序代码解释根据题目要求噪声的信噪比为0.1，输入信号的功率为1，故噪声信号的功率为0.1。v = sqrt(0.1)*randn(2100,1);这条语句的作用是产生了均值为0，方差为1，长度为2000的高斯白噪声序列。参数辨识的主要过程为公式（17）的两个迭代式，题目中一共给出了2000个数据，程序中的主要计算部分就是利用这2000个

20、数据按式（17）进行迭代。主要代码如下所示：for k=6:L; %开始求K h1=-z(k-1,1),-z(k-2,1),-z(k-3,1),-z(k-4,1),-z(k-5,1),u1(k-1,1),u1(k-2,1),u1(k-3,1),u1(k-4,1),u1(k-5,1),u2(k-1,1),u2(k-2,1),u2(k-3,1),u2(k-4,1),u2(k-5,1),v(k),v(k-1),v(k-2),v(k-3),v(k-4),v(k-5);%为求R(k)作准备r1 = r0 + rou(k)(h1*h1-r0) ;c1=c0+rou(k)inv(r1)h1(z(k)-h1*

21、c0);r0=r1;c0=c1; end%循环结束迭代过程结束后，c矩阵的第2100列的元素即为辨识出的参数，利用辨识出的参数可以计算观测输出和模型输出的差，进一步可以求出残差平方和函数，改变的值可以求得不同阶次下的残差平方和函数，由函数可以判断实际系统的阶次 （22）5 残差输出及参数估计5.1 残差输出观测输出与模型输出的曲线图如图2所示。图2 观测输出与模型输出5.2 参数估计输出图3是a1-a5的估计输出图图3 a1-a5的参数估计结果输出图4是b1-b5的参数估计输出。图4 b1-b5的参数估计输出图5是c1-c5的参数估计输出图5 c1-c5的参数估计输出图6是d0-d5的参数估计

22、输出图6 d0-d5的参数估计输出由图3-图6参数估计结果可得到系统辨识的最终结果。系统辨识的最终结果如下：6 预报系统输出误差由辨识结果来预报系统输出误差e(k)=y(k)-(k)，其e(101)-e(400)的曲线图如下图7所示。图7 e(101)-e(400)的曲线图题目二： 数据文件heating_system.dat来自网上，以下是它的描述：1. Contributed by: Roy Smith Dept. of Electrical & Computer Engineering University of California,Santa Barbara, CA 93106U.S

23、.A. royece.ucsb.edu2. Process/Description:The experiment is a simple SISO heating system.The input drives a 300 Watt Halogen lamp, suspendedseveral inches above a thin steel plate. The outputis a thermocouple measurement taken from the back ofthe plate.3. Sampling interval: 2.0 seconds4. Number of s

24、amples：8015. Inputs: u: input drive voltage6. Outputs: y: temperature (deg. C)7. References:The use of this experiment and data for robustcontrol model validation is described in: Sampled Data Model Validation: an Algorithm and Experimental Application, Geir Dullerud & Roy Smith, International Journ

25、al of Robust and Nonlinear Control, Vol. 6, No. 9/10, pp. 1065-1078, 1996.请自选算法，辨识结构与参数（注意它的数据排列为：t,u,y）。1问题描述一个单输入单输出加热系统，输入驱动一盏300w的碘钨灯。碘钨灯悬挂在离薄钢板上方几英寸。系统的输出是薄钢板背面的热电偶温度计（单位摄氏度），温度计的采样间隔为2秒，采样数据数目为801个。请请自选算法，辨识系统结构与参数。2辨识原理仍考虑系统（1）的参数辨识问题，选取准则函数 (23) 极小化，求得参数的估计值，可以使模型的输出最好地预报系统的输出。可以解得使极小化的 (24)

26、由公式（32）可一次完成系统参数辨识，但是该方法需要较大的内存，另外对输入信号要求比较苛刻，所以仅在理论分析中使用，在实际操作中常使用递推最小二乘法，递推最小二乘法的公式为 （25）3.仿真程序假设系统模型为，其中为零均值，方差为1 的高斯白噪声。首先对模型进行定阶，采用F检验法计算模型各阶次的残差平方和函数，数值如表2所示表2 F检验法数值123模型阶次每增加1，待辨识参数数目就增加4。查询F分布表可知，。观察表1中t的值，当阶次为2时， ，当阶次为3时， 。根据F检验法的思想，在模型阶次为3时，残差平方函数值显著下降的现象中止，所以实际系统的阶次为2。3.1流程图与代码清单仿真程序流程图如

27、图8所示结束画图分离参数数据是否用完按式（25）第三式计算按式（25）第一式计算按式（25）第二式计算给参数赋初值读入数据清零NY图8程序流程图程序主要代码清单如下：clear%清理工作间变量n = 2;size= 2*n;L=801;% M序列的周期data = load(heating_system.dat);z(:,1) = data(:,3);u(:,1) = data(:,2);%RLS递推最小二乘辨识c0=0.001*ones(size,1);%直接给出被辨识参数的初始值,即一个充分小的实向量p0=106*eye(size,size);%直接给出初始状态P0，即一个充分大的实数单位

28、矩阵E=0.000000005;%相对误差E=0.000000005c=c0,zeros(size,L-1);%被辨识参数矩阵的初始值及大小e=zeros(size,L);%相对误差的初始值及大小zm = zeros(size,1);for k=n+1:L; %开始求K h1=-z(k-1),-z(k-2),u(k-1),u(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;%求被辨识参数c e1=c1-c0;%求参数当前值与上一次的值的差值 e2=e1./c0;%求参数的

29、相对变化 e(:,k)=e2; %把当前相对变化的列向量加入误差矩阵的最后一列 c0=c1;%新获得的参数作为下一次递推的旧参数 c(:,k)=c1;%把辨识参数c 列向量加入辨识参数矩阵的最后一列 p1=p0-k1*k1*h1*p0*h1+1;%求出 p(k)的值 p0=p1;%给下次用end%大循环结束for k =n+1:L h2 = -zm(k-1),-zm(k-2),u(k-1),u(k-2); zm(k) = h2*c(:,L);endfor k=1:L ek(k,1) = z(k)-zm(k);endfigure(1);k= 1:L;plot(k,z,r,k,zm,b);xlab

30、el(k);ylabel(幅值);title(观测输出与模型输出);figure(2);k= 1:L;plot(k,ek,b);xlabel(k);ylabel(幅值);title(残差);程序的主要计算部分就是根据式（33）进行迭代，这部分代码如下所示：for k=n+1:L; %开始求K h1=-z(k-1),-z(k-2),u(k-1),u(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;%求被辨识参数c e1=c1-c0;%求参数当前值与上一次的值的差值 e2=e

31、1./c0;%求参数的相对变化 e(:,k)=e2; %把当前相对变化的列向量加入误差矩阵的最后一列 c0=c1;%新获得的参数作为下一次递推的旧参数 c(:,k)=c1;%把辨识参数c 列向量加入辨识参数矩阵的最后一列 p1=p0-k1*k1*h1*p0*h1+1;%求出 p(k)的值 p0=p1;%给下次用end%大循环结束3.2结果与讨论系统观测输出和模型输出的对比图如图9所示，残差图如图10所示图9 观测输出和模型输出图10 残差图最终系统的辨识结果为:各系数的仿真结果如图11-12所示图11 系数a1-a2仿真图图12 系数b1-b2仿真图各待辨识系数在辨识过程中相对误差的变化图如图13-14所示图13 系数a1-a2相对误差图17系数b1-b2相对误差图