系统辨识试验

1、2.用普通最小二乘法(OLS法辨识对象数学模型选择的仿真对象的数学模型如下z(k) 1.5z(k 1)0.7z(k 2) u(k 1)0.5u(k 2) v(k)其中，v(k)是服从正态分布的白噪声N(0,1)。输入信号采用4阶M序列，幅度为1。选择如下形式的辨识模型z(k) a1z(k 1)a2z(k 2) bk 1)b2u(k 2) v( k)设输入信号的取值是从k =1到k =16的M序列，则待辨识参数Js为?ls=(HlTHl) 1HlTzl。其中，被辨识参数Js、观测矩阵z L、Hl的表达式为31z(3)z(2)z(1)u(2)u(1)? a2 ,JLSz(4)uZl， Hlz(3)

2、z(2)u(3)u(2)b1b2z(16)zez(14)ugu(14程序框图如下所示:图2最小二乘一次完成算法程序框图参考程序：%ols u=-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1; % 系统辨识的输入信号为一个周期的 M 序列z=zeros(1,16); % 定义输出观测值的长度for k=3:16z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); % 用理想输出值作为观测值 endsubplot(3,1,1) % 画三行一列图形窗口中的第一个图形stem(u) %画出输入信号u的经线图形subplot(3,1,2)

3、% 画三行一列图形窗口中的第二个图形 i=1:1:16; % 横坐标范围是 1到 16，步长为 1plot(i,z) % 图形的横坐标是采样时刻 i, 纵坐标是输出观测值 z, 图形格式为连续曲线 subplot(3,1,3) % 画三行一列图形窗口中的第三个图形stem(z),grid on%画出输出观测值z的经线图形，并显示坐标网格u,z%显示输入信号和输出观测信号%L=14数据长度HL=-z(2) -z(1) u(2) u(1);-z(3)-z(2) u(3) u(2);-z(4)-z(3) u(4) u(3);-z(5)-z(4)u(5) u(4);-z(6) -z(5) u(6) u

4、(5);-z(7) -z(6) u(7) u(6);-z(8) -z(7) u(8) u(7);-z(9)-z(8) u(9) u(8);-z(10)-z(9)u(10) u(9);-z(11)-z(10) u(11) u(10);-z(12)-z(11) u(12)u(11);-z(13)-z(12)u(13)u(12);-z(14)-z(13)u(14) u(13);-z(15)-z(14)u(15)u(14) %给样本矩阵HL赋值ZL=z(3);z(4);z(5);z(6);z(7);z(8);z(9);z(10);z(11);z(12);z(13);z(14);z(15);z(16)%

5、 给样本矩阵zL赋值%calculating parameters% 计算参数c1=HL*HL; c2=inv(c1); c3=HL*ZL; c=c2*c3 % 计算并显示%DISPLAY PARAMETERSa1=c(1), a2=c(2), b1=c(3), b2=c(4) % 从 中分离出并显示 a1 、a2、 b1、 b2%End注：由于输出观测值没有任何噪音成分， 所以辨识结果也无任何误差， 同学们可以在输 出观测值中添加噪音，观察 ols 的辨识效果。同时，可以尝试增加输入信号的数量，看辨识 结果有何变化。实验二 基于RLS法的系统辨识数字仿真实验一、实验目的1、深入理解系统辨识中

6、相关分析法及最小二乘法的相关内容2、学会用Matlab或C语言等进行系统辨识的仿真研究二、实验设备装有相应软件的计算机。三、实验原理1. 考虑如下图所示的仿真对象：N(z1)e(k)G(z )u(k)y(k)v(k)图1递推最小二乘算法辨识z(k)图中，v(k)是服从N(0,1)分布的不相关随机噪声。且(1)1 1G (z 1 )BtZ L ， N (z 1 ) D (z IA(z 1)C(z 1)A(z 1)1 1.5a1z 10.7z 2 C(z 1)B(z 1)1.0z 10.5z 2D(z 1) 1选择上图所示的辨识模型。仿真对象选择如下的模型结构：z(k) a1z(k 1) a2z(

7、k 2) b1u(k 1)pu(k 2) v(k)( 2)其中，v(k)是服从正态分布的白噪声N(0,1)。输入信号采用 4位移位寄存器产生的M序列，幅度为0.03。按式(3)z(k) 1.5z(k 1)0.7z(k 2) u(k 1) 0.5u(k 2) v(k) (3)构造h(k)；加权阵取单位阵 A,I ；利用如下公式计算K(k)、$k)和P(k)，计算各次参数辨识的相对误差，精度满足要求后停机。递推最小二乘法的推导公式如下：(k 1) (k) k(k 1) y(k 1) xT(k 1) (k)K(k 1) P(k)x(k 1)1 x (k 1)P(k)x(k 1) (4)P(k 1)

8、P(k) k(k 1)x (k 1)P(k)2. 阶的辨识前面所讨论的系统辨识方法， 都是假定模型的阶次是已知的， 因此仅仅要求估计差分方 程的系数。 但实际上， 系统的阶次是很难被准确知道的。 因为对阶次的理解程度是直接与一 个线性差分方程的准确结构有关的， 所以有关阶次的确定也可以称为系统结构的确定。 经验 指出， 一个模型的阶次不准， 就可能在控制系统设计时发生严重问题。 故在辨识过程中， 模 型的阶次是否合适是必须加以检验的。 一般阶的方法中， 常用的有这么几种： 零极点相消法、 目标函数法和 F 检验法。下面只介绍其中的目标函数法。当我们用不同阶的模型给系统的输入输出观测数据进行最小

