系统辨识试验

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) b1u(k 1) b2u(k 2) v(k) 设输入信号的取值是从 k =1 到 k =16 的 M序列，则待辨识参数?LS 为?LS =(H LH L ) 1 H LzL 。其中，被辨识参数 ?LS 、观测矩阵 z L、 HL的表达式为a1z(3)z(2)z

2、(1)u(2)u(1)?a2 ,LSz(4) ，zL， HLz(3)z(2)u(3)u(2)LSb1b2z(16)z(15)z(14)u(15)u(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) % 画三行一列图形窗口

3、中的第一个图形stem(u) % 画出输入信号 u的经线图形 subplot(3,1,2) % 画三行一列图形窗口中的第二个图形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(

4、3) u(4) u(3);-z(5) -z(4) u(5) u(4);-z(6) -z(5) u(6) u(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)

5、;z(5);z(6);z(7);z(8);z(9);z(10);z(11);z(12);z(13);z(14);z(15);z(16)% 给样本矩阵 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注：由于输出观测值没有任何噪音成分， 所以辨识结果也无任何误差， 同学们可以在输 出观测值中添加噪音，观

6、察 ols 的辨识效果。同时，可以尝试增加输入信号的数量，看辨识 结果有何变化。实验二 基于 RLS法的系统辨识数字仿真实验一、实验目的1、深入理解系统辨识中相关分析法及最小二乘法的相关内容2、学会用 Matlab 或 C 语言等进行系统辨识的仿真研究二、实验设备装有相应软件的计算机。三、实验原理1. 考虑如下图所示的仿真对象：N(z1)v(k)e(k)1+G(z )u(k)y(k)图 1 递推最小二乘算法辨识z(k)图中， v(k)是服从N(0,1) 分布的不相关随机噪声。且111)G(z 1) B(z ) ，N (z 1)D(z )A(z 1) C (z 1)A(z 1) 1 1.5a1z

7、 1 0.7z 2 C(z 1)B(z 1) 1.0z 1 0.5z 2D(z 1) 1选择上图所示的辨识模型。仿真对象选择如下的模型结构：z(k) a1z(k 1) a2 z(k 2) b1u(k 1) b2u(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) ；加权阵取单位阵 L I ；利用如下公式计算 K(k)、?(k)和 P( k) ，计算各次参数辨识的相

8、对误差，精度满足要求后停机。递推最小二乘法的推导公式如下：(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) P(k) k(k 1)x (k 1)P(k)2. 阶的辨识前面所讨论的系统辨识方法， 都是假定模型的阶次是已知的， 因此仅仅要求估计差分方 程的系数。 但实际上， 系统的阶次是很难被准确知道的。 因为对阶次的理解程度是直接与一 个线性差分方程的准确结构有关的， 所以有关阶次的确定也可以称为系统结构的确定。 经验 指出， 一个模型的阶次不准， 就可能在控制系统设计时发生严重问题

9、。 故在辨识过程中， 模 型的阶次是否合适是必须加以检验的。 一般阶的方法中， 常用的有这么几种： 零极点相消法、 目标函数法和 F 检验法。下面只介绍其中的目标函数法。当我们用不同阶的模型给系统的输入输出观测数据进行最小二乘拟合时，会得到 不同的估计误差：N2TJe2(k) ET Ek1因此利用 J 极小化确定阶是很自然的。实验表明，假设模型具有大于 1 而小于 Nmax 的 阶 n，当 n 取 1,2, 时 , 若随着 n 的增加 , 在 n? ( 阶的估计量 )-1 时， J 最后一次出现陡峭 的下降，往后 J就近似地保持不变或者只有微小的下降 (见下图), 则取n n?。也就是说， 模

10、型阶次的确定可以直接依次计算阶次n 1，2，时的最小二乘估计 n?以及相应的损失函数 J ，然后选择当 J 下降不明显时的阶次作为合适的模型阶次n ，这种方法也叫确定阶的估计准则方法，有很广的应用。四、实验内容1. 用递推最小二乘法（ RLS）法辨识对象数学模型。2. 对象阶的辨识。五、实验要求1. 熟悉系统辨识中的相关内容。2. 掌握 Matlab 或 C 语言等进行系统辨识仿真研究的一般步骤。3. 实验前基本应完成相关的编程任务，实验时调试相应程序。4. 修改相应参数与随机噪声幅度，观察并分析结果。5. 软件包人机界面的开发与设计。 （选做）六、实验步骤1. 首先要熟悉一下 MATLAB的

11、运行环境 :1）File->New->M-File 打开 M文件编辑窗口2）输入自己编写的程序3）点击 run 按钮，如果程序出错则调试程序，如果运行正常的话则观察程序的 运行结果2. 用递推最小二乘法（ RLS）法辨识对象数学模型在这个实验中，我们采用以下模型进行仿真：y=1.5*yk-1-0.7*yk-2+0*uk+1.0*uk-10.5*uk-2+ek(5)图2 递推最小二乘法算法辨识的 Malab程序流程图其中 uk 是幅值为 1 的 PRBS信号输入， ek 是白噪声，即（ 0， 1）的正态分布序列，它的方差时可以调整的。这个系统的采样值y（ k）和 u（k） 作为已知数

12、据，采用实验一的最 小二乘法估计这个系统的参数。面是递推最小二乘法的流程图：相关程序如下：%RLSclear% 清理工作间变量L=15;% M 序列的周期y1=1;y2=1;y3=1;y4=0;% 四个移位寄存器的输出初始值for i=1:L;% 开始循环，长度为 Lx1=xor(y3,y4);% 第一个移位积存器的输入是第 3 个与第 4 个移位积存器的输出的 “或”x2=y1;%第二个移位积存器的输入是第3 个移位积存器的输出x3=y2;%第三个移位积存器的输入是第2 个移位积存器的输出x4=y3;%第四个移位积存器的输入是第3 个移位积存器的输出y(i)=y4;% 取出第四个移位积存器幅

13、值为"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;% 为下一次的输入信号做准备 end%大循环结束，产生输入信号 ufigure(1);% 第 1 个图形stem(u),grid on% 以径的形式显示出输入信号并给图形加上网格z(2)=0;z(

14、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.001 0.001 0.001 0.001'% 直接给出被辨识参数的初始值 , 即一个充分小的实向 量p0=106*eye(4,4);% 直接给出初始状态 P0，即一个充分大的实数单位矩阵E=0.000000005;% 相对误差 E=0.000000005c=c0,zeros(4,14);% 被辨识参数矩阵的初始值及大小e=zeros(4,

15、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; % 把当前相对变化的列向量加入误差矩阵的最后一列c0=c1;% 新获得的参数作为下一次递推的旧参数c(:,k)=c1;% 把辨识参数

16、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 个图形i=1:15;% 横坐标从 1 到 15plot(i,a1,'

17、;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

18、 的各次辨 识结果的收敛情况title('Identification Precision') %图形标题注：同样这个程序使用的输出信号也没有噪音， 所以辨识的结果没有误差， 请同学们在输出信号中加入噪音，再使用 RLS对其辨识，观察辨识结果，进行分析。2. 同样采用这个模型采用 MATLAB或者 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） 作为已知数据，采用实验一的最小二乘法估计这个系统的参数。令模型的阶次分别为n1， 2，3，其 J 值如下表所示。从表中看出，由 n? 1到n? 2 ，J发生显著的变化。而由 n? 2到n? 3，J 没有多大的变化， 故可确定这系统的阶次 n? 2辨识的结果如下，具体编程由同学们来实现。噪声水平损失函数 Jn1n2n30.0265.8630.0000.1248.1470.9870.9830.5337.8482