2003版系统辨识最小二乘法大作业
系统辨识相关分析最小二乘
相关分析法辨识系统单位脉冲响应1辨识原理对于下图示的单输入单输出线性系统，其输入输出的因果关系可用卷积公式描述。
公式为：0()()()y t g x t d λλλ∞=-⎰把变量t 换成t +τ，上式两边同乘以x (t )，取时间的平均值，得11lim()(+)(){lim()(+)}22TTTTT T x t y t dt g x t x t dt d TTτλτλλ∞--→∞→∞=-⎰⎰⎰即 0()()()x y x R g R d τστλλ∞=-⎰上式即为维纳-霍夫方程，其给出了输入的自相关函数，输入、输出的互相关函数及脉冲响应函数三者之间的关系。
令x (t )为白噪声信号，则其自相关函数为：()(), ()()x x R k R k τδττλδτλ=-=-代入维纳-霍夫方程得：()()()()xy x R g R d kg τλτλλτ∞=-=⎰则有：()()xy R g kττ=这样，只要记录x(t)、y(t)的值，并计算它们的互相关函数，即可求得脉冲响应函数g(τ)。
在系统有正常输入的情形下，辨识脉冲响应的原理图如下图所示。
2辨识过程2.1预备实验以二阶系统22()2G s s s ++=作为辨识对象。
在实验前首先要进行预备实验，以了解系统特性。
通过简单阶跃响应确定系统过度过程时间T s 大约为11s ，如下图所示。
给系统施加不同周期的正弦信号，系统输出为输入的0.707倍时，确定截止频率f M 大约为0.318Hz 。
2.2选择二位式伪随机序列的参数（1）选择t 和N 由0.3Mt f ∆≤，得0.94t s ∆≤。
因为系统的时间常数1nT s ζω=，根据时间常数可按照0.050.1t T ∆= （）选择t ∆。
由二位式伪随机序列周期要大于系统过渡过程时间，若t ∆选择0.94s ，则由(1)s N t T -⨯∆≥，得12.7021N ≥；若t ∆选择0.195s ，则由(1)s N t T -⨯∆≥，得57.4103N ≥。
系统辨识之最小二乘法
系统辨识之最小二乘法方法一、最小二乘一次性算法：首先对最小二乘法的一次性辨识算法做简要介绍如下：过程的黑箱模型如图所示：其中u(k)和z(k)分别是过程的输入输出，)(1-z G 描述输入输出关系的模型，成为过程模型。
过程的输入输出关系可以描述成以下最小二乘格式：)()()(k n k h k z T +=θ （1）其中z(k)为系统输出，θ是待辨识的参数，h(k)是观测数据向量，n(k)是均值为0的随机噪声。
利用数据序列{z （k ）}和{h （k ）}极小化下列准则函数：∑=-=Lk T k h k z J 12])()([)(θθ （2）使J 最小的θ的估计值^θ，成为最小二乘估计值。
具体的对于时不变SISO 动态过程的数学模型为 )()()()()(11k n k u z B k z z A +=-- （3）应该利用过程的输入、输出数据确定)(1-z A 和)(1-Z B 的系数。
对于求解θ的估计值^θ，一般对模型的阶次a n ，b n 已定，且b a n n >；其次将（3）模型写成最小二乘格式)()()(k n k h k z T +=θ （4）式中=------=T n n T b a b a b b b a a a n k u k u n k z k z k h ],,,,,,,[)](,),1(),(,),1([)(2121 θ （5）L k ,,2,1 =因此结合式（4）（5）可以得到一个线性方程组L L L n H Z +=θ （6）其中==T L TL L n n n n L z z z z )](),2(),1([)](),2(),1([ （7）对此可以分析得出，L H 矩阵的行数为),max(b a n n L -，列数b a n n +。
在过程的输入为2n 阶次，噪声为方差为1，均值为0的随机序列，数据长度)(b a n n L +>的情况下，取加权矩阵L Λ为正定的单位矩阵I ，可以得出：L T L L T L z H H H 1^)(-=θ （8）其次，利用在Matlab 中编写M 文件，实现上述算法。
（完整）系统辨识—最小二乘法汇总-推荐文档
（完整）系统辨识—最⼩⼆乘法汇总-推荐⽂档(完整)系统辨识—最⼩⼆乘法汇总-推荐⽂档-CAL-FENGHAI-(2020YEAR-YICAI)_JINGBIAN最⼩⼆乘法参数辨识201403027摘要:系统辨识在⼯程中的应⽤⾮常⼴泛,系统辨识的⽅法有很多种,最⼩⼆乘法是⼀种应⽤极其⼴泛的系统辨识⽅法.阐述了动态系统模型的建⽴及其最⼩⼆乘法在系统辨识中的应⽤,并通过实例分析说明了最⼩⼆乘法应⽤于系统辨识中的重要意义.关键词:最⼩⼆乘法;系统辨识;动态系统Abstract: System identification in engineering is widely used, system identification methods there are many ways， least squares method is a very wide range of application of system identification method and the least squares method elaborated establish a dynamic system models in System Identification applications and examples analyzed by the least squares method is applied to illustrate the importance of system identification.Keywords: Least Squares; system identification; dynamic system引⾔随着科学技术的不断发展,⼈们认识⾃然、利⽤⾃然的能⼒越来越强,对于未知对象的探索也越来越深⼊.我们所研究的对象,可以依据对其了解的程度分为三种类型:⽩箱、灰箱和⿊箱.如果我们对于研究对象的内部结构、内部机制了解很深⼊的话,这样的研究对象通常称之为“⽩箱”;⽽有的研究对象,我们对于其内部结构、机制只了解⼀部分,对于其内部运⾏规律并不⼗分清楚,这样的研究对象通常称之为“灰箱”;如果我们对于研究对象的内部结构、内部机制及运⾏规律均⼀⽆所知的话,则把这样的研究对象称之为“⿊箱”.研究灰箱和⿊箱时,将研究的对象看作是⼀个系统,通过建⽴该系统的模型,对模型参数进⾏辨识来确定该系统的运⾏规律.对于动态系统辨识的⽅法有很多,但其中应⽤最⼴泛,辨识效果良好的就是最⼩⼆乘辨识⽅法,研究最⼩⼆乘法在系统辨识中的应⽤具有现实的、⼴泛的意义.1．1 系统辨识简介系统辨识是根据系统的输⼊输出时间函数来确定描述系统⾏为的数学模型。
系统辨识最小二乘法大作业 (2)
系统辨识大作业最小二乘法及其相关估值方法应用学院：自动化学院学号：姓名：日期：基于最小二乘法的多种系统辨识方法研究一、实验原理1.最小二乘法在系统辨识中用得最广泛的估计方法是最小二乘法(LS)。
设单输入-单输出线性定长系统的差分方程为(5.1.1)式中:为随机干扰；为理论上的输出值。
只有通过观测才能得到，在观测过程中往往附加有随机干扰。
的观测值可表示为(5.1.2)式中:为随机干扰。
由式(5.1.2)得(5.1.3)将式(5.1.3)带入式(5.1.1)得(5.1.4)我们可能不知道的统计特性，在这种情况下，往往把看做均值为0的白噪声。
设(5.1.5)则式(5.1.4)可写成(5.1.6)在观测时也有测量误差，系统内部也可能有噪声，应当考虑它们的影响。
因此假定不仅包含了的测量误差，而且包含了的测量误差和系统内部噪声。
假定是不相关随机序列(实际上是相关随机序列)。
现分别测出个随机输入值，则可写成个方程，即上述个方程可写成向量-矩阵形式(5.1.7) 设则式(5.1.7)可写为(5.1.8)式中：为维输出向量；为维噪声向量；为维参数向量；为测量矩阵。
因此式(5.1.8)是一个含有个未知参数，由个方程组成的联立方程组。
如果，方程数少于未知数数目，则方程组的解是不定的，不能唯一地确定参数向量。
如果，方程组正好与未知数数目相等，当噪声时，就能准确地解出(5.1.9)如果噪声，则(5.1.10)从上式可以看出噪声对参数估计是有影响的，为了尽量较小噪声对估值的影响。
在给定输出向量和测量矩阵的条件下求系统参数的估值，这就是系统辨识问题。
可用最小二乘法来求的估值，以下讨论最小二乘法估计。
2.最小二乘法估计算法设表示的最优估值，表示的最优估值，则有(5.1.11)写出式(5.1.11)的某一行，则有(5.1.12) 设表示与之差，即-(5.1.13)式中成为残差。
把分别代入式(5.1.13)可得残差。
设则有(5.1.14) 最小二乘估计要求残差的平方和为最小，即按照指数函数(5.1.15) 为最小来确定估值。
系统辨识大作业1201张青
《系统辨识》大作业学号：********班级：自动化1班姓名：**信息与控制工程学院自动化系2015-07-11第一题模仿index2，搭建对象，由相关分析法，获得脉冲响应序列ˆ()g k，由ˆ()g k，参照讲义，获得系统的脉冲传递函数()G z和传递函数()G s;应用最小二乘辨识,获得脉冲响应序列ˆ()g k;同图显示两种方法的辨识效果图；应用相关最小二乘法，拟合对象的差分方程模型；构建时变对象，用最小二乘法和带遗忘因子的最小二乘法，(可以用辨识工具箱) 辨识模型的参数，比较两种方法的辨识效果差异;答：根据index2搭建结构框图：相关分析法：利用结构框图得到UY 和tout其中的U就是题目中要求得出的M序列，根据结构框图可知序列的周期是1512124=-=-=npN。
在command window中输入下列指令，既可以得到脉冲相应序列()g k：aa=5;NNPP=15;ts=2; RR=ones(15)+eye(15); for i=15:-1:1UU(16-i,:)=UY(16+i:30+i,1)'; endYY=[UY(31:45,2)];GG=RR*UU*YY/[aa*aa*(NNPP+1)*ts]; plot(0:2:29,GG) hold onstem(0:2:29,GG,'filled') Grid;title('脉冲序列g(τ)')最小二乘法建模的响应序列由于是二阶水箱系统，可以假设系统的传递函数为221101)(sa s a sb b s G +++=，已知)(τg ，求2110,,,a a b b已知G （s ）的结构，用长除法求得G(s)的s 展开式，其系数等于从 )( g 求得的各阶矩，然后求G(s)的参数。
得到结果： a1 =-1.1561 a2 =0.4283 b0 =-0.0028 b1=0.2961在command window 中输入下列指令得到传递函数：最小二乘一次算法相关参数%最小二乘法一次完成算法 M=UY(:,1); z=UY(:,2); H=zeros(100,4); for i=1:100 H(i,1)=-z(i+1); H(i,2)=-z(i); H(i,3)=M(i+1); H(i,4)=M(i); endEstimate=inv(H'*H)*H'*(z(3:102)) %结束得到相关系数为：Estimate =-0.7866 0.1388 0.5707 0.3115带遗忘因子最小二乘法：%带遗忘因子最小二乘法程序M=UY(:,1);z=UY(:,2);P=1000*eye(5); %Theta=zeros(5,200); %Theta(:,1)=[0;0;0;0;0];K=zeros(4,400); %K=[10;10;10;10;10];lamda=0.99;%遗忘因数for i=3:201h=[-z(i-1);-z(i-2);M(i);M(i-1);M(i-2)];K=P*h*inv(h'*P*h+lamda);Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h'*Theta(:,i-2));P=(eye(5)-K*h')*P/lamda;endi=1:200;figure(1)plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:),i,Theta(5,:) )title('带遗忘因子最小二乘法')grid%结束Estimate 可由仿真图得出，可知两种方法参数确定十分接近。
系统辨识与参数估计大作业
系统辨识与参数估计大作业第一题 递推最小二乘估计参数考虑如上图所示的仿真对象，选择模型结构为：)()2()1()2()1()(2121k v k u b k u b k z a k z a k z +-+-=-+-+，其中)(k v 是服从)1,0(N 正态分布的不相关随机噪声；输入信号)(k u 采用4阶逆M 序列，特征多项式取41)(s s s F ⊕⊕=，幅度为1，循环周期为bit N p 62=；控制λ值，使数据的噪信比分别为10％，73％，100％三种情况。
加权因子1)(=Λk ；数据长度L=500;初始条件取001.0)0(ˆ,10)0(6==θI P ， （1） 利用递推最小二乘算法在线估计参数，（2） 利用模型阶次辨识方法（AIC 准则），确定模型的阶次。
（3） 估计噪声)(k v 的方差和模型静态增益K （4） 作出参数估计值随时间的变化图 答：设过程的输入输出关系可以描述成()()()Tz k h k n k θ=+()z k 是输出量，()h k 是可观测的数据向量，n (k)是均值为0的随机噪声[]()(1),(2),(1),(2)Th k z k z k u k u k =------[]1212,,,Ta ab b θ=选取的模型为结构是1212()(1)(2)(1)(2)z k a z k a z k bu k bu k =----+-+- 12121.5,0.7, 1.0,0.5a a b b =-===加权最小二乘参数估计递推算法RWLS 的公式如下，11()(1)()()(1)()()()(1)()()()(1)()()()(1)T T TK k p k h k h k p k h k k k k K k z k h k k p k I K k h k p k θθθ-⎡⎤=--+⎢⎥Λ⎣⎦⎡⎤=-+--⎣⎦⎡⎤=--⎣⎦为了把p(k)的对称性，可以把p(k)写成1()(1)()()()(1)()()T T p k p k K k K k h k p k h k k ⎡⎤=---+⎢⎥Λ⎣⎦如果把()k Λ设成1的时候，加权最小二乘法就退化成最小二乘法。
最小二乘法系统辨识
最小二乘法一次完成算法的MATLAB仿真一、二阶系统的最小二乘法一次完成算法的辨识程序及结果程序源代码：clear % 工作间清零u = [-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1]; % 系统辨识的输入信号为一个周期的M序列z = zeros( 1, 16); % 定义输出观测值的长度for k = 3 : 16% 用理想输出值作为观测值z(k) = 1.8*z(k-1) - 0.8*z(k-2) + 1.1*u(k-1) + 0.4*u(k-2);ZL( k, : ) = z(k); % 给样本矩阵ZL赋值Hl = [ -z(k-1), -z(k-2), u(k-1), u(k-2) ]; % 用理想输出值作为观测值HL( k, : ) = Hl;endsubplot(3,1,1) % 画三行一列图形窗口中的第一个图形stem(u) % 画出输入信号u的经线图形subplot(3,1,2) % 画三行一列图形窗口中的第二个图形i=1:1:16; % 横坐标范围是1到16，步长为1% 图形的横坐标是采样时刻i, 纵坐标是输出观测值z, 图形格式为连续曲线plot(i,z)subplot(3,1,3) % 画三行一列图形窗口中的第三个图形stem(z),grid on % 画出输出观测值z的经线图形，并显示坐标网格u,z % 显示输入信号和输出观测信号c1=HL'*HL; % 计算参数并显示c2=inv(c1);c3=HL'*ZL;c=c2*c3a1=c(1), a2=c(2), b1=c(3), b2=c(4) % 从中分离出并显示a1 、a2、b1、b2运算结果如下：u =-1 1 -1 1 1 1 1 -1 -1 -1 1 -1 -1 1 1z =Columns 1 through 90 0 0.7000 0.5600 1.1480 3.1184 6.1947 10.1558 12.6246 Columns 10 through 1613.0997 11.9798 11.7838 10.9270 8.7416 7.6933 8.3546c =-1.80000.80001.10000.4000a1 =-1.8000a2 =0.8000b1 =1.1000b2 =0.4000所画图形如下：二、实际压力系统的最小二乘辨识程序及结果程序源代码：clear % 工作间清零V = [ 54.3, 61.8, 72.4, 88.7, 118.6, 194.0]' % 赋初值V，并显示P = [ 61.2, 49.5, 37.6, 28.4, 19.2, 10.1]' % 赋初值P，并显示% P、V之间的关系:logP=-alpha*logV% +logbeita=[-logV,1][alpha,log(beita)]'=HL*sitafor i = 1 : 6; % 循环变量的取值为从1到6 Z(i) = log( P (i) ); % 赋系统的输出采样值ZL( i, : ) = Z(i); % 给样本矩阵ZL赋值Hl = [ log( V(i) ),1 ]; % 给HL赋值HL( i, : ) = Hlend % 循环结束c1=HL'*HL; % 计算参数并显示c2=inv(c1);c3=HL'*ZL;c4=c2*c3alpha = c4(1) % 为c4的第1个元素beita = exp(c4(2)) % 为以自然数为底的c4的第2个元素的指数运算结果如下：V =54.300061.800072.400088.7000118.6000194.0000P =61.200049.500037.600028.400019.200010.1000HL =3.9945 1.0000HL =3.9945 1.00004.1239 1.0000HL =3.9945 1.00004.1239 1.00004.2822 1.0000HL =3.9945 1.00004.1239 1.00004.2822 1.00004.4853 1.0000HL =3.9945 1.00004.1239 1.00004.2822 1.00004.4853 1.00004.7758 1.0000HL =3.9945 1.00004.1239 1.00004.2822 1.00004.4853 1.00004.7758 1.00005.2679 1.0000 c4 =-1.40429.6786alpha =-1.4042beita =1.5972e+004。
系统辨识作业及答案
一． 问答题1. 介绍系统辨识的步骤。
答：（1）先验知识和建模目的的依据；（2）实验设计；（3）结构辨识；（4）参数估计；（5）模型适用性检验。
2. 考虑单输入单输出随机系统，状态空间模型[])()(11)()(11)(0201)1(k v k x k y k u k x k x +=⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡=+ 转换成ARMA 模型。
答：ARMA 模型的特点是u(k)=0,[])()(11)()(0201)1(k v k x k y k x k x +=⎥⎦⎤⎢⎣⎡=+3. 设有一个五级移位寄存器，反馈取自第2级和第3级输出的模2加法和。
试说明：（1） 其输出序列是什么？ （2） 是否是M 序列？（3） 它与反馈取自第4级与第3级输出模2加法和所得的序列有何不同？ （4） 其逆M 序列是什么？ 答：（1）设设输入序列1 1 1 1 1111018110107101006010015100114001113011112111111）（）（）（）（）（）（）（）（()()()()()()()01110161110115110101410100)13(010011210011110011110011109()()()()()()()001112401110)23(111012211010211010020010011910011180011117()()()()()()()()10011320011131011103000111291101028101002701001261001125 其输出序列为：1 1 1 1 1 0 0 1 0 1⑵不是M 序列⑶第4级与第3级模2相加结果100108001007010006100015000114001113011112111111）（）（）（）（）（）（）（）（()()()()()()()11110161110115110101410101)13(010111210110110110010110019()()()()()()()110012410010)23(001002201000211000120000111900111180111117()()()()()()()()01111321111031111013011010291010128010112710110260110025 不同点：第2级和第3级模二相加产生的序列，是从第4时刻开始，每隔7个时刻重复一次；第4级与第3级模2相加产生的,序列，是从第2时刻开始每隔15个时刻重复一次。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
西北工业大学系统辩识大作业题目：最小二乘法系统辨识一、 问题重述：用递推最小二乘法、加权最小二乘法、遗忘因子法、增广最小二乘法、广义最小二乘法、辅助变量法辨识如下模型的参数离散化有z^4 - 3.935 z^3 + 5.806 z^2 - 3.807 z + 0.9362---------------------------------------------- =z^4 - 3.808 z^3 + 5.434 z^2 - 3.445 z + 0.8187噪声的成形滤波器离散化有4.004e-010 z^3 + 4.232e-009 z^2 + 4.066e-009 z + 3.551e-010----------------------------------------------------------------------------- =z^4 - 3.808 z^3 + 5.434 z^2 - 3.445 z + 0.8187采样时间0.01s要求：1.用Matlab 写出程序代码；2.画出实际模型和辨识得到模型的误差曲线；3.画出递推算法迭代时各辨识参数的变化曲线；最小二乘法:在系统辨识领域中 ,最小二乘法是一种得到广泛应用的估计方法 ,可用于动态 ,静态 , 线性 ,非线性系统。
在使用最小二乘法进行参数估计时 ,为了实现实时控制 ,必须优化成参数递推算法 ,即最小二乘递推算法。
这种辨识方法主要用于在线辨识。
MATLAB 是一套高性能数字计算和可视化软件 ,它集成概念设计 ,算法开发 ,建模仿真 ,实时实现于一体 ,构成了一个使用方便、界面友好的用户环境 ,其强大的扩展功能为各领域的应用提供了基础。
对4324326.51411.5320120232320Y s s s s G U s s s s ++++==++++432120120232320E N W s s s s ==++++于一个简单的系统,可以通过分析其过程的运动规律 ,应用一些已知的定理和原理,建立数学模型 ,即所谓的“白箱建模 ”。
但对于比较复杂的生产过程 ,该建模方法有很大的局限性。
由于过程的输入输出信号一般总是可以测量的 ,而且过程的动态特性必然表现在这些输入输出数据中 ,那么就可以利用输入输出数据所提供的信息来建立过程的数学模型。
这种建模方法就称为系统辨识。
把辨识建模称作“黑箱建模”。
系统辨识又分为参数辨识和阶次辨识 ,在本文中只讨论参数辨识问题最小二乘递推算法所用的模型：Z(k)=B()u(k)+v(k)最小二乘递推算法为:)(k v 是服从N )1,0(分布的不相关随机噪声。
)(1-z G )()(11--=z A z B ，)(1-z N )()(11--=z C z D ，考虑如图 1示仿真对象 ,系统的差分方程为z(k)=3.808*z(k-1)-5.434*z(k-2)+3.445*z(k-3)-0.8187*z(k-4)+u(k)-3.935*u(k-1)+5.806*u(k-2)-3.807*u(k-3)+0.9362*u(k-4)+)(k v (3.69) 仿真对象选择如下的模型结构)()4(5)3()2(3)1()()4(4)3(3)2()1()(42121k v k u b k u b k u b k u b k u b k z a k z a k z a k z a k z +-+-+-+-+=-+-+-+-+(3.68)+ e (k ) 图1 最小二乘递推算法辨识实例结构图y (k )u (k )z (k )v (k ))(1-z N)(1-z G其中，)(k v 是服从正态分布的白噪声N )1,0(。
输入信号采用4位移位寄存器产生的M 序列，幅度为1。
按式(3.69)构造h (k )；加权阵取单位阵I Λ L ；利用式(3.61)计算K (k )、)(ˆk θ和P (k )，计算各次参数辨识的相对误差，精度满足要求式(3.67)后停机。
最小二乘递推算法辨识的Malab6.0程序流程图：最小二乘递推算法辨识程序clear %清理工作间变量L=35; % M序列的周期y1=1;y2=1;y3=1;y4=0; %四个移位寄存器的输出初始值for i=1:L;%开始循环，长度为Lx1=xor(y3,y4); %第一个移位寄存器的输入是第三个与第四个移位寄存器的输出的“或”x2=y1; %第二个移位寄存器的输入是第一个移位寄存器的输出x3=y2; %第三个移位寄存器的输入是第二个移位寄存器的输出x4=y3; %第四个移位寄存器的输入是第三个移位寄存器的输出y(i)=y4; %取出第四个移位寄存器的幅值为"0"和"1"的输出信号，即M序列if y(i)>0.5,u(i)=-1; %如果M序列的值为"1", 辨识的输入信号取“-1”else u(i)=1; %如果M序列的值为"0", 辨识的输入信号取“1”end %小循环结束y1=x1;y2=x2;y3=x3;y4=x4; %为下一次的输入信号做准备end %大循环结束，产生输入信号ufigure(1); %第一个图形stem(u),grid on %显示出输入信号径线图并给图形加上网格z(2)=0;z(1)=0; z(3)=0;z(4)=0;%设z的前四个初始值为零for k=5:25; %循环变量从5到25z(k)=3.808*z(k-1)-5.434*z(k-2)+3.445*z(k-3)-0.8187*z(k-4)+u(k)-3.935*u(k-1)+5.8 06*u(k-2)-3.807*u(k-3)+0.9362*u(k-4); %输出采样信号end%RLS递推最小二乘辨识c0=[0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001]'; %直接给出被辨识参数的初始值,即一个充分小的实向量p0=10^6*eye(9,9); %直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.000000005; %取相对误差E=0.000000005c=[c0,zeros(9,24)]; %被辨识参数矩阵的初始值及大小e=zeros(9,25); %相对误差的初始值及大小for k=5:25; %开始求Kh1=[-z(k-1),-z(k-2),-z(k-3),-z(k-4),u(k),u(k-1),u(k-2),u(k-3),u(k-4)]';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,:);a3=c(3,:); a4=c(4,:);b1=c(5,:); b2=c(6,:); b3=c(7,:); b4=c(8,:); b5=c(9,:);ea1=e(1,:); ea2=e(2,:);ea3=e(3,:); ea4=e(4,:);eb1=e(5,:); eb2=e(6,:); eb3=e(7,:); eb4=e(8,:); eb5=e(9,:);figure(2); %第二个图形i=1:25; %横坐标从1到25plot(i,a1,'r',i,a2,':',i,a3,'r',i,a4,'k',i,b1,'y',i,b2,':',i,b3,'m',i,b4,':',i, b5,'g') %画出a1，a2，a3,a4,b1，b2,b3,b4,b5的各次辨识结果title('参数变化曲线') %图形标题figure(3); %第三个图形i=1:25; %横坐标从1到25plot(i,ea1,'r',i,ea2,'k:',i,ea3,'y',i,ea4,'m',i,eb1,'g',i,eb2,'c:',i,eb3,'r',i, eb4,':',i,eb5,'g') %画出a1，a2，a3,a4,b1，b2,b3,b4,b5的各次辨识结果的收敛情况title('误差曲线') %图形标题参数变化图相对误差图仿真结果表明，大约递推到第十五步时，参数辨识的结果基本达到稳定状态，即a1=-3.808, a2=5.434, a3=-3.445, a4= -0.8187, b1=1, b2=-3.935 b3=5.806 b4=-3.807 b5=0.9362。
此时，参数的相对变化量E 0.000000005。
从整个辨识过程来看，精度的要求直接影响辨识的速度。
虽然最终的精度可以达到很小，但开始阶段的相对误差会很大，从相对误差图形中可看出，参数的最大相对误差会达到三位数增广最小二乘法算法所用的模型：Z(k)=B()u(k)+D()e(k)增广最小二乘法算法为：模型结构选用如下形式：)4)-v(k *4d 3)-v(k *3d 2)-v(k *2d 1)-v(k *1d )4(5)3()2(3)1()()4(4)3(3)2()1()(42121++++-+-+-+-+=-+-+-+-+k u b k u b k u b k u b k u b k z a k z a k z a k z a k z增广最小二乘算法流程图如下页图所示增广最小二乘辨识程序clearL=55; %4位移位寄存器产生的M序列的周期y1=1;y2=1;y3=1;y4=0; %4个移位寄存器的输出初始值for i=1:L;x1=xor(y3,y4); %第一个移位寄存器的输入信号x2=y1; %第二个移位寄存器的输入信号x3=y2; %第三个移位寄存器的输入信号x4=y3; %第四个移位寄存器的输入信号y(i)=y4; %第四个移位寄存器的输出信号，M序列, 幅值"0"和"1"，if y(i)>0.5,u(i)=-1; %M序列的值为"1"时,辨识的输入信号取“-1”else u(i)=1; %M序列的值为"0"时,辨识的输入信号取“1”endy1=x1;y2=x2;y3=x3;y4=x4; %为下一次的输入信号准备endfigure(1); %画第一个图形subplot(2,1,1); %画第一个图形的第一个子图stem(u),grid on %画出M序列输入信号v=randn(1,60); %产生一组60个正态分布的随机噪声subplot(2,1,2); %画第一个图形的第二个子图plot(v),grid on; %画出随机噪声信号u ,v %显示输入信号和噪声信号z=zeros(13,55); %输出采样矩阵z(2)=0; z(1)=0;z(3)=0; z(4)=0;%输出采样的初值c0=[0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001]'; %直接给出被辨识参数的初始值,即一个充分小的实向量p0=10^6*eye(13,13); %直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.00000000005; %相对误差E=0.000000005c=[c0,zeros(13,54)]; %被辨识参数矩阵的初始值及大小e=zeros(13,55); %相对误差的初始值及大小for k=5:55; %开始求Kz(k)=3.808*z(k-1)-5.434*z(k-2)+3.445*z(k-3)-0.8187*z(k-4)+u(k)-3.935*u(k-1)+5.8 06*u(k-2)-3.807*u(k-3)+0.9362*u(k-4)+4.004e-010*v(k-1)+4.232e-009*v(k-2)+4.066e -009*v(k-3)+3.551e-010*v(k-4); %系统在M序列输入下的输出采样信号h1=[-z(k-1),-z(k-2),-z(k-3),-z(k-4),u(k),u(k-1),u(k-2),u(k-3),u(k-4),v(k-1),v(k -2),v(k-3),v(k-4)]'; %为求K(k)作准备x=h1'*p0*h1+1; x1=inv(x); k1=p0*h1*x1; %Kd1=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]; %find p(k)p0=p1; %给下次用if e2<=E break; %若收敛情况满足要求，终止计算end %判断结束end %循环结束c, e, %显示被辨识参数及参数收敛情况%分离变量a1=c(1,:); a2=c(2,:);a3=c(3,:); a4=c(4,:);b1=c(5,:); b2=c(6,:); b3=c(7,:); b4=c(8,:); b5=c(9,:); %分离出a1，a2，a3,a4,b1，b2,b3,b4,b5d1=c(10,:); d2=c(11,:); d3=c(12,:); d4=c(13,:);%分离出d1、 d2、 d3 d4ea1=e(1,:); ea2=e(2,:);ea3=e(3,:); ea4=e(4,:);eb1=e(5,:); eb2=e(6,:); eb3=e(7,:); eb4=e(8,:); eb5=e(9,:); %分离出a1，a2，a3,a4,b1，b2,b3,b4,b5的收敛情况ed1=e(10,:); ed2=e(11,:); ed3=e(12,:);ed4=e(13,:); %分离出d1 、d2 、d3 d4的收敛情况figure(2); %画第二个图形i=1:55;plot(i,a1,'r',i,a2,':',i,a3,'r',i,a4,'k',i,b1,'y',i,b2,':',i,b3,'m',i,b4 ,':',i,b5,'g',i,d1,'k',i,d2,'g:',i,d3,'m',i,d4,'r') %画出各个被辨识参数title('参数变化曲线') %标题figure(3); i=1:55; %画出第三个图形plot(i,ea1,'r',i,ea2,'k:',i,ea3,'y',i,ea4,'m',i,eb1,'g',i,eb2,'c:',i,eb3,'r',i, eb4,':',i,eb5,'g',i,ed1,'y',i,ed2,'g:',i,ed3,'k',i,ed4,'m') %画出各个参数收敛情况title('误差曲线') %标题参数变化曲线相对误差曲线仿真结果表明，递推到第十步时，参数辨识的结果基本达到稳定状态，即a 1=-3.808, a 2=5.434,a 3=-3.445, a 4= -0.8187,b 1=1, b 2=-3.935 b 3=5.806 b 4=-3.807 b 5=0.9362，d 1=0.0000, d 2=0.0000, d 3=0.0000，d4=0.0000。