9、二乘拟合时，会得到 不同的估计误差：N2TJe2(k) ET Ek1因此利用 J 极小化确定阶是很自然的。实验表明，假设模型具有大于1 而小于 Nmax 的阶n,当n取1,2,时，若随着n的增加，在?(阶的估计量)-1时，J最后一次出现陡峭 的下降，往后 J 就近似地保持不变或者只有微小的下降 ( 见下图 ), 则取 n n? 。也就是说， 模型阶次的确定可以直接依次计算阶次n= 1, 2，时的最小二乘估计 ？以及相应的损失函数J ,然后选择当J下降不明显时的阶次作为合适的模型阶次 n,这种方法也叫确定阶的估 计准则方法,有很广的应用。J四、实验内容1. 用递推最小二乘法（RLS法辨识对象数学

10、模型。2. 对象阶的辨识。五、实验要求1. 熟悉系统辨识中的相关内容。2. 掌握Matlab或C语言等进行系统辨识仿真研究的一般步骤。3. 实验前基本应完成相关的编程任务，实验时调试相应程序。4. 修改相应参数与随机噪声幅度，观察并分析结果。5. 软件包人机界面的开发与设计。（选做）六、实验步骤1. 首先要熟悉一下 MATLAB的运行环境：1）File-New-M-File 打开M文件编辑窗口2）输入自己编写的程序3）点击run按钮，如果程序出错则调试程序，如果运行正常的话则观察程序的 运行结果2. 用递推最小二乘法（RLS法辨识对象数学模型在这个实验中，我们采用以下模型进行仿真：+ 0.5*

11、uk-2+eky=1.5*yk-1-0.7*yk-2+0*uk+1.0*uk-1图2递推最小二乘法算法辨识的Malab程序流程图其中uk是幅值为1的PRBS言号输入，ek是白噪声，即(0, 1)的正态分布序列,它的方差时可以调整的。这个系统的采样值y (k)和u(k)作为已知数据，采用实验一的最 小二乘法估计这个系统的参数。面是递推最小二乘法的流程图：相关程序如下：%RLSclear% 清理工作间变量L=15;% M 序列的周期y1=1;y2=1;y3=1;y4=0;% 四个移位寄存器的输出初始值for i=1:L;% 开始循环，长度为 Lx1=xor(y3,y4);% 第一个移位积存器的输入

12、是第 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;% 为下一次的输入信号做准备en d%大循环

13、结束，产生输入信号ufigure(1);% 第 1 个图形stem(u),grid on% 以径的形式显示出输入信号并给图形加上网格 z(2)=0;z(1)=0;% 取 z 的前两个初始值为零for k=3:15;% 循环变量从 3 到 15z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2);%给出理想的辨识输出采样信号end%RLS递推最小二乘辨识c0=0.0010.001 0.001 0.001;% 直接给出被辨识参数的初始值 ,即一个充分、的实向量p0=10A6*eye(4,4);%直接给出初始状态 P0，即一个充分大的实数单位矩阵E=0.00000

14、0005;% 相对误差 E=0.000000005c=c0,zeros(4,14);% 被辨识参数矩阵的初始值及大小e=zeros(4,15);% 相对误差的初始值及大小for k=3:15; % 开始求 Kh1=-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;% 求被辨识参数 ce1=c1-c0;% 求参数当前值与上一次的值的差值e2=e1./c0;% 求参数的相对变化e(:,k)=e2; % 把当前相对变化的列

15、向量加入误差矩阵的最后一列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,:);ea1=e(1,:);ea2=e(2,:);eb1=e(3,:);eb2=e(4,:);figure(2);% 第 2 个图

16、形i=1:15;% 横坐标从 1 到 15plot(i,a1,r,i,a2,:,i,b1,g,i,b2,:)%画出 a1，a2，b1，b2 的各次辨识结果title(Parameter Identification with Recursive Least Squares Method)%图形标题figure(3); % 第 3 个图形i=1:15; % 横坐标从 1 到 15plot(i,ea1,r,i,ea2,g,i,eb1,b,i,eb2,r:) %画出 a1，a2，b1，b2 的各次辨识结果的收敛情况title(Identification Precision) %图形标题注：同样这个

17、程序使用的输出信号也没有噪音， 所以辨识的结果没有误差， 请同学们在输出信号中加入噪音，再使用RLS对其辨识，观察辨识结果，进行分析。2. 同样采用这个模型采用MATLA或者C语言进行阶的辨识：y=1.5*yk-1-0.7*yk-2+0*uk+1.0*uk-1+ 0.5*uk-2+ek其中uk是幅值为1的PRBS言号输入，ek是白噪声，即(0, 1)的正态分布序列， 它的方差时可以调整的。这个系统的采样值y (k)和u(k)作为已知数据，采用实验一的最小二乘法估计这个系统的参数。令模型的阶次分别为n= 1, 2, 3，其J值如下表所示。从表中看出，由n? 1到n? 2, j发生显著的变化。而由 ?2到n? 3, j没有多大的变化,故可确定这系统的阶次 ? 2辨识的结果如下，具体编程由同学们来实现。噪声水平损失函数Jn = 1n= 2n= 3= 0.0265.8630.000=0.1248.1470.9870.983= 0.