系统辨识与自适应控制--大作业

1 辨识的对象模型假设有一理想数学模型，它的离散化方程如下式所示：() 1.8(1)0.3(2) 1.2(1)(2)()y k y k y k u k u k e k +-+-=-+-+式中，()e k 是服从正态分布的白噪声)1,0(N ，()k u 为系统输入，()k y 为系统输出。

现在输入信号采用4阶M 序列，其幅值为1。

假设系统的模型阶次是已知的，即1212()(1)(2)(1)(2)()y k a y k a y k bu k b u k e k +-+-=-+-+。

下面采用递推最小二乘参数辨识。

2 递推最小二乘参数辨识方法简单的最小二乘参数辨识一次性方法计算复杂，不能够进行在线辨识，而且所需要的计算存储空间很大，而很多计算都是重复的计算。

为了解决这个问题，并实现在线的实时辨识，引入递推的最小二乘参数辨识。

递推最小二乘参数辨识的整体思想是，最新辨识出来的参数是建立在上次辨识的参数基础上，根据最新得到的辨识数据，对辨识的参数添加了一个参数增量。

下面利用数学语言对递推最小二乘参数辨识方法进行描述。

根据最小二乘原理，用N 次观测数据，得出参数向量θ的最小二乘估计lθˆ1()()T T N N NH H H Y N θ-= (1) 其中，ˆNθ表示根据N 次观测数据所得到的最小二乘值计量，下表N 表示该符号代表N 次观测数据构成的矩阵。

()[(1),(2),...,()]T Y N y y y N = (2)N H =(0).....(1)(0).....(1)(1).....(2)(1).....(2)..(1).....()(1).....()y y n u u n y y n u u n y N y N n u N u N n ----⎡⎤⎢⎥----⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥------⎣⎦（3）然后令1()TN N N P H H -=，且N P 是一方阵，它的维数取决于未知数的个数，而与观测次数无关。

则11111TN N N N P P h h --+++⎡⎤=+⎣⎦ （4） 式中1N h +表示第1N +次观测数据。

利用矩阵反演公式计算（4）式1111111111T T TN N N N N N N N N N N N P P h h P P h I h P h h P ---+++++++⎡⎤⎡⎤=+=-+⎣⎦⎣⎦ （5） 式中11TN N N I h P h ++⎡⎤+⎣⎦是一个标量，其求逆只是一个简单的除法。

令1111TN N N N I h P h γ-+++⎡⎤=+⎣⎦，于是有1111T N N N N N N P I P h h P γ++++⎡⎤=-⎣⎦，则最小二乘估计量ˆN θ的表达式可改写为1ˆ()T N N N P h Y N θ+=，于是第1N +次观测数据所得到的最小二乘值计量为111111ˆˆ(1)(1)T T N N N N N N P h Y N P P h Y N θθ++++++⎡⎤=+=++⎣⎦由（4）式可知111T T TN N N N P P h h +++=- （6）于是111ˆˆˆ(1)T N N N N N K Y N h θθθ+++⎡⎤=++-⎣⎦ （7）111N N N K P h +++= （8）将上述的推导过程总结如下：1111111111ˆˆˆy(1)1T N N N N N N N N TN N N T N N N N K N h P h K h P h P I K h Pθθθ++++++++++⎧⎡⎤=++-⎣⎦⎪⎪=⎨+⎪⎪⎡⎤=-⎣⎦⎩ 3 辨识过程整体上的辨识过程流程图如下图所示：图-1 递推最小二乘参数辨识流程示意图详细步骤如下所示：1 在此选择4阶的M 序列。

这样，M 序列的最大循环周期长度为15（bit ）。

首先利用MATLAB 程序生成4阶幅值为1的M 序列和服从正态分布的白噪声序列。

2 根据上述的输入序列，带入待辨识的系统模型中，得到一系列对应的输出信号序列。

然后根据输入和输出序列，组建辨识所需的观测矩阵。

本系统的观测矩阵形式如下：(2)(1)(2)(1)(3)(2)(3)(2)(14)(13)(14)(13)m y y u u y y u u H y y u u --⎡⎤⎢⎥--⎢⎥=⎢⎥⎢⎥--⎣⎦3选取待辨识的参数向量为模极小的数值，选取迭代矩阵P 的初值为极大的矩阵，如20,100.P I σσ==4 计算增益矩阵G(1)()()1()(1)()TP k k G k k P k k φφφ-=+- 其中T ()[(1) ... () (1) ... u() ]q b k y k y k n u k k n φ=------5 将上述的初值和增益矩阵G(k)代入上述的递推最小二乘参数辨识方程组中进行循环计算。

直到计算的结果精度符合要求或者迭代的次数达到最大为止。

4 辨识结果与分析利用matlab 对上述系统进行辨识方针，输入的M 序列如下图所示：图-2 M序列信号图辨识的结果如下图所示：图-3 递推最小二乘参数辨识结果图辨识结果与真实的系统参数对比：对结果的分析：1 利用递推最小二乘参数辨识方法，能够很好的辨识出系统的参数，并且最后收敛稳定。

其实，对于稳定的系统，递推最小二乘算法随着递推次数的增加，所得到的辨识结果将会越来越接近真实值。

但是对于其他系统，则会存在数据饱和的现象，随着递推次数的增加，辨识结果不但不接近真实数值，反而出现发散的情况。

对于这种情况，需要在递推最小二乘参数辨识法的基础上进一步进行修改，得到限定记忆法和加权最小二乘参数辨识法。

2 本次仿真所采用的M序列阶次很小，明显导致了参数仿真得到的递推估计曲线在开始的5秒钟内变化非常剧烈，而且非常不光滑。

5 仿真源程序%利用递推最小二乘法进行参数辨识%%%%%%%%%%产生M序列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)=-1;%如果M序列的值为"1"时,辨识的输入信号取“-0.03”else u(i)=1;%当M序列的值为"0"时,辨识的输入信号取“0.03”end %小循环结束y1=x1;y2=x2;y3=x3;y4=x4; %为下一次的输入信号做准备end %大循环结束，产生输入信号ufigure(1); stem(u);grid on;title('输入信号')%%%%%%%%%%% 产生白噪声A=6;x0=1;M=255;for k=1:10000x2=A*x0;x1=mod (x2,M);v1=x1/256;v(:,k)=(v1-0.5)*2;x0=x1;v0=v1；endnum=v;k1=k;%%%%%%%%%% 递推最小二乘辨识程序lamt=1;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)+0*num(1,k);%给出理想辨识输出采样信号endc0=[0.001 0.001 0.001 0.001]';%直接给出被辨识参数的初始值,即一个充分小的实向量p0=10^3*eye(4,4);%直接给出协方差阵初始状态P0，即一个充分大的实数单位矩阵E=0.000000005; %相对误差E=0.000000005c=[c0,zeros(4,14)]; %被辨识参数矩阵的初始值及大小e=zeros(4,15); %相对误差的初始值及大小for k=3:15;h1=[-z(k-1),-z(k-2),u(k-1),u(k-2)]';x=h1'*p0*h1+1*lamt;x1=inv(x); %开始求K(k)k1=p0*h1*x1; %求出K的值d1=z(k)-h1'*c0;c1=c0+k1*d1; %求被辨识参数cp1=1/lamt*(eye(4)-k1*h1')*p0;e1=c1-c0; %求参数当前值与上一次的值的差值e2=e1./c0; %求参数的相对变化e(:,k)=e2; %当前相对变化列向量加入误差矩阵的最后一列c(:,k)=c1; %把辨识参数c列向量加入辨识参数矩阵c0=c1; %新获得的参数作为下一次递推的旧参数p0=p1; %给下次用if e2<=Ebreak; %若参数收敛满足要求，终止计算endend%% %%%%%%%%%得到辨识结果，并作图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);i=1:15;plot(i,a1,'r:',i,a2,'b-.',i,b1,'k--',i,b2,'g','LineWidth',2)legend('a1','a2','b1','b2');title('递推最小二乘参数辨识')（在得到辨识结果之后，对辨识结果图进行坐标注视和图形的修改。

记得到报告中的结果图。

）。