2.用普通最小二乘法（OLS ）法辨识对象数学模型 选择的仿真对象的数学模型如下)()2(5.0)1()2(7.0)1(5.1)(k v k u k u k z k z k z +-+-=-+--其中，)(k v 是服从正态分布的白噪声N )1,0(。

输入信号采用4阶M 序列，幅度为1。

选择如下形式的辨识模型)()2()1()2()1()(2121k v k u b k u b k z a k z a k z +-+-=-+-+设输入信号的取值是从k =1到k =16的M 序列，则待辨识参数LSθˆ为LS θˆ=L τL 1L τL z H )H H -(。

其中，被辨识参数LSθˆ、观测矩阵z L 、H L 的表达式为 ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=2121ˆb b a a LSθ , ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=)16()4()3(z z z L z ， ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡------=)14()2()1()15()3()2()14()2()1()15()3()2(u u u u u u z z z z z z LH 程序框图如下所示：参考程序：%olsu=[-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) %画三行一列图形窗口中的第二个图形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(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)]% 给样本矩阵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、深入理解系统辨识中相关分析法及最小二乘法的相关内容。

2、学会用Matlab 或C 语言等进行系统辨识的仿真研究二、实验设备装有相应软件的计算机。

三、实验原理1. 考虑如下图所示的仿真对象：图中， )(k v 是服从N )1,0(分布的不相关随机噪声。

且)(1-z G )()(11--=z A z B ，)(1-z N )()(11--=z C z D ， （1）⎪⎪⎩⎪⎪⎨⎧=+==+-=--------1)(5.00.1)()(7.05.11)(121112111z D z z z B z C z z a z A 选择上图所示的辨识模型。

仿真对象选择如下的模型结构：)()2()1()2()1()(2121k v k u b k u b k z a k z a k z +-+-=-+-+ （2）其中，)(k v 是服从正态分布的白噪声N )1,0(。

输入信号采用4位移位寄存器产生的M 序列，幅度为0.03。

按式(3))()2(5.0)1()2(7.0)1(5.1)(k v k u k u k z k z k z +-+-=-+-- (3)构造h (k )；加权阵取单位阵I Λ=L ；利用如下公式计算K (k )、)(ˆk θ和P (k )，计算各次参数辨识的相对误差，精度满足要求后停机。

递推最小二乘法的推导公式如下：(1)()()(1)(1)1(1)()(1)(1)()(1)(1)()(1)()(1)[][1](1)T k k k k k T k k k k k k T k k k k k k y xK P x x P x P P k k x P θθθ∧∧∧+++-++++++=++-=+=-+ （4）2.阶的辨识前面所讨论的系统辨识方法，都是假定模型的阶次是已知的，因此仅仅要求估计差分方程的系数。

但实际上，系统的阶次是很难被准确知道的。

因为对阶次的理解程度是直接与一个线性差分方程的准确结构有关的，所以有关阶次的确定也可以称为系统结构的确定。

经验指出，一个模型的阶次不准，就可能在控制系统设计时发生严重问题。

故在辨识过程中，模型的阶次是否合适是必须加以检验的。

一般阶的方法中，常用的有这么几种：零极点相消法、目标函数法和F 检验法。

下面只介绍其中的目标函数法。

当我们用不同阶的模型给系统的输入——输出观测数据进行最小二乘拟合时，会得到不同的估计误差：因此利用J 极小化确定阶是很自然的。

实验表明，假设模型具有大于1而小于maxN 的阶n ，当n 取1,2,…时, 若随着n 的增加, 在ˆn(阶的估计量)-1时，J 最后一次出现陡峭的下降，往后J 就近似地保持不变或者只有微小的下降(见下图),则取ˆn n=。

也就是说，模型阶次的确定可以直接依次计算阶次n ＝1，2，…时的最小二乘估计ˆn以及相应的损失函数J ，然后选择当J 下降不明显时的阶次作为合适的模型阶次n ，这种方法也叫确定阶的估计准则方法，有很广的应用。

∑===Nk T EE k e J 12)(Jn四、实验内容1. 用递推最小二乘法（RLS）法辨识对象数学模型。

2. 对象阶的辨识。

五、实验要求1. 熟悉系统辨识中的相关内容。

2. 掌握Matlab或C语言等进行系统辨识仿真研究的一般步骤。

3. 实验前基本应完成相关的编程任务，实验时调试相应程序。

4. 修改相应参数与随机噪声幅度，观察并分析结果。

5. 软件包人机界面的开发与设计。

（选做）六、实验步骤1.首先要熟悉一下MATLAB的运行环境:1)File->New->M-File打开M文件编辑窗口2)输入自己编写的程序3)点击run按钮，如果程序出错则调试程序，如果运行正常的话则观察程序的运行结果2. 用递推最小二乘法（RLS）法辨识对象数学模型在这个实验中，我们采用以下模型进行仿真：y=1.5*y[k-1]-0.7*y[k-2]+0*u[k]+1.0*u[k-1]＋0.5*u[k-2]+e[k] (5)其中u[k]是幅值为1的PRBS信号输入，e[k]是白噪声，即（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);%第一个移位积存器的输入是第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;%为下一次的输入信号做准备end%大循环结束，产生输入信号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.001 0.001 0.001 0.001]';%直接给出被辨识参数的初始值,即一个充分小的实向量p0=10^6*eye(4,4);%直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.000000005;%相对误差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; %把当前相对变化的列向量加入误差矩阵的最后一列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个图形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') %图形标题注：同样这个程序使用的输出信号也没有噪音，所以辨识的结果没有误差，请同学们在输出信号中加入噪音，再使用RLS 对其辨识，观察辨识结果，进行分析。