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个时刻重复一次。
系统辩识作业题
系统辨识大作业
一.设SlSO系统差分方程为
y(k)=—α1y(k-1)-a2y(k-2)+bλu(k-1)+b2u(k-2)+ξ{k)
辨识参数向量为θ=[q a2b l b2]r,输入输出数据详见数据文件UyLtXt—uy3.txtoξ(k)为噪声方差各异的白噪声或有色噪声。
试求解:
1)用n元一次方程解析法,再求其平均值方法估计。
2)用最小二乘及递推最小二乘法估计。
;
3)用辅助变量法及其递推算法估计
4)用广义最小二乘法及其递推算法估计
5)用夏氏偏差修正法、夏氏改良法及其递推算法估计
6)用增广矩阵法估计
7)分析噪声父攵)特性;
二.用极大似然法估计6。
三.以上题的结果为例,进行:
1.分析比较各种方法估计的精度;
2.分析其计算量;
3.分析噪声方差的影响;
4.比较白噪声和有色噪声对辨识的影响。
四.系统模型阶次的辨识:
1.用三种方法确定系统的阶次并辨识;
2.分析噪声对定阶的影响;
3.比较所用三种方法的优劣及有效性;
五.给出由正弦输入求取系统开环频率响应特性曲线的辨识方法。
六.提出一种自己创造的辨识新方法,并用所给数据进行辨识验证。
注:闭卷考试时提交大作业报告。
系统辨识作业_梯度法最小二乘法
系统辨识一、最小二乘算法仿真对象为:)()2(5.0)1()2(7.0)1(5.1)(k v k k u k k z k z +-+-=-+-+,)(k v 为服从N(0,1)分布的白噪声,输入信号采用4阶M 序列。
辨识模型为:)()2()1()2()1()(2121k v k z b k z b k z a k z a k z +-+-+----=,观测数据长度取400=L ,加权矩阵取I =Λ。
则matlab 程序如下:% M array constructionNp=15;r=4;X1=1;X2=1;X3=1;X4=1;m_length = r*Np;a=1;for i=1:1:m_lengthY4=X4;Y3=X3;Y2=X2;Y1=X1;X4=Y3;X3=Y2;X2=Y1;X1=xor(Y3,Y4);if Y4==0M(i)=-a;elseM(i)=a;endendfigure;i=1:1:m_length;plot(i,M);% 白噪声noise = zeros(1,m_length);for i=1:1:m_lengthtemp = noise + 0.5*rands(1,m_length);noise = temp;endnoise = noise/12;%noise = temp;% parameter of systemn=2;d=1;a1=-1;a2=0.5;b1=1;b2=0.5;S_U0=0.2;S_Y0=0.2;% generate u,yu0=ones(1,m_length)*S_U0;U = M + u0 + noise;figure;i=1:1:m_length;plot(i,U);%formulationy(1)=0;y(2)=0;y(3)=0;Y(1)=S_Y0+y(1)+noise(1);Y(2)=S_Y0+y(2)+noise(2);Y(3)=S_Y0+y(3) +noise(3); for k=4:m_lengthy(k) = b1*U(k-1-d)+b2*U(k-2-d)-a1*y(k-1)-a2*y(k-2);Y(k)=S_Y0+y(k)+noise(k);endfigure;i=1:1:m_length;plot(i,Y);%辨识% premitive valuec=2;P = (c^3)*eye(3*n+d);sita(:,3) = [0,0,0,0,0,0,0]';alf = 0.995;% M%compute U0,Y0sum_U = 0;sum_Y = 0;for k=1:1:m_lengthsum_U = sum_U + U(k);sum_Y = sum_Y + Y(k);endt_U0 = sum_U/m_length;t_Y0 = sum_Y/m_length;for k=1:1:m_lengtht_u(k) = U(k) - t_U0;t_y(k) = Y(k) - t_Y0;endfigure;i=1:1:m_length;plot(i,t_u);figure;i=1:1:m_length;plot(i,t_y);v(1)=0;v(2)=0;v(3)=0;for k=4:1:m_lengthfai = [-t_y(k-1),-t_y(k-2),t_u(k-1-d),t_u(k-2-d),v(k-1),v(k-2),v(k-3)]';v(k) = t_y(k) - fai'*sita(:,k-1);W = P*fai;dan = 1/(alf + fai'*W);sita(:,k) = sita(:,k-1) + dan*W*v(k);P = ( P - dan * W*W')/alf;end仿真结果如下图所示:则各个参数的估计值为:⎪⎪⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛∧∧∧∧5048.00062.17003.04961.1121bbaa二、梯度校正辨识方法仿真对象为)()2(5.0)1()2(7.0)1(5.1)(kvkkukkzkz+-+-=-+--,其中)(kv为服从N(0,1)分布的白噪声,取初值0)0(=∧θ,1=α,1.0=c。
系统辨识大作业
系统辨识大作业专业班级:自动化09-3学号:09051325姓名:吴恩作业一:设某物理量Y与X满足关系式2=++,实验获得一批数据Y aX bX c如下表,试辨识模型参数,,a b c。
(15分)解答:问题描述:由题意知,这是一个已知模型为Y=aX2+bX+c,给出了10组实验输入输出数据,要求对模型参数a,b,c进行辨识。
问题求解:这里对该模型参数辨识采用最小二乘法的一次算法(LS)求解。
2=++可以写成矩阵形式Y=AE+e;其中A=[X^2,X,1]构成, Y aX bX c利用matlab不难求解出结果。
运行结果:利用所求的的参数,求出给定的X对应的YE值,列表如下做出上表的图形如下12345678910xyy=ax 2+bx+c 参数求解结果分析:根据运行结果可以看出,拟合的曲线与真是观测的数据有误差,有出入,但是误差较小,可以接受。
出现误差的原因,一方面是由于给出的数据只有十个点,数据量太少,难以真正的充分的计算出其参数,另外,该问题求解采用的是LS 一次算法,因此计算方法本身也会造成相应的误差。
作业二:模仿实验二,搭建对象,由相关分析法,获得脉冲相应序列()g k,由()G z;和传递函数g k,参照讲义,获得系统的脉冲传递函数()G s及应用相关最小二乘法,拟合对象的差分方程模型;加阶跃()扰动,用最小二乘法和带遗忘因子的最小二乘法,辨识二阶差分方程的参数,比较两种方法的辨识差异;采用不少于两种定阶方法,确定对象的阶次。
对象模型如图:利用相关分析法,得到对象的脉冲相应序列。
如下图:(1).由脉冲相应序列,求解系统的脉冲传递函数G(z)Transfer function:0.006072 z^2 + 0.288 z + 0.1671-------------------------------z^2 + 0.1018 z - 0.7509Sampling time: 2(2).由脉冲相应序列求解系统的传递函数G(s)Transfer function:(0.04849+2.494e-018i)-----------------------s^2 + 0.1315 s + 0.6048(3).利用相关最小二乘法拟合系统的差分方程模型如下:(4).在t=100,加入一个0.5的阶跃扰动,,利用RLS求解差分方程模型:RLS加入遗忘因子之后与未加之前的曲线情况如下:未加遗忘因子之前参数以及残差的计算过程加入0.99的遗忘因子得到的参数辨识过程与残差的变化过程根据上面两种方法所得到的误差曲线和参数过渡过程曲线,我们可以看出来利用最小二乘法得到的参数最终趋于稳定,为利用带遗忘因子的最小二乘算法,曲线参数最终还是有小幅度震荡。
系统辨识大作业.
一、 问题描述考虑仿真对象:()0.9(1)0.15(2)0.02(3)0.7(1)0.15(2)()z k z k z k z k u k u k e k +-+-+-=---+ e() 1.0e(1)0.41e(2)(),~(0,1)k k k v k v N λ+-+-=式中,u(k)和z(k)是输入输出数据,v(k)是零均值、方差为1的不相关的随机噪声;u(k)采用与e(k)不相关的随机序列。
1. 设计实验,产生输入输出数据;2. 使用基本最小二乘法估计参数;3. 考虑其他适用于有色噪声的辨识方法估计参数;4. 模型验证。
二、 问题分析对于单输入单输出系统(Single Input Single Output, SISO ),如图 1所示,将待辨识的系统看作“灰箱”,它只考虑系统的输入输出特性,而不强调系统的内部机理。
图 1中,输入u(k)和输出z(k)是可以测量的,1()G z -是系统模型,用来描述系统的输入输出特性,y(k)是系统的实际输出。
1()N z -是噪声模型,v(k)是均值为零的不相关随机噪声,e(k)是有色噪声。
图 1 SISO 系统的“灰箱”结构对于SISO 随机系统,被辨识模型()G z 为:12121212()()()1n n nn b z b z b z y z G z u z a z a z a z ------+++==++++ 其相应的差分方程为11()()()n ni i i i y k a y k i b u k i ===--+-∑∑若考虑被辨识系统或观测信息中含有噪声,被辨识模型可改写为11()()()()n ni i i i z k a y k i b u k i v k ===--+-+∑∑式中,z(k)为系统输出量的第k 次观测值;y(k)为系统输出量的第k 次真值,y(k-1)为系统输出量的第k-1次真值,以此类推;u(k)为系统的第k 个输入值,u(k-1)为系统的第k-1个输入值;v(k)为均值为0的不相关随机噪声。
(完整word版)2003版系统辨识最小二乘法大作业
西北工业大学系统辩识大作业题目:最小二乘法系统辨识一、 问题重述:用递推最小二乘法、加权最小二乘法、遗忘因子法、增广最小二乘法、广义最小二乘法、辅助变量法辨识如下模型的参数离散化有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 是一套高性能数字计算和可视化软件 ,它集成概念设计 ,算法开发 ,建模仿真 ,实时实现于一体 ,构成了一个使用方便、界面友好的用户环境 ,其强大的扩展功能为各领域的应用提供了基础。
对于一个简单的系统 ,可以通过分析其过程的运动规律 ,应用一些已知的定理和原理,建立数学模型 ,即所谓的“白箱建模 ”。
但对于比较复杂的生产过程 ,该建模方法有很大的局限性。
由于过程的输入输出信号一般总是可以测量的 ,而且过程的动态特性必然表现在这些输入输出数据中 ,那么就可以利用输入输出数据所提供的信息来建立过程的数学模型。
(完整)系统辨识大作业汇总,推荐文档
参数递推估计是指被辨识的系统,每取得一次新的测量数据后,就在前一 次估计结果的基础上,利用新引入的测量数据对前一次估计的结果进行修正, 从而递推地得出新的参数估计值。这样,随着新测量数据的引入,一次接一次 地进行参数估计,直到估计值达到满意的精确程度为止。最小二乘递推算法的 基本思想可以概括为:
当前的估计值ˆ(k) =上次估计值ˆ(k 1) +修正项 即新的估计值ˆ(k) 是在旧的估计值ˆ(k 1) 的基础上,利用新的观测数据对旧的 估计值进行修正而成的。
可以看出,取 (k) 1的时候,加权最小二乘估计就退化成了最小二乘参数 估计的递推算法(Recursive Least Squares, RLS)。加权参数 1 可以在
(0,1]范围内选择,如果选 1 1,所有的采样数据都是等同加权的,如果
(k)
1 1,则表示对新近获得的数据给予充分大的加权因子,而削弱历史观测 (k)
可以根据生成的白噪声序列和输入序列,以及必要的 0 初始值,带入表 达式即可得到采样输出数据。
2. 差分模型阶检验 在实际场景中,辨识模型的阶数和纯时延往往是未知的,在很多情况下仅
仅依靠猜测。在模型的阶数和纯时延不确定时,设系统模型为
n
n
y(t) ai y(t i) bj y(t i) (t)
数据的影响。 实际计算时,需要首先确定初始参数ˆ(0) 和 P(0) 。
P(0) 2I 为充分大实数
一般说来选取
(0)
为充分小的向量
对于这样的系统,使用最小二乘法参数估计的递推算法进行辨识可以得到 无偏估计,但是如果噪声模型必须用 C(z1)v(k) 表示时,此时就无法得到无偏估 计了,因为该方法没有把噪声模型考虑进去。
K (k) P(k 1)h(k)[hT (k) p(k 1)h(k) 1 ]1
系统辨识相关分析和最小二乘
对于一个脉冲响应函数为g(t)的线性系统,当以 输入信号u(t)的自相关函数Ruu(τ)作为系统输 入时,则系统的输出即为输入u和输出z之间的互相 关函数Ruz(τ)。
如果输入信号是白噪声的话,由白噪声的性质得:
求解脉冲响应函数问题简化成了计算互相关函数问 题
需要解决积分时间长和白噪声的物理实现这两个问 题。 为了解决第一个问题,即在有限的时间内,完成互 相关函数的计算,可以采用周期白噪声;为了解决 第二个问题,可采用近似白噪声信号。
1 n
ˆ n N y
b0 ˆ bn
设e(k)=y(k)- ŷ(k), e(k)称为残差,则有e=y- ŷ=y-Φθ 最小二乘估计要求残差的平方和最小,即按照指数 T T ˆ ˆ J e e y y 函数: ˆ 的偏导数并令其等于0可得: 则求J对
y N u n N
u N
令:
a1 y n 1 n 1 y n 2 a n 2 , n , y b0 y n N n N bn y 1 u n 1 y n y n 1 y 2 u n 2 y N u n N y n N 1
PN 1 PN PN N 1 1 P
1 T T N 1 N N 1 N 1 N
本次试验差分方程的参数真值为:
a1 1.5; a2 0.7; b1 1; b2 0.50;
实验结果图:
递推最小二乘法相关辨识方法(97-2003)
定义 式中
k 1 T T P k H k H k h i h i i 1 k 1 P 1 k 1 H T H h i hT i k 1 k 1 i 1
hT 1 T h 2 Hk T h k
把P(k)写成
利用矩阵反演公式 可得
P k P
1
1
k 1 h k h k
T
1
1
A CCT A1 A1C I CT A1C CT A1
1
P k P k 1 P k 1 h k hT k P k 1 hT k P k 1 h k 1 P k 1 h k hT k I T P k 1 h k P k 1 h k 1
z k a1z k 1 ana z k na bu k 1 bnb u k nb e k 1
令k=1,2,3,…,L,可构成线性方程组 zL H L eL 式中
z 1 e 1 z 2 e 2 zL , eL z L e L u 0 u 1 nb z 0 z 1 na z 1 z 2 na u 1 u 2 nb HL z L 1 z L na u L 1 u L nb
hT 1 T h 2 H k 1 T h k 1
式中,h(i)是一个列向量,也就是HL 的第i行的倒置, P(k)是一个方阵,它的维数取决于未知参数的个数。
系统辨识作业及答案
系统辨识作业及答案一.问答题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(0100112100111 10011110011109()()()()()()()001112401110)23(111012********* 010020010011910011180011117()()()()()()()()10011320011131011103000111291101028101002701001261001125 其输出序列为:1 1 1 1 1 0 0 1 0 1⑵不是M 序列⑶第4级与第3级模2相加结果100108001007010006100015000114001113011112111111)()()()()()()()(()()()()()()()11110161110115110101410101)13(0101112101101 10110010110019()()()()()()()110012410010)23(001002201000211 000120000111900111180111117()()()()()()()()01111321111031111013011010291010128010112710110260110025 不同点:第2级和第3级模二相加产生的序列,是从第4时刻开始,每隔7个时刻重复一次;第4级与第3级模2相加产生的,序列,是从第2时刻开始每隔15个时刻重复一次。
系统辨识第5章 线性动态模型参数辨识 最小二乘法
度函数
,则称uS(uk()为) “持续激励”信号。
● 定义4 一个具有谱密度 Fn (为z 1的) 平f1z稳1 信f2号z 2u(k)称fn为z nn 阶
“持续激励”Fn信(e号j ),2 S若u (对) 一0 切形如 Fn (e j ) 0
的滤波器,关系式
,意味着
。
● 定理2 设输入信号u(kR)u是(0)平稳R随u (1机) 信号,Ru (如n 果1)相关函数矩阵
式中
zL H L nL
nzHLLL[[zn(h(hh11TT)T),((,(zL12n())()22)),,,,znz(((LzLzL)(()]10]))1)
z(1 na ) z(2 na )
z(L na )
u(0) u(1)
u(L 1)
u(1 nb )
u(2
nb
)
u(L nb )
5.2 最小二乘法的基本概念
● 两种算法形式
① 批处理算法:利用一批观测数据,一次计算或经反复迭代,
以获得模型参数的估计值。
②
递推算法:在上次模型参数估计值
ˆ
(k
1)的基础上,根据当
前获得的数据提出修正,进而获得本次模型参数估计值ˆ (k ),
广泛采用的递推算法形式为
(k ) (k 1) K (k )h(k d )~z (k )
z(k ) h (k ) n(k )
式中z(k)为模型输出变量,h(k)为输入数据向量, 为模型参
数向量,n(k)为零均值随机噪声。为了求此模型的参数估计值, 可以利用上述最小二乘原理。根据观测到的已知数据序列
和{z(k)} ,{h极(k小)} 化下列准则函数
L
J ( ) [z(k ) h (k ) ]2
- 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。