系统辨识大作业1201张青
上机实习Ⅰ：系统辨识方法初识
上机实习Ⅰ：系统辨识方法初识12自动化许天野12350068指导老师：王国利摘要：系统辨识、状态估计和控制理论是现代控制论中相互渗透的三个领域。
控制理论的应用离不开系统辨识技术，实际中，许多控制系统的模型在工作中是变化的，为了实现自适应控制，需要系统辨识技术不断更新模型参数。
通过学习使用MATLAB软件，初步体验系统辨识方法。
关键字：系统辨识，控制理论，MATLAB。
Practice 1：Practice the method of System identification Abstract: System identification, State estimation and The Principle of Automatic Control are three different disciplines of the modern control theory, which are interpenetrated with one another. In practice, the model of system is changing all the time. To control adaptively, the system model should be update its parameters, by the method of System identification. By learning the using of MATLAB, we are supposed to practice the method of system identification.Key Words: System identification, System identification, MATLAB目录一、引言 (3)1.1介绍 (3)1.2实验目的 (3)二、实验内容和方法 (3)2.1实验内容 (3)2.2实验步骤 (3)2.2.1输入信号选择 (4)2.2.2 数据收集 (4)2.2.3 实验步骤 (4)三、实验结果 (5)四、实验分析探究 (7)4.1分析 (7)4.2探究 (7)4.3结果分析 (10)一、引言1.1介绍在自然科学和社会科学的许多领域中，人们越来越重视对系统进行定量的系统分析、系统综合、仿真、控制和预测。
系统辨识基础--经典辨识方法
其中
Lh i 1t s1c1 sc2s2 1 ci 1 si 1
h
9
进一步利用下式
e s t 1 s ts2 t2 si ti
1 ! 2 !
i!
可得 得
L1h*t 1h*testdt Misi
0
i0
Mi
0
1h*t
ti
i!
dt
1Aisi1Misi11
i1 i0
a4 4 4.1207
b1 7.5 7.50402
b2 17.5 17.5233
h
15
4.3 脉冲响应法
ut
yt
1 ut
过 程 yt
0
t
0
t
gk1hkhk1
T0
h
16
ut
过程
yt
gt,0
模型参数 调整机构
~yt +
-
模型
gt,
图4.6 “学习法”原理
h
17
由脉冲响应求过程的传递函数-一阶过程
条件
增益K
a1
噪声 情况
无测量噪 声
有测量噪 声（方差 为0.01）
采样时间 4秒 1.5秒
1.5秒
数据长度 12 30
30
1.0 0.999984 0.999965
1.00204
10.0 11.7097 10.2171
11.5776
a2
6.5 6.52053 6.49897
6.47451
参数
真值 估 计 值
h
14
例 4.2
G s4s41 1s.5 5 3 7s 21 .7 5 7 .s 52 s 7 1.5s1
安大系统辨识作业
2012系统辨识实验姓名：周自飞学号：P4*******年级专业：09自动化线性系统参数估计的最小二乘法一、实验目的：1、掌握线性离散系统的数学模型；2、掌握线性离散系统在无噪声时的单位脉冲响应；3、掌握线性系统参数估计的最小二乘法。
二、实验原理：1、基本最小二乘算法2、递推最小二乘算法3、渐消记忆递推最小二乘算法4、辅助变量递推最小二乘算法5、增广递推最小二乘算法三、实验内容：设带有噪声的离散系统模型为：()()()()()()k213.0-13.0-+.028.11-412+ukkukyk--yξky=现给定系统的一列长度为50的输入序列：u(k)= 1,0.3,-0.5,0.9,-0.5,-0.3,-0.2,0.4,0.8,-0.6,-0.1,0,0.1,0.5,-0.6,-0.2,0.3,0.9,0.5,0.2,-0.6,-0.3,-0.1,0.7,0,0.3,-0.6,0.3,0.1,0.5,-0.7,-0.4,-0.9,-0.6,-0.2,-0.4,-0.2,0.1,-0.1,0.1,0.9,0.5,0.3,0.7,0.4,-0.2,-0.7,-0.2,0.1,01. 构造离散系统模型函数，给出无噪声时系统在单位脉冲输入下的响应序列；%zhouzifei.mb=[0.3 -0.213];a=[1 -1.28 0.41];impz(b,a,30);title(‘系统单位脉冲响应’);axis([-1 30 -0.2 0.5]);2. 给出系统在给定输入u(k)和噪声()kξ下的输出信号y(k)；%zifei2.m3. 利用已有输入信号u(k)和在第2步中得出的输出信号y(k)，使用下列方法辨识系统的参数：(1) 基本最小二乘算法%LS.m参数：a1 =-1.2803a2 =0.4138b0 = 0.2867b1 = -0.2104(2) 递推最小二乘算法%RLS.m参数：a1 =-1.2800 a2 =0.4100 b0 = 0.3000 b1 = -0.2130(3) 渐消记忆递推最小二乘算法%FMRLS参数：a1 =-1.2800 a2 =0.4100 b0 = 0.3000 b1 = -0.2130(4) 辅助变量递推最小二乘算法%IVLS.m参数：a1 =-1.2800 a2 =0.4100b0 = 0.3000 b1 = -0.2130(5) 增广递推最小二乘算法（可选做，有加分）%ARLS.m参数：a1 =-1.2800 a2 =0.4100 b0 = 0.3000 b1 = -0.2130四、实验结果：最小二乘法思想是使各次观测值和计算值之间差值的平方乘以度量其精确度的数值后的和为最小。
系统辨识大作业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 可由仿真图得出，可知两种方法参数确定十分接近。
系统辨识步骤及内容
系统辨识步骤及内容系统辨识是研究如何用实验研究分析的办法来建立待求系统数学模型的一门学科。
Zadeh(1962)指出：“系统辨识是在输入和输出数据的基础上，从一类模型中确定一个与所观测系统等价的模型”。
Ljung(1978)也给出如下定义：“系统辨识有三个要素——数据、模型类和准则，即根据某一准则，利用实测数据，在模型类中选取一个拟合得最好的模型”。
实际上，系统的数学模型就是对该系统动态本质的一种数学描述，它向人们提示该实际系统运行中的有关动态信息。
但系统的数学模型总比真实系统要简单些，因此，它仅是真实系统降低了复杂程度但仍保留其主要特征的一种近似数学描述。
建立数学模型通常有两种方法，即机理分析建模和实验分析建模。
机理分析建模就是根据系统内部的物理和化学过程，概括其内部变化规律，导出其反映系统动态行为并表征其输入输出关系的数学方程（即机理模型）。
但有些复杂过程，人们对其复杂机理和内部变化规律尚未完全掌握（如高炉和转炉的冶炼过程等）。
因此，用实验分析方法获得表征过程动态行为的输入输出数据，以建立统计模型，实际上是系统辨识的主要方面，它可适用于任何结构的复杂过程。
系统辨识的主要步骤和内容有以下几个方面。
1、辨识目的根据对系统模型应用场合的不同，对建模要求也有所不同。
例如，对理论模型参数的检验及故障检测和诊断用的模型则要求建得精确些。
而对于过程控制和自适应控制等用的模型的精度则可降低一些，因为这类模型所关心的主要是控制效果的好坏，而不是所估计的模型参数是否收敛到真值。
2、验前知识验前知识是在进行辨识模型之前对系统机理和操作条件、建模目的等了解的统称。
有些场合为了获得足够的验前知识还要对系统进行一些预备性的实验，以便获得一些必要的系统参数，如系统中主要的时间常数和纯滞后时间，是否存在非线性，参数是否随时间变化，允许输入输出幅度和过程中的噪声水平等。
3、实验设计实验设计的主要内容是选择和决定：输入信号的类型、产生方法、引入点、采样周期、在线或离线辨识、信号的滤波等。
系统辨识大作业
系统辨识大作业专业班级：自动化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—N采样时刻的输入输出数据，即
u (0), u (1)," , u ( N − 1), u ( N ) y (0), y (1)," , y ( N − 1), y ( N )
假定系统的模型属于如下的模型类：
y ( k ) + ay ( k − 1) = bu (k − 1) + v(k )
k =1
N
∂V (θ ) N = ∑ 2ay 2 (k − 1) + 2 y (k ) y (k − 1) − 2by (k − 1)u (k − 1) ∂a k =1 ∂V (θ ) N = ∑ 2bu 2 (k − 1) − 2 y (k )u (k − 1) − 2ay (k − 1)u (k − 1) ∂b k 等：子空间辨识
1990年代，为了克服PEM针对多变量系统辨识
时需要进行非线性优化，以及IV不能同时辨识 出噪声模型的缺点。Bart De Moor, Verhaegen 等提出了针对多变量系统的subspace identification methods。该类方法不是基于优化 某个criterion，主要用到矩阵的奇异值分解， 无需非线性优化，因而计算量较小。
1.2 模型
数学模型是用来描述系统行为的数学语
言。 非线性系统的数学模型是非线性状态方 程和输出方程。线性系统的数学模型可 以有多种相互等价的形式：状态空间方 程、传递函数、阶跃响应、差分方程等。
扰 动 输入
系统
输出
1.3 建模的两大类方法
机理分析法(first principles modeling）或称为白
何求取参数估计值。least-squares, prediction error, instrumental variable 参数估计算法的统计性质：无偏性、一致性。 如何验证所得模型的有效性？如何选择模型阶数？
系统辨识大作业.
一、 问题描述考虑仿真对象：()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的不相关随机噪声。
系统辨识的基本概念
系统正确描述系统动态性能的数学摸型——就成了自 动控制 理论 和工程实践的重要组成部分。
系统辨识就是从对系统进行观察和测量所获得的信
息重提取系统数学模型的一种理论和方法。日渐成熟。
建模——成为各门学科的共同语言。
系统辨识的基本概念
2
1.1 系统和模型
1.1.1 系统
（system/process）
到95％时的调节时间。
26
系统辨识的基本概念
4、数据的零值化处理
•差分法(Isermann,1981)
•平均法
•剔除高频成分（一般采用低通滤波器）
5、模型结构辨识
模型验前结构的假定、模型结构参数的确定。
6、模型参数辨识（本课程的主要内容）
当模型结构确定后，进行的就是模型参数辨识
7、模型检验
模型检验是辨识不可缺少的步骤。常用的有“白色度”检验
18
系统辨识的基本概念
● 误差准则
L
J() f ((k))
k1
也叫等价准则、误差准则、损失函数或准则函数。
用的最多的是： f((k))2(k)
● 输出误差准则： ( k ) z ( k ) z m ( k ) z ( k )[ u ( k )]
● 输入误差准则： ( k ) u ( k ) u m ( k ) u ( k ) 1 [z ( k )]
12
系统辨识的基本概念
又置：
logP(k) logV (k) logc
令
y(k) z(k)
logP(k),1 logV (k),2
logc
h(k) [z(k),1]t
[1,2]
则y(k和 ) h(k)都是可观测的变量应，的对最小二乘格式
(完整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 是一套高性能数字计算和可视化软件 ,它集成概念设计 ,算法开发 ,建模仿真 ,实时实现于一体 ,构成了一个使用方便、界面友好的用户环境 ,其强大的扩展功能为各领域的应用提供了基础。
对于一个简单的系统 ,可以通过分析其过程的运动规律 ,应用一些已知的定理和原理,建立数学模型 ,即所谓的“白箱建模 ”。
但对于比较复杂的生产过程 ,该建模方法有很大的局限性。
由于过程的输入输出信号一般总是可以测量的 ,而且过程的动态特性必然表现在这些输入输出数据中 ,那么就可以利用输入输出数据所提供的信息来建立过程的数学模型。
系统辨识基础教学大纲
主要内容： 1. 阶跃响应法 2. 脉冲响应法 3. 频率响应法 4. 相关分析法
基本要求： 掌握阶跃响应法、脉冲响应法、频率响应法 了解相关分析法基本原理
学时分配： 8 学时 第三章 最小二乘法
主要内容： 1. 最小二乘法问题的解 2. 最小二乘法问题的递推算法 3. 偏差补偿最小二乘法 4. 增广最小二乘法
基本要求： 掌握最小二乘法的基本概念、最小二乘法问题的解 掌握最小二乘法问题的递推算法 掌握偏差补偿最小二乘法 了解增广最小二乘法的基本原理
了解最小二乘法问题的统计性质 学时分配：8 学时
数辨识程序 六、教学参考书 1．方崇智主编，《过程辨识》（第一版），清华大学出版社，1998 年 2．刘 豹主编，《系统辨识》（第二版），机械工业出版社，1996 年 3．王秀峰主编，《系统建模与辨识》（第一版），电子工业出版社，2004 年
《系统辨识基础》教学大纲
课程名称：系统辨识基础（Basic of System Discrimination） 课程编码：152050 学 分：2 分 总 学 时：32 学时，理论学时：22 学时；上机学时：10 学时 适用专业：自动化专业 先修课程：自动控制理论、现代控制理论 一、课程的性质、目的与任务
第四章 梯度校正参数辨识方法 主要内容：
1. 确定性问题的梯度校正参数辨识方法 2. 随机问题的梯度校正参数辨识方法 基本要求： 掌握确定性问题的梯度校正参数辨识方法 了解随机问题的梯度校正参数辨识方法的基本原理 学时分配：4 学时 三、上机内容与学时分配 上机内容： 给定一个模拟或数字对象，编写计算机程序对其参数进行辨识 学时分配：10 学时 四、大纲说明 1． 本课程是自动化专业的一门专业选修课程。 2． 本课程的先修课程为自动控制理论、现代控制理论。 3． 本课程安排上机大作业，要求学生结合课程学习，编写调试数字、模拟对象的参
(完整)系统辨识大作业汇总,推荐文档
一、 问题描述考虑仿真对象：()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 n n b z b z b z y z G z u z a z a z a z------+++==++++L L 其相应的差分方程为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的不相关随机噪声。
《系统辨识》实验手册
《系统辨识》实验手册第二炮兵工程大学控制工程系2015年5月目录实验1 白噪声和M序列的产生---------------------------------------------------------- 2 实验2 相关分析法辨识脉冲响应------------------------------------------------------- 5 实验3 最小二乘法的实现--------------------------------------------------------------- 9 实验4 递推最小二乘法的实现---------------------------------------------------------- 12 附录实验报告模板---------------------------------------------------------------------- 16实验1 白噪声和M 序列的产生一、实验目的1、熟悉并掌握产生均匀分布随机序列方法以及进而产生高斯白噪声方法2、熟悉并掌握M 序列生成原理及仿真生成方法二、实验原理1、混合同余法混合同余法是加同余法和乘同余法的混合形式，其迭代式如下：111(*)mod /n n n n x a x b MR x M +++=+⎧⎨=⎩ 式中a 为乘子，0x 为种子，b 为常数，M 为模。
混合同余法是一种递归算法，即先提供一个种子0x ，逐次递归即得到一个不超过模M 的整数数列。
2、正态分布随机数产生方法由独立同分布中心极限定理有：设随机变量12,,....,,...n X X X 相互独立，服从同一分布，且具有数学期望和方差：2(),()0,(1,2,...)k k E X D X k μσ==>=则随机变量之和1nk i X =∑的标准化变量:()nnnkk kXE X Xn Y μ--==∑∑∑近似服从(0,1)N 分布。
系统辨识最小二乘法大作业
系统辨识大作业最小二乘法及其相关估值方法应用学号：2012302259姓名：王家琦基于最小二乘法的多种系统辨识方法研究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)可得残差，，，。
11班张青期末管理反思
_11班学生教育管理班主任：张青问题一：当前班级学生的现状与全年级比存在的优势、差距分别是什么？与平行班比存在的优势、差距分别是什么？（精神状态、行为表现、层次划分、突出问题等方面）2012级11班有60名学生，33名男生，27名女生；1人借读；1人即将出国呢，2人长期生病。
班上大多数同学纪律意识不强，在校学习习惯很一般，学习态度不很端正，尤其中下游同学学习习惯有待进一步培养与加强。
尤其从南面科技楼搬到北楼，尤其是一楼，学生很不适应，仿佛从与世无扰的净土跌落到人间，变得爱说爱动，下课就往外窜，上课铃响才匆匆赶回教室。
这是一个比较严峻的问题。
经过两个月全体师生的奋战，本次期中考试，优生状况有所提升：前100名9人，前150名14人（加上体育分15人）。
与平行班级相比，优生数量基本能超过平均数，但中游学生很少，有20名同学位于学困生行列，且600名以后占了14人，导致班级总平均分不很理想。
问题二：进入初三以后，学生管理面临哪些问题？采取了哪些措施？效果如何？下一步有什么计划？接过初三11班，明显感觉部分男同学自制能力特差，时常不能控制自己，上课时有8人爱随便说话或者做小动作。
女同学有三人比较叛逆，一人爱和外班同学交往，经常撒谎，成绩一直是班里倒数第一；一人天天带零食，买惊悚小说，后有所改观。
积极回答问题的孩子集中在前3分之1孩子身上。
个别学生的基础很差，作业完成拖拉、懒散，完成不够主动，质量不够好。
有相当一部分学生不让家长过问，拒绝让家长签字。
措施：以班规约束学生行为，从而激励同学们激流勇进，目光长远；严格班级纪律，从出勤、仪容仪表、卫生值日、自习课秩序，到教室内的卫生行为规范等，都进行了细致的量化管理，班里都有明确要求和记录，培养学生的主人公意识、大局意识。
秉持“对尖子不捧不惯，对差生不离不弃，严而有度，爱而不纵”的原则，尽量让学生有序、自主管理并养成习惯；学习上主要培养优生的竞争意识、中下游学生的进步意识。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
《系统辨识》大作业学号：********班级：自动化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 可由仿真图得出，可知两种方法参数确定十分接近。
两种方法的辨识效果差异:采用相关分析法和最小二乘法辨识出的脉冲响应图形可看出，用最小二乘法辨识出的图形更像脉冲响应，他在最后无偏差，衰减为零，其可实现无偏差估计。
而相关分析法图形虽与相关分析的相近，但它最后有偏差。
带遗忘因子的递推最小二乘辨识的参数平均值可随着实际参数变化而突变。
但他对噪声比较敏感，只是参数围绕参数实际值上下波动，而辨识出参数的平均值接近实际值，而且比其他方法更加接近，并且可以防止数据饱和。
第二题设SISO系统差分方程为（40分）y(k) =)()2()1()2()1(2121kkubkubkyakyaξ+-+-+----辨识参数向量为θ＝1[a2a1b2b ]T，输入输出数据详见数据文件uy01.txt—uy03.txt。
)(kξ为噪声方差各异的白噪声或有色噪声。
试求解：1）用最小二乘及递推最小二乘法估计θ；2）用辅助变量法及其递推算法估计θ；3）用广义最小二乘法及其递推算法估计θ；4）用夏氏偏差修正法、夏氏改良法及其递推算法估计θ；5）用增广矩阵法估计θ；6) 用极大似然法估计θ；7）分析噪声)(kξ特性；答：1.基本最小二乘法：%首先将TXT中的数据装载在Matlab中：clear;clc;Wbs = 0;while Wbs<1disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt') filename = input('输入所需打开文件名:\','s');[Wbs,Message] = fopen(filename,'r');%打开文件；endData = fscanf(Wbs,'%g %g',[2 inf]);%生成数据矩阵Data = Data';%最小二乘法，取b0=0n = 2;%根据题目，本系统为2阶系统L = length(Data);%辨识所需数据的总长度N = L-n;%构造测量矩阵的数据长度FIA = zeros(N,2*n);%构造测量矩阵for i = 1:N %测量矩阵赋值for l = 1:n*2if(l>n)FIA(i,l) = Data(n+2+i-l,1);elseFIA(i,l) = -Data(n+i-l,2);endendendY = Data(n+1:n+N,2);%输出数据矩阵thita = inv(FIA'*FIA)*FIA'*Y;%计算参数矩阵disp('使用最小二乘算法所得结果如下：')%辨识参数数据输出如下：lsa1 = thita(1),lsa2 = thita(2),lsb1 = thita(3),lsb2 = thita(4)实验结果：Uy1.txt uy2.txt uy3.txt2.递推最小二乘法：%首先将TXT中的数据装载在Matlab中：clear;clc;Wbs = 0;while Wbs<1disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt') filename = input('输入所需打开文件名:\','s');[Wbs,Message] = fopen(filename,'r');%打开文件；endData = fscanf(Wbs,'%g %g',[2 inf]);%生成数据矩阵Data = Data';%递推最小二乘法n = 2;%根据题目已知该系统阶数为2L = length(Data);%为辨识参数向量赋初值，这里均取为0.001for a = 1:1:n*2thita0(a,1) = 0.001;end%直接给出矩阵Pn的初始状态P0P0 = 10^9*eye(n*2,n*2);thita = [thita0,zeros(n*2,L)];%需辨识参数矩阵的初值及大小for i = n+1:1:L %这里i从第n+1个数开始for m = 1:n*2 %输入矩阵赋值if(m<=n)FIA1(m,1) = -Data(i-m,2);elseFIA1(m,1) = Data(i-m+2,1);endendX = FIA1'*P0*FIA1+1;%引入中间变量X为K递推公式中的求逆项,1维K1 = P0*FIA1/X;D1 = Data(i,2)-FIA1'*thita0;P1 = P0-K1*K1'*X;%求出P(K)矩阵P0 = P1;%作为下次迭代初值thita1 = thita0+K1*D1;%估值补偿公式，thita1为求被辨识的参数向量thita0 = thita1;%所得的参数向量作为下次迭代初值enddisp('递推最小二乘法所得结果如下：')Dta1=thita1(1),Dta2=thita1(2),Dtb1=thita1(3),Dtb2=thita1(4) %显示被辨识参数实验结果：Uy1.txt uy2.txt uy3.txt3.辅助变量法：%首先将TXT中的数据装载在Matlab中：clear;clc;Wbs = 0;while Wbs<1disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt')filename = input('输入所需打开文件名:\','s');[Wbs,Message] = fopen(filename,'r');%打开文件；endData = fscanf(Wbs,'%g %g',[2 inf]);%生成数据矩阵Data = Data';%使用递推辅助变量法对2阶系统进行辨识，取b0=0n = 2;L = length(Data);N = L-n;%分别取出输入、输出数据u,yU = Data(1:L,1);Y = Data(1:L,2);%首先使用最小二乘法，粗略得到初始的被辨识参数向量fzOL = [-Y(2:L-1),-Y(1:L-2),U(2:L-1),U(1:L-2)];Y1 = Data(3:L,2);%构造输出向量，也用于辅助变量法的输出向量thita = inv(fzOL'*fzOL)*fzOL'*Y1;%得到初始参数向量%获得初值，以下使用辅助变量法for i=1:400Yf = fzOL*thita;%辅助模型输出%构造辅助变量矩阵，改变新的矩阵中的输出量，输入量不变Zfz = fzOL;Zfz(2:N,1) = -Yf(1:(N-1),1);Zfz(3:N,2) = -Yf(1:(N-2),1);thita = inv(Zfz'*fzOL)*Zfz'*Y1;%求取被估计参数的辅助变量估值enddisp('使用辅助变量法，得到如下辨识结果：')Fza1=thita(1),Fza2=thita(2),Fzb1=thita(3),Fzb2=thita(4)实验结果：Uy1.txt uy2.txt uy3.txt4.递推辅助变量法：%首先将TXT中的数据装载在Matlab中：clear;clc;Wbs = 0;while Wbs<1disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt') filename = input('输入所需打开文件名:\','s');[Wbs,Message] = fopen(filename,'r');%打开文件；endData = fscanf(Wbs,'%g %g',[2 inf]);%生成数据矩阵Data = Data';%使用递推辅助变量法对2阶系统进行辨识，取b0=0n = 2;L = length(Data);%分别取出输入、输出数据u,yU = Data(1:L,1);Y = Data(1:L,2);N = 498;n0 = 100;% 取前100个数据利用最小二乘辨识FIA = [-Y(2:n0-1),-Y(1:n0-2),U(2:n0-1),U(1:n0-2)];Y1 = Data(3:n0,2);%构造输出向量，也用于辅助变量法的输出向量thita = inv(FIA'*FIA)*FIA'*Y1;%得到初始参数向量Z = [-Y(2:n0-1) -Y(1:n0-2) U(2:n0-1) U(1:n0-2)];thita = inv(Z'*FIA)*Z'*Y1;P0 = inv(Z'*FIA);%后面的数据计算，使用递推辅助变量法for n = n0+1:NY11 = [Y(1); Y(2); Z*thita];Zn = [-Y11(n-1) -Y11(n-2) U(n-1) U(n-2)];Z = [Z; Zn];kesy = [-Y(n-1); -Y(n-2); U(n-1); U(n-2)];K = P0*Zn'/(1+kesy'*P0*Zn');P1 = P0 - K*kesy'*P0;P0 = P1;thita = thita + K*(Y(n) - kesy'*thita);enddisp('递推辅助变量法辨识结果：');Dfa1=thita(1),Dfa2=thita(2),Dfb1=thita(3),Dfb2=thita(4)实验结果：Uy1.txt uy2.txt uy3.txt5.广义最小二乘法：%首先将TXT中的数据装载在Matlab中：clear;clc;Wbs = 0;while Wbs<1disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt')filename = input('输入所需打开文件名:\','s');[Wbs,Message] = fopen(filename,'r');%打开文件；endData = fscanf(Wbs,'%g %g',[2 inf]);%生成数据矩阵Data = Data';%使用广义最小二乘法，辨识2阶系统的参数n = 2;L = length(Data);N = L-n;for m=1:1:300%分别取出输入、输出数据u,yU = Data(1:L,1); Y = Data(1:L,2);%首先使用最小二乘法，求系统参数向量初值gyOL = [-Y(2:L-1),-Y(1:L-2),U(2:L-1),U(1:L-2)];%LS测量矩阵Zgy1 = Data(3:L,2);%输出矩阵ythita = inv(gyOL'*gyOL)*gyOL'*Zgy1;%参数向量初值Gya1=thita(1);Gya2=thita(2);Gyb1=thita(3);Gyb2=thita(4);%得到初值后，以下使用广义最小二乘法进行辨识，由残差代替噪声误差E(1) = Y(1);%k=1时，残差的值E(2) = Y(2)+Gya1*Y(1)-Gyb1*U(1);%k=2时，残差的值for i = 3:N+n %残差公式如下：E(i)= Y(i)+Gya1*Y(i-1)+Gya2*Y(i-2)-Gyb1*U(i-1)-Gyb2*U(i-2);endE1 = E(n+1:n+N)';omiga(1:N,1) = -E(n:n+N-1)';omiga(1:N,2) = -E(n-1:n+N-2)';%fl表示由残差代替噪声项，而辨识得到的滤波器的估计参数fl=inv(omiga'*omiga)*omiga'*E1;%得到fl后，计算经过滤波的u，y数据；其中，k=1,2时：Data(1,1)=U(1);Data(1,2)=Y(1);Data(2,2)=Y(2)+fl(1)*Y(1);Data(2,1)=U(2)+fl(1)*U(1);%k>=3时：for k=3:N+nData(k,1) = U(k)+fl(1)*U(k-1)+fl(2)*U(k-2);Data(k,2) = Y(k)+fl(1)*Y(k-1)+fl(2)*Y(k-2);endenddisp('由广义最小二乘法得到的辨识结果如下：')Gya1=thita(1),Gya2=thita(2),Gyb1=thita(3),Gyb2=thita(4)实验结果：Uy1.txt uy2.txt uy3.txt6.递推广义最小二乘法：%首先将TXT中的数据装载在Matlab中：clear;clc;Wbs = 0;while Wbs<1disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt') filename = input('输入所需打开文件名:\','s');[Wbs,Message] = fopen(filename,'r');%打开文件；endData = fscanf(Wbs,'%g %g',[2 inf]);%生成数据矩阵Data = Data';u1=Data(:,1);y1=Data(:,2);n=2;N=500-n;n0=50;%用前50组数据利用lsm估计出初值%构造矩阵phiphi(:,1)=-y1(n:n+n0-1);phi(:,2)=-y1(n-1:n+n0-2);phi(:,3)=u1(n:n+n0-1);phi(:,4)=u1(n-1:n+n0-2);y(:,1)=y1(n+1:n+n0);seta=(inv(phi'*phi))*phi'*y;seta0=seta;f0=[0 0]';p10=inv(phi'*phi);p20=10^6*eye(2,2);y(1:i,1)=y1(n+1:n+i);u(1:i,1)=u1(n+1:n+i);%计算残差e=y-phi*seta0;%f 估计omega=[-e(n+i-2) -e(i+1-1)]';K2=p20*omega*inv(1+omega'*p20*omega);f1=f0+K2*(e(n+i-2)-omega'*f0);p2=p20-K2*omega'*p20;p20=p2;f0=f1;%滤波yy(1:i,1)=y1(n:n+i-1);yy(1:i,2)=y1(n-1:n+i-2);uu(1:i,1)=u1(n:(n+i-1));uu(1:i,2)=u1((n-1):(n+i-2));yg=y+yy*f1; %yg(k)=y(k)+f(1)*y(k-1)+f(2)*y(k-2)ug=u+uu*f1; %ug(k)=u(k)+f(1)*u(k-1)+f(2)*u(k-2)% seta 估计ksai=[-yg(n+i-2) -yg(n+i-3) ug(n+i-2) ug(n+i-3)]';K1=p10*ksai*inv(1+ksai'*p10*ksai);% K(N+1)值p1=p10-K1*ksai'*p10; % P(N+1)值p10=p1;seta1=seta0+K1*(y1(n+i+1)-ksai'*seta0); % seta(N+1)值seta0=seta1;% 构造新phi 阵phi(1:i+1,1)=-y1(n:n+i-1+1);phi(1:i+1,2)=-y1(n-1:n+i-2+1);phi(1:i+1,3)=u1(n:n+i-1+1);phi(1:i+1,4)=u1(n-1:n+i-2+1);enddisp('递推广义最小二乘法辨识结果：');a1=seta0(1),a2=seta0(2),b1=seta0(3),b2=seta0(4)实验结果：Uy1.txt uy2.txt uy3.txt7.夏氏偏差修正法：%首先将TXT中的数据装载在Matlab中：clear;clc;Wbs = 0;while Wbs<1disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt') filename = input('输入所需打开文件名:\','s');[Wbs,Message] = fopen(filename,'r');%打开文件；endData = fscanf(Wbs,'%g %g',[2 inf]);%生成数据矩阵Data = Data';%使用夏氏修正法，对2阶系统进行参数辨识n = 2;L = length(Data);N = L-n;U = Data(:,1);Y = Data(:,2);%首先使用最小二乘法，求系统参数向量初值glOL =[-Y(2:L-1),-Y(1:L-2),U(2:L-1),U(1:L-2)];Zgl1 = Data(3:L,2);Sgl1 = glOL'*glOL;Sgl2=inv(Sgl1);Sgl3=glOL'*Zgl1;Xsthita = Sgl2*Sgl3;%计算参数矩阵%得到初值后，以下使用夏氏修正法进行参数辨识：thitab0 = 0;%设偏差项的偏差初值为0Fa = Sgl2*glOL';M = eye(N)-glOL*Sgl2*glOL';F = eye(N);if(F>=10^-6*eye(N))E1 = Zgl1-glOL*Xsthita;%计算残差Eomiga(2:N,1) = -E1(1:N-1);omiga(3:N,2) = -E1(1:N-2);D = omiga'*M*omiga;Fx = inv(D)*omiga'*M*Zgl1;thitab = Fa*omiga*Fx;Xsthita = Xsthita - thitab;F = thitab - thitab0;thitab0 = thitab;enddisp('夏氏修正法')Xsa1 = Xsthita(1),Xsa2 = Xsthita(2),Xsb1 = Xsthita(3),Xsb2 = Xsthita(4)实验结果：Uy1.txt uy2.txt uy3.txt%首先将TXT中的数据装载在Matlab中：clear;clc;Wbs = 0;while Wbs<1disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt') filename = input('输入所需打开文件名:\','s');[Wbs,Message] = fopen(filename,'r');%打开文件；endData = fscanf(Wbs,'%g %g',[2 inf]);%生成数据矩阵Data = Data';%最小二乘法，取b0=0n = 2;%根据题目，本系统为2阶系统L = length(Data);%辨识所需数据的总长度N = L-n;%构造测量矩阵的数据长度FIA = zeros(N,2*n);%构造测量矩阵for i = 1:N %测量矩阵赋值for l = 1:n*2if(l>n)FIA(i,l) = Data(n+2+i-l,1);elseFIA(i,l) = -Data(n+i-l,2);endendendY = Data(n+1:n+N,2);%输出数据矩阵thita = inv(FIA'*FIA)*FIA'*Y;%计算参数矩阵thitaLS = thita;%最小二乘估值m = inv(FIA'*FIA)*FIA';for i = 1:1:1000%构造OmegaErr = Y-FIA*thitaLS;%计算残差Eomiga(2:N,1) = -Err(1:N-1);omiga(3:N,2) = -Err(1:N-2);F = inv(omiga'*omiga)*omiga'*Err;thita = thitaLS-m*omiga*F;enddisp('使用夏氏改良法辨识结果为：')Xga1 = thita(1),Xga2 = thita(2),Xgb1 = thita(3),Xgb2 = thita(4)实验结果：Uy1.txt uy2.txt uy3.txt%首先将TXT中的数据装载在Matlab中：clear;clc;Wbs = 0;while Wbs<1disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt')filename = input('输入所需打开文件名:\','s');[Wbs,Message] = fopen(filename,'r');%打开文件；endData = fscanf(Wbs,'%g %g',[2 inf]);%生成数据矩阵Data = Data';%这里选用450组数据进行递推迭代n = 2;N = 450-n;U = Data(:,1);Y = Data(:,2);%初始赋值如下：beta = [0 0 0 0 0 0]';P0= 10^10*eye(6,6);E1 = 0;E2 = 0;%循环迭代for i = 3:NE1 = Y(n+i)-beta(1)*Y(n+i-1)-beta(2)*Y(n+i-2)-beta(3)*Y(n+i-1)-beta(4)*U(n+i-2);E2 = Y(n+i-1)-beta(1)*Y(n+i-2)-beta(2)*Y(n+i-3)-beta(3)*U(n+i-2)-beta(4)*U(n+i-3); FIA = [-Y(n+i) -Y(n+i-1) U(n+i) U(n+i-1) -E1 -E2]';K = P0*FIA*inv(1+FIA'*P0*FIA);P1 = P0-K*FIA'*P0;P0 = P1;beta1 = beta+K*(Y(n+i+1)-FIA'*beta);beta = beta1;enddisp('取450组数据进行递推夏氏辨识的结果为:');Dxa1=beta(1),Dxa2=beta(2),Dxb1=beta(3),Dxb2=beta(4)实验结果：Uy1.txt uy2.txt uy3.txt10.增广矩阵法：%首先将TXT中的数据装载在Matlab中：clear;clc;Wbs = 0;while Wbs<1disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt') filename = input('输入所需打开文件名:\','s');[Wbs,Message] = fopen(filename,'r');%打开文件；endData = fscanf(Wbs,'%g %g',[2 inf]);%生成数据矩阵Data = Data';%使用增广矩阵法进行2阶系统的参数辨识n = 2;L = length(Data);U = Data(1:L,1);Y = Data(1:L,2);for i = 1:3*n %首先给辨识参数向量赋初值，这里取thita0为0thita0(i,1) = 0;endP0 = 10^10*eye(3*n,3*n);%然后设定初始状态P0为充分大的对角阵ERR0 = 0.00000000001;%收敛情况指标，相对误差ERR0thita = [thita0,zeros(3*n,L)];%设定参数向量的初值及大小zs = zeros(L,1);%构造初始噪声矩阵hn = [0,0,0,0,0,0]';for j = n+1:1:Lzs(j-1) = Y(j-1)-hn'*thita0;hn = [-Y(j-1),-Y(j-2),U(j-1),U(j-2),zs(j-1),zs(j-2)]';K = P0*hn*inv(hn'*P0*hn+1);%计算修正矩阵Kthita0 = thita0+K*[Y(j)-hn'*thita0];P0 = P0-K*hn'*P0;%计算状态矩阵P(N)thita(:,j) = thita0;E = abs((thita(:,j-1)-thita0)./thita0);if E<=ERR0 break;%根据设定的误差限，判断何时跳出endendZgma1=thita0(1),Zgma2=thita0(2),Zgmb1=thita0(3),Zgmb2=thita0(4), Zgmc1=thita0(5),Zgmc2=thita0(6)Uy1.txt uy2.txt uy3.txt11.极大似然法：%首先将TXT中的数据装载在Matlab中：clear;clc;Wbs = 0;while Wbs<1disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt')filename = input('输入所需打开文件名:\','s');[Wbs,Message] = fopen(filename,'r');%打开文件；endData = fscanf(Wbs,'%g %g',[2 inf]);%生成数据矩阵Data = Data';n = 2;L = length(Data);N = L-n;U = Data(:,1);Y = Data(:,2);mthita = zeros(3*n,1);%给被辨识参数矩阵赋初始值%用最小二乘法求初始的被辨识参数矩阵OL = [-Y(2:L-1),-Y(1:L-2),U(2:L-1),U(1:L-2)];Z1 = Data(3:L,2);thita = inv(OL'*OL)*OL'*Z1;mthita(1:2*n,1) = thita;%得到出之后，使用牛顿—拉布森极大似然法编程E = zeros(L,1);J = 0;sgma0 = 1;Ea1 = zeros(L,1);Ea2 = zeros(L,1);Eb1 = zeros(L,1);Eb2 = zeros(L,1);Ec1 = zeros(L,1);Ec2 = zeros(L,1);%给J的梯度和Hassian矩阵赋初始值Ethita = zeros(3*n,1);fd = zeros(3*n,1);Hassian = zeros(3*n,3*n);for p = 1:500Ja1 = mthita(1);Ja2 = mthita(2);Jb1 = mthita(3);Jb2 = mthita(4);Jc1 = mthita(5);Jc2 = mthita(6);for i = n+1:Ly2(i) = -Ja1*Y(i-1)-Ja2*Y(i-2)+Jb1*U(i-1)+Jb2*U(i-2)+Jc1*E(i-1)+Jc2*E(i-2);E(i) = Y(i)-y2(i);J = 0.5*(E(i)^2);sgma = (1/N)*(E(i)^2);endfor j = n+1:1:LEa1(j) = Y(j-1)-[Jc1,Jc2]*[Ea1(j-1),Ea1(j-2)]'; Ea2(j) = Y(j-2)-[Jc1,Jc2]*[Ea2(j-1),Ea2(j-2)]';Eb1(j) = -U(j-1)-[Jc1,Jc2]*[Eb1(j-1),Eb1(j-2)]';Eb2(j) = -U(j-2)-[Jc1,Jc2]*[Eb2(j-1),Eb2(j-2)]';Ec1(j) = -E(j-1)-[Jc1,Jc2]*[Ec1(j-1),Ec1(j-2)]';Ec2(j) = -E(j-2)-[Jc1,Jc2]*[Ec2(j-1),Ec2(j-2)]';Ethita = [Ea1(j),Ea2(j),Eb1(j),Eb2(j),Ec1(j),Ec2(j)];fd = fd+E(j)*Ethita';%迭代计算J的梯度矩阵Hassian = Hassian + Ethita'*Ethita;%迭代计算Hassian矩阵end%用牛顿-拉卜森法计算被辨识的参数的估值mthita = mthita - inv(Hassian)*fd;JFIA = (sgma-sgma0)/sgma0;sgma0 = sgma;%跳出迭代的判断if JFIA<0.000000001 breakendenddisp('极大似然法（牛顿-拉卜森方法），得到的结果为：')Ja1 = mthita(1),Ja2 = mthita(2),Jb1 = mthita(3),Jb2 = mthita(4),Jc1 = mthita(5),Jc2 = mthita(6)实验结果：Uy1.txt uy2.txt uy3.txt12.分析噪声)(k 特性:RLS 法：当噪声比较小时，辨识精度较高；当噪声比较大时，收敛点可能不唯一，参数估计值往往是有偏的；但是运用此方法时，数据要充分多，否则辨识精度明显下降；另外噪声模型阶次不宜过高，也会影响精度，所以误差的原因主要是数据量的大小以及噪声模型阶次